Archive

Posts Tagged ‘PDE’

On the first eigenfunction of the Dirichlet-Laplacian on the regular polygon

October 7, 2021 Leave a comment

Suppose that {P_n} is the regular polygon inscribed in the unit circle with one vertex at {(1,0)}. Consider the Dirichlet-Laplace eigenvalue problem

\displaystyle \left\{ \begin{array}{rcll} -\Delta u &= &\lambda u& \text{ in }P_n \\ u & = & 0 & \text{ on } \partial P_n \end{array}\right.

This problem has solutions for a sequence of eigenvalues

\displaystyle 0<\lambda_1(P_n)<\lambda_2(P_n)\leq ... \rightarrow +\infty.

Denote by {u_i} the eigenfunction associated to {\lambda_i}, {i \geq 1} normalized such that {\int_{P_n} u_i^2 = 1}. The case of the first eigenvalue is particular. Since the first eigenvalue is simple, there exists a unique normalized eigenfunction such that {u_1 \geq 0}. This implies that the function {u_1} has all the symmetries of the regular polygon. Below you can see an image of the eigenfunction {u_1} in the case of the regular pentagon.

Consider now the triangles {T_j} given by vertices {(0,0),(\cos (j\theta),\sin (j\theta)),(\cos((j+1)\theta),\sin((j+1)\theta)} for {j=0,..,n-1}, where {\theta = 2\pi/n}. These triangles decompose {P_n} into equal slices.

Furthermore, let {A_{xx} = \int_{T_0} (\partial_x u)^2, A_{yy} = \int_{T_0} (\partial_y u)^2, A_{xy} = \int_{T_0} \partial_x u \partial_y u}. Then we have, by the symmetry of the eigenfunction, that {A_{xx}+A_{yy} = \lambda_1/n}.

In the following, due to the symmetry we can infer that {\partial_x u} is an even function with respect to {y} and {\partial_y u} is odd. The fact that the gradients undergo a rotation when transferred from {T_{-1}} to {T_0} we have the matrix equality

\displaystyle \begin{pmatrix} \cos\theta& -\sin \theta \\ \sin \theta& \cos \theta \end{pmatrix}\begin{pmatrix} A_{xx} & -A_{xy} \\ -A_{xy} & A_{yy} \end{pmatrix} \begin{pmatrix} \cos\theta& \sin \theta \\ -\sin \theta& \cos \theta \end{pmatrix} = \begin{pmatrix} A_{xx} & A_{xy} \\ A_{xy} & A_{yy} \end{pmatrix}which implies that {-A_{xx}\sin \frac{2\pi}{n} + A_{yy} \sin \frac{2\pi}{n} + 2A_{xy} \cos \frac{2\pi}{n} = 0}.

Up until now we have two relations between {A_{xx},A_{yy},A_{xy}}. Any supplementary independent relation would completely determine these quantities. 

Question: Can you find such a relation?

FreeFem++ Tutorial – Part 2

May 1, 2018 2 comments

Click here for the first part. Some other posts related to FreeFem: link 1, link 2.

Here are a few tricks, once you know the basics of FreeFem. If your plan is straightforward: define the domain, build the mesh, define the problem, solve the problem, plot the result… then things are rather easy. If you want to do stuff in a loop, like an optimization problem, things may get more complicated, but FreeFem still has lots of tricks up its sleeves. I’ll go through some of them.

  1. Defining the geometry of the domain using border might be tricky at the beginning. Keep in mind that the domain should be on the left side of the curves definining the boundaries. This could be acheived in the parametrization chosen or you could reverse the parametrization of a particular part of the domain in the following way. Suppose you have the pieces of boundary C1,C2,C3,C4. You wish to build a mesh with the command mesh Th = buildmesh(C1(100)+C2(100)+C3(100)+C4(100));but you get an error concerning a bad sens on one of the boundaries. If you identify, for example, that you need to change the orientation of C3, you can acheive this by changing the sign of the integer defining the number of points on that part of the boundary: mesh Th = buildmesh(C1(100)+C2(100)+C3(-100)+C4(100));You could make sure that you have the right orientations and good connectivities for all the boundaries by running a plot command before trying mesh: something likeplot(C1(100)+C2(100)+C3(-100)+C4(100));should produce a graph of all your boundaries with arrows showing the orientations. This is good as a debug tool when you don’t know where the error comes from when you define the domain.
  2. Keep in mind that there are also other ways to build a mesh, like square, which meshes a rectangular domain with quite a few options and trunc which truncates or modifies a mesh following some criteria. Take a look in the documentation to see all the options.
  3. Let’s say that you need to build a complex boundary, but with many components with similar properties. Two examples come to mind: polygons and domains with many circular holes. In order to do this keep in mind that it is possible to define a some kind of “vectorial boundary”. Let’s say that you want to mesh a polygon and you have the coordinates stored in the arrays xs,ys. Furthermore, you have another array of integers ind which point to the index of the next vertex. Then the boundary of the polygon could be defined with the following syntax:
    border poly(t=0,1; i){
    x=(1-t)*xx[i]+t*xx[ind(i)];
    y=(1-t)*yy[i]+t*yy[ind(i)];
    label=i;
    }
    Now the mesh could be constructed using a vector of integers NC containing the number of desired points on each of the sides of the polygon:

    mesh Th = buildmesh (poly(NC));
  4. You can change a 2D mesh using the command adaptmesh. There are various options and you’ll need to search in the documentation for a complete list. I use it in order to improve a mesh build with buildmesh. An example of command which gives good results in some of my codes is:

    Th = adaptmesh(Th,0.02,IsMetric=1,nbvx=30000);

    The parameters are related to the size of the triangles, the geometric properties of the griangles and the maximal number of vertices you want in your mesh. Experimenting a bit might give you a better idea of how this command works in practice.
  5. There are two ways of defining the problems in FreeFem. One is with solve and the other one is with problem. If you use solve then FreeFem solves the problem where it is defined. If you use problem then FreeFem remembers the problem as a variable and will solve it whenever you call this variable. This is useful when solving the same problem multiple times. Note that you can modify the coefficients of the PDE and FreeFem will build the new problem with the updated coefficients.
  6. In a following post I’ll talk about simplifying your code using macros and functions. What is good to keep in mind is that macros are verbatim code replacements, which are quite useful when dealing with complex formulas in your problem definition or elsewhere. Functions allow you to run a part of the code with various parameters (just like functions in other languages like Matlab).

I’ll finish with a code which computes the eigenvalue of the Laplace operator on a square domain with multiple holes. Try and figure out what the commands and parameters do.

int N=5; //number of holes
int k=1; //number of eigenvalue

int nb = 100;  // parameter for mesh size
verbosity = 10; // parameter for the infos FreeFem gives back

real delta = 0.1;
int nbd = floor(1.0*nb/N*delta*2*pi);

int bsquare = 0;
int bdisk   = 1;

// vertices of the squares
real[int] xs(N^2),ys(N^2);

real[int] initx = 0:(N^2-1);
for(int i=0;i&amp;amp;amp;lt;N^2;i++){
  xs[i] = floor(i/N)+0.5;
  ys[i] = (i%N)+0.5;
}


int[int] Nd(N^2);

Nd = 1;
Nd = -nbd*Nd;


// sides of the square
border Cd(t=0,N){x = t;y=0;label=bsquare;}
border Cr(t=0,N){x = N;y=t;label=bsquare;}
border Cu(t=0,N){x = N-t;y=N;label=bsquare;}
border Cl(t=0,N){x = 0;y=N-t;label=bsquare;}

border disks(t=0,2*pi; i){
    x=xs[i]+delta*cos(t);
    y=ys[i]+delta*sin(t);
    label=bdisk;
}

plot(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));

mesh Th = buildmesh(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));

plot(Th);

int[int] bc = [0,1]; // Dirichlet boundary conditions

load "Element_P3";
fespace Vh(Th,P2);
// variables on the mesh
Vh u1,u2;

// Define the problem in weak form
varf a(u1,u2) = int2d(Th)
(dx(u1)*dx(u2) + dy(u1)*dy(u2))+on(1,u1=0)//on(C1,C2,C3,C4,u1=1)
                               +on(bc,u1=0);
varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ;
// define matrices for the eigenvalue problem
matrix A= a(Vh,Vh,solver=Crout,factorize=1);
matrix B= b(Vh,Vh,solver=CG,eps=1e-20); 

// we are interested only in the first eigenvalue
int eigCount = k;
real[int] ev(eigCount); // Holds eigenvalues
Vh[int] eV(eigCount);   // holds eigenfunctions
// Solve Ax=lBx
int numEigs = EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eV);

plot(eV[k-1],fill=1,nbiso=50,value=1);

Here is the mesh and the solution given by the above program:

Mesh

Eig5

FreeFem++ Tutorial – Part 1

February 24, 2016 2 comments

First of all, FreeFem is a numerical computing software which allows a fast and automatized treatment of a variety of problems related to partial differential equations. Its name, FreeFem, speaks for itself: it is free and it uses the finite element method. Here are a few reasons for which you may choose to use FreeFem for a certain task:

  1. It allows the user to easily define 2D (and 3D) geometries and it does all the work regarding the construction of meshes on these domains.
  2. The problems you want to solve can be easily written in the program once we know their weak forms.
  3. Once we have variables defined on meshes or solutions to some PDE, we can easily compute all sorts of quantities like integral energies, etc.

Before showing a first example, you need to install FreeFem. If you are not familiar with command line work or you just want to get to work, like me, you can install the visual version of FreeFem which is available here. Of course, you can find example programs in the FreeFem manual or by making a brief search on the internet.

I’ll present some basic stuff, which will allow us in the end to solve the Laplace equation in a circular domain. Once we have the structure of the program, it is possible to change the shape of the domain in no time.

Read more…

Simple triangle mesh – Matlab code

April 20, 2015 Leave a comment

There are cases when you need to define a simple, regular mesh on a nice set, like a triangle, or a polygon. There are software which are capable of creating such meshes, but most of them will require some computation time (Delaunay triangulation, plus optimization, edge swapping, etc.). I propose to you a simple algorithm for defining a mesh on a triangle. Of course, this can be extended to polygons or other cases (I used a variant of it to create a mesh of a spherical surface).

The idea is very simple: start with three points {P_1,P_2,P_3}, and define the (only) initial triangle as {T = [1,\ 2,\ 3]}. Then, at each step, take a triangle {t_i} from the list {T}, add the midpoints of {t_i} to the list of points and replace the triangle {t_i} with the four smaller triangles determined by the vertices and midpoints of {t_i}. This is a basic mesh refining strategy. Thus, in short, we start with a triangle and do a number of mesh refinements, until we reach the desired side-length. Of course, at each step we should be careful not to add the midpoint of an edge two times (once for every neighbouring triangle), but in what follows, I will always remove duplicates by bruteforce in the end (not the most economic solution).

Read more…

Categories: Uncategorized Tags: , ,

Regularity of the weak solution Part 1

October 27, 2012 Leave a comment

I will present here how to recover the regularity of a weak solution for the Dirichlet problem. The arguments can easily be adapted to most of the weak formulations involving the Laplace operator, the essential tool being the estimate of the {L^2} norm of the derivatives of the solution {u}. The arguments are adapted from H. Brezis, Functional Analysis, Sobolev Spaces and Partial Differential Equations, Chapter 9, and the tool named the method of translations is due to L. Niremberg.

Consider the problem

\displaystyle \begin{cases} -\Delta u= f & \text{ in }\Omega \\ \hfill u=0 & \text{ on }\partial \Omega. \end{cases}

whose weak variational form is

\displaystyle \int_\Omega \nabla u \nabla v = \int_\Omega fv,\ \forall v \in H_0^1(\Omega).

Note that in the variational form we only assume {u \in H^1(\Omega)}, and to be able to recover the PDE we need to use Green’s formula, which is valid if {u \in H^2(\Omega)}.

Read more…

Weak formulation for Laplace Equation with Robin boundary conditions

October 22, 2012 15 comments

Consider {\Omega \subset \Bbb{R}^N} an open set with Lipschitz boundary and consider on {\Omega} the following problem

\displaystyle \begin{cases} -\Delta u =f &\text{ in }\Omega \\ \frac{\partial u}{\partial n}+\beta u =0 & \text{ on }\partial \Omega. \end{cases}

where {\beta>0} is a constant. This is the Laplace equation with Robin boundary conditions. I will prove that the problem is well posed and for each {f \in L^2(\Omega)} there exists a solution {u \in H^2(\Omega)}.

Read more…

Finite Difference Method for 2D Laplace equation

October 19, 2012 21 comments

[Edit: This is, in fact Poisson’s equation.]

[For solving this equation on an arbitrary region using the finite difference method, take a look at this post.]

I will present here how to solve the Laplace equation using finite differences 2-dimensional case:

\displaystyle \begin{cases} - \Delta u=1 & \text{ in }\Omega=[0,1]^2 \\ u=0 &\text{ on }\partial \Omega. \end{cases}

Read more…

Lax Milgram application

October 14, 2012 Leave a comment

Let {I=(0,2)} and {V=H^1(I)}. Consider the bilinear form

\displaystyle a(u,v)= \int_0^2 u'(t)v'(t)dt +\left( \int_0^1 u(t)dt\right)\left( \int_0^1 v(t)dt\right).

1. Check that {a(u,v)} is a continuous symmetric bilinear form and that {a(u,u)=0} implies {u=0}.

2. Prove that {a} is coercive.

3. Deduce that for every {f \in L^2(I)} there exists a unique {u \in H^1(I)} satisfying

\displaystyle (1) \ \ \ \ a(u,v)=\int_0^2 fv,\ \forall v \in H^1(I).

What is the corresponding minimization problem?

4. Show that the solution of {(1)} belongs to {H^2(I)} (and in particular {u \in C^1(\overline{I})}). Determine the equation and the boundary conditions satisfied by {u}.

5. Assume that {f \in C(\overline{I})}, and let {u} be the solution of {(1)}. Prove that {u} belongs to {W^{2,p}(I)} for every {p<\infty}. Show that {u \in C^2(\overline{I})} if and only if {\int_If=0}.

6. Determine explicitly the solution {u} of {(1)} when {f} is a constant.

7. Set {u=Tf}, where {u} is the solution of {(1)} and {f \in L^2(I)}. Check that {T} is a self-adjoint compact operator from {L^2(I)} into itself.

8. Study the eigenvalues of {T}.

H. Brezis, Functional Analysis

Read more…

One dimensional trace of a Sobolev function

October 7, 2012 Leave a comment

Check that the mapping {u\mapsto u(0)} from {H^1(\Bbb{R})} to {\Bbb{R}} is a continuous linear functional on {H^1(0,1)}. Deduce that there exists a unique {v_0 \in H^1(0,1)} such that

\displaystyle u(0)=\int_0^1(u'v_0'+uv_0).

Show that {v_0} is the solution of some differential equation with appropriate boundary conditions and compute {v_0} explicitly.

H. Brezis, Functional Analysis, Ex 8.18

Read more…

Cinetic Equations Problem

September 13, 2012 Leave a comment

Cinetic Models Application

Consider the following BGK type model

 

\displaystyle (1) \ \ \ \ \begin{cases} \partial_tf+v\cdot \nabla_x f+f=\chi_{n(x,t)}(v), & x,v \in \Bbb{R}\\ f(x,v,0)=f_0(x,v) \geq 0 & f_0 \in L^1(\Bbb{R}^2) \end{cases}

where

\displaystyle n(x,t)=\int_\Bbb{R} f(x,v,t)dv

and

\displaystyle \chi_n(v)= \begin{cases} 1 & \text{if } 0 \leq v\leq n(x,t), \\ 0 & \text{otherwise} \end{cases}

I – Let {f} be the solution of {(1)} with {\chi_{n(x,t)}(v) \in L^1_{t,x,v}}.

1. Prove that {f\geq 0}.

2. Prove that {\|f(x,v,t)\|_{L^1(\Bbb{R}^2)}=\|f_0\|_{L^1(\Bbb{R}^2)}}.

3. Prove that if {f_0(x,v) \leq 1} then {f(x,v,t)\leq 1}.

II – Let {\mathcal{T}} be the operator defined on {E=\{\rho(x,t) \geq 0,\ \rho \in L^\infty([0,T];L^1(\Bbb{R}))\}} which assigns to {\rho} the function {n(x,t)=\displaystyle \int_\Bbb{R} f(x,v,t)dv} where {f} is the solution of

\displaystyle (2) \ \ \ \ \begin{cases} \partial_tf+v\cdot \nabla_x f+f=\chi_{\rho(x,t)}(v), & \\ f(x,v,0)=f_0(x,v) \geq 0 & f_0 \in L^1(\Bbb{R}^2) \end{cases}

4. Let {C=\displaystyle\left\{ \rho \in E : \int_\Bbb{R} \rho(x,t)dx \leq \|f_0\|_{L^1}\right\}}. Prove that {\mathcal{T}} maps {C} to {C} and that {C} is a convex closed subset of {E}.

5. Prove that for every {\rho_1,\rho_2} in {C} we have

\displaystyle \|\chi_{\rho_1}-\chi_{\rho_2}\|_{L^1(\Bbb{R}^2)}=\|\rho_1-\rho_2\|_{L^1(\Bbb{R})}

6. Prove that if we denote by {f_1,f_2} the solutions of {(2)} corresponding to {\rho_1,\rho_2} we have for every {t \in [0,T]} that

\displaystyle \|f(x,v,t)-f_2(x,v,t)\|_{L^1(\Bbb{R}^2)} \leq (1-e^{-t})\sup_{t \in [0,T]} \|\rho_1-\rho_2\|_{L^1(\Bbb{R}_x)}

7. Prove that {\mathcal{T}} is a contraction from {C} to {C} and that {(1)} has a unique solution.

III Choose {f_0=\chi_{R_\varepsilon(x)}} bounded in {L^1(\Bbb{R})\cap L^\infty(\Bbb{R})}, and denote {f_\varepsilon} the solution of {(1)} which corresponds to {f_0}.

8. Prove that if {v \geq \sup_\varepsilon \|R_\varepsilon\|_{L^\infty}=r}, then {f_\varepsilon(x,v,t)=0,\ \forall x \in \Bbb{R}, t \geq 0}. (we can prove that {f_\varepsilon \leq \chi_R(v)}).

9. Prove that {f_\varepsilon} and {\chi_{n_\varepsilon}} are bounded in {L^\infty([0,T];L^2(\Bbb{R}^2))}.

10. Prove that {(n_\varepsilon)_{\varepsilon \geq 0}} is compact in {L^2([0,T]\times B_K)} for every ball {B_K} of {\Bbb{R}_x}.

11. Prove that, up to a subsequence, we ca pass to the limit in the sense of distributions in {(1)} with {f_0=\chi_{R_\varepsilon(x)}} as {\varepsilon \rightarrow 0}.

Uniqueness and Error estimates via Kinetic Entropy Defect Measure

March 7, 2012 Leave a comment

Here are a few thoughts from my preparation for the exam of Kinetic Equations at Universite de Savoie, France. The teachers of the course were Christian Bourdarias and Stephane Gerbi. I had to study an article of Benoit Perthame entitled Uniqueness and Error estimates in First Order Quasilinear Conservation Laws via the Kinetic Entropy Defect Measure.

This was a very nice article to study, since it used many things like distribution theory, measures and regularization. It showed the power of these tools, and motivated me to learn more about them.

As the title of the article says, a relatively new proof of the uniqueness of the solution for a scalar conservation law coupled with some entropy inequalities is given. The only known proof at the time the article was published was due to Kruzkov and was more intricate and difficult to understand than the one provided in the article. The estimates on the entropy defect measure, which will be introduced can yield some error term approximation for approximate equation, which in particular imply unicity at once.

Here are my detailed notes on the article. They are handwritten, but I think they are readable. Perthame-Uniqueness and Error Estimates

Read more…

Traian Lalescu Student Contest 2011 Problem 4

May 25, 2011 Leave a comment

Let D=(0,\infty)\times (0,\infty),\ u \in C^1(D) and \varepsilon>0 fixed.

1) Prove that x \frac{\partial u}{\partial x}(x,y)+y\frac{\partial u}{\partial y}(x,y)=u(x,y),\ \forall (x,y) \in D if and only if there exists \phi \in C^1(0,\infty) such that u(x,y)=x\phi(y/x),\ \forall (x,y) \in D.

2) Prove that \left|x \frac{\partial u}{\partial x}(x,y)+y\frac{\partial u}{\partial y}(x,y)-u(x,y)\right|\leq \varepsilon,\ \forall (x,y) \in D then there is a unique function \phi \in C^1(0,\infty) such that \left|u(x,y)-x\phi(y/x),\right| \leq \varepsilon \ \forall (x,y) \in D.

Weakly harmonic implies harmonic

Suppose u \in C^0(\Omega) is weakly harmonic on an open set \Omega, i.e. the relation \displaystyle \int_\Omega u \Delta \phi dx =0 holds for all \phi \in C_0^2(\Omega). Show that u is harmonic in \Omega.

PHD Iowa (6202)

See this blog post for a proof.

Divergent integral

November 23, 2010 Leave a comment

Let u satisfy u \in \mathcal{C}^2(\Bbb{R}^n),\ \Delta u=0 on \Bbb{R}^n. Show that the integral \displaystyle \int_{\Bbb{R}^n}u^2 dx is convergent if and only if u \equiv 0.
PHD 6201

Read more…

Compacts convergence implies characteristic convergence

November 22, 2010 Leave a comment

Prove that if the sequence of of open sets (\Omega_n) converges to \Omega in the compact way, and if the boundary of \Omega is of measure zero, then \Omega_n converges in the way of characteristic functions.

Read more…

Categories: Uncategorized Tags: ,

Frontier convergence

November 21, 2010 3 comments

Suppose that a sequence of open sets \Omega_n \subset \Bbb{R}^N converges in the way of compact sets to an open set \Omega. Is it true that \forall x \in \partial \Omega we nave d(x,\partial \Omega_n) \to 0?

Same question for the convergence in the way of characteristic functions and for the Hausdorff convergence.

Read more…

Categories: Uncategorized Tags: ,

Open sets convergence

November 21, 2010 3 comments

Consider the sequence of open sets of [0,1] defined by
a) \displaystyle \Omega_n:= \bigcup_{k=0}^{2^{n-1}-1} \left(\frac{2k}{2^n},\frac{2k+1}{2^n}\right).

b) \displaystyle \Omega_n:= \bigcup_{k=0}^{2^{n}-1} \left(\frac{k}{2^n},\frac{k+1}{2^n}\right).

Study the convergence of the sets \Omega_n in the following types of convergence:

  • weak convergence of characteristic functions;
  • convergence in the topology defined by the Hausdorff distance;
  • convergence in the way of compact sets.

(Each of the previous types of convergence will be presented below.)

Read more…

Harmonic function

April 30, 2010 Leave a comment

Let e_1,e_2,...,e_n be semilines on the plane starting from a common point. Prove that if there doesn’t exist any function u \not\equiv 0 harmonic on the whole plane that vanishes on the set e_1\cup e_2\cup...\cup e_n, then there exists a pair (i,j) such that there is no function u\not\equiv 0 harmonic on the whole plane such that u vanishes on e_i \cup e_j.
Miklos Schweitzer 2001
Read more…