Archive

Archive for the ‘Analysis’ Category

FreeFem++ Tutorial – Part 2

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 like

    plot(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<N^2;i++){
  xs[i] = floor(i/N)+0.5;
  ys[i] = (i%N)+0.5;
}

cout << xs << endl;
cout << ys << endl;

int[int] Nd(N^2);

Nd = 1;
Nd = -nbd*Nd;

cout << Nd << endl;

// 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);

cout &lt;&lt; ev[k-1] &lt;<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;"></span>&lt; endl;

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

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

MeshEig5

Advertisements

SEEMOUS 2018 – Problems

March 1, 2018 Leave a comment

Problem 1. Let {f:[0,1] \rightarrow (0,1)} be a Riemann integrable function. Show that

\displaystyle \frac{\displaystyle 2\int_0^1 xf^2(x) dx }{\displaystyle \int_0^1 (f^2(x)+1)dx }< \frac{\displaystyle \int_0^1 f^2(x) dx}{\displaystyle \int_0^1 f(x)dx}.

Problem 2. Let {m,n,p,q \geq 1} and let the matrices {A \in \mathcal M_{m,n}(\Bbb{R})}, {B \in \mathcal M_{n,p}(\Bbb{R})}, {C \in \mathcal M_{p,q}(\Bbb{R})}, {D \in \mathcal M_{q,m}(\Bbb{R})} be such that

\displaystyle A^t = BCD,\ B^t = CDA,\ C^t = DAB,\ D^t = ABC.

Prove that {(ABCD)^2 = ABCD}.

Problem 3. Let {A,B \in \mathcal M_{2018}(\Bbb{R})} such that {AB = BA} and {A^{2018} = B^{2018} = I}, where {I} is the identity matrix. Prove that if {\text{tr}(AB) = 2018} then {\text{tr}(A) = \text{tr}(B)}.

Problem 4. (a) Let {f: \Bbb{R} \rightarrow \Bbb{R}} be a polynomial function. Prove that

\displaystyle \int_0^\infty e^{-x} f(x) dx = f(0)+f'(0)+f''(0)+...

(b) Let {f} be a function which has a Taylor series expansion at {0} with radius of convergence {R=\infty}. Prove that if {\displaystyle \sum_{n=0}^\infty f^{(n)}(0)} converges absolutely then {\displaystyle \int_0^{\infty} e^{-x} f(x)dx} converges and

\displaystyle \sum_{n=0}^\infty f^{(n)}(0) = \int_0^\infty e^{-x} f(x).

Source: official site of SEEMOUS 2018 

Hints: 1. Just use 2f(x) \leq f^2(x)+1  and xf^2(x) < f^2(x). The strict inequality comes from the fact that the Riemann integral of strictly positive function cannot be equal to zero. This problem was too simple…

2. Use the fact that ABCD = AA^t, therefore ABCD is symmetric and positive definite. Next, notice that (ABCD)^3 = ABCDABCDABCD = D^tC^tB^tA^t = (ABCD)^t = ABCD. Notice that ABCD  is diagonalizable and has eigenvalues among -1,0,1. Since it is also positive definite, -1 cannot be an eigenvalue. This allows to conclude.

3. First note that the commutativity allows us to diagonalize A,B  using the same basis. Next, note that A,B both have eigenvalues of modulus one. Then the trace of AB is simply the sum \sum \lambda_i\mu_i where \lambda_i,\mu_i are eigenvalues of A and B, respectively. The fact that the trace equals 2018  and the triangle inequality shows that eigenvalues of A are a multiple of eigenvalues of B. Finish by observing that they have the same eigenvalues.

4. (a) Integrate by parts and use a recurrence. (b) Use (a) and the approximation of a continuous function by polynomials on compacts to conclude.

I’m not sure about what others think, but the problems of this year seemed a bit too straightforward.

Putnam 2017 A3 – Solution

December 4, 2017 Leave a comment

Problem A3. Denote {\phi = f-g}. Then {\phi} is continuous and {\int_a^b \phi = 0}. We can see that

\displaystyle I_{n+1}-I_n = \int_a^b (f/g)^n \phi = \int_{\phi\geq 0} (f/g)^n \phi+ \int_{\phi<0} (f/g)^n \phi

Now note that on {\{ \phi>=0\}} we have {f/g \geq 1} so {(f/g)^n \phi \geq \phi}. Furthermore, on {\{\phi<0\}} we have {(f/g)^n <1} so multiplying with {\phi<0} we get {(f/g)^n \phi \geq \phi}. Therefore

\displaystyle I_{n+1}-I_n \geq \int_{\phi \geq 0} \phi + \int_{\phi<0} \phi = \int \phi = 0.

To prove that {I_n} goes to {+\infty} we can still work with {I_{n+1}-I_n}. Note that the negative part cannot get too big:

\displaystyle \left|\int_{ \phi <0 } (f/g)^n \phi \right| \leq \int_{\phi<0} |\phi| \leq \int_a^b |f-g|.

As for the positive part, taking {0<\varepsilon< \max_{[a,b]} \phi} we have

\displaystyle \int_{\phi\geq 0} (f/g)^n \phi \geq \int_{\phi>\varepsilon}(f/g)^n \varepsilon.

Next, note that on {\{ \phi \geq \varepsilon\}}

\displaystyle \frac{f}{g} = \frac{g+\phi}{g} \geq \frac{g+ \varepsilon}{g}.

We would need that the last term be larger than {1+\delta}. This is equivalent to {g\delta <\varepsilon}. Since {g} is continuous on {[a,b]}, it is bounded above, so some delta small enough exists in order for this to work.

Grouping all of the above we get that

\displaystyle I_{n+1}-I_n \geq \int_{\phi \geq 0} (f/g)^n \phi \geq \int_{\phi>\varepsilon} \varepsilon (1+\delta)^n.

Since {|\phi>\varepsilon|>0} we get that {I_{n+1}-I_n} goes to {+\infty}.

Project Euler Problem 285

March 25, 2017 Leave a comment

Another quite nice problem from Project Euler is number 285. The result of the problem depends on the computation of a certain probability, which in turn is related to the computation of a certain area. Below is an illustration of the computation for k equal to 10.

pb285_10

To save you some time, here’s a picture of the case k=1 which I ignored and spent quite a lot of time debugging… Plus, it only affects the last three digits or so after the decimal point…

pb285_1

Here’s a Matlab code which can construct the pictures above and can compute the result for low cases. To solve the problem, you should compute explicitly all these areas.


function problem285(k)

N = 100000;

a = rand(1,N);
b = rand(1,N);

ind = find(abs(sqrt((k*a+1).^2+(k*b+1).^2)-k)<0.5);

plot(a(ind),b(ind),'.');
axis equal

M = k;
pl = 1;

for k=1:M
if mod(k,100)==0
k
end
r1 = (k+0.5)/k;
r2 = (k-0.5)/k;

f1 = @(x) (x<=(-1/k+r1)).*(x>=(-1/k-r1)).*(sqrt(r1^2-(x+1/k).^2)-1/k).*(x>=0).*(x<=1); f1 = @(x) f1(x).*(f1(x)>=0);
f2 = @(x) (x<=(-1/k+r2)).*(x>=(-1/k-r2)).*(sqrt(r2^2-(x+1/k).^2)-1/k).*(x>=0).*(x<=1); f2 = @(x) f2(x).*(f2(x)>=0);

if k == pl
thetas = linspace(0,pi/2,200);
hold on
plot(-1/k+r1*cos(thetas),-1/k+r1*sin(thetas),'r','LineWidth',2);
plot(-1/k+r2*cos(thetas),-1/k+r2*sin(thetas),'r','LineWidth',2);
plot([0 1 1 0 0],[0 0 1 1 0],'k','LineWidth',3);
hold off
axis off
end

A(k) = integral(@(x) f1(x)-f2(x),0,1);

end

xs = xlim;
ys = ylim;

w = 0.01;
axis([xs(1)-w xs(2)+w ys(1)-w ys(2)+w]);

sum((1:k).*A)

Some of the easy Putnam 2016 Problems

December 11, 2016 Leave a comment

Here are a few of the problems of the Putnam 2016 contest. I choose to only list problems which I managed to solve. Most of them are pretty straightforward, so maybe the solutions posted here may be very similar to the official ones. You can find a complete list of the problems on other sites, for example here.

A1. Find the smallest integer {j} such that for every polynomial {p} with integer coefficients and every integer {k}, the number

\displaystyle p^{(j)}(k),

that is the {j}-th derivative of {p} evaluated at {k}, is divisible by {2016}.

Hints. Successive derivatives give rise to terms containing products of consecutive numbers. The product of {j} consecutive numbers is divisible by {j!}. Find the smallest number such that {2016 | j!}. Prove that {j-1} does not work by choosing {p = x^{j-1}}. Prove that {j} works by working only on monomials…

Read more…

IMC 2016 – Day 2 – Problem 7

July 28, 2016 Leave a comment

Problem 7. Today, Ivan the Confessor prefers continuous functions {f:[0,1]\rightarrow \Bbb{R}} satisfying {f(x)+f(y) \geq |x-y|} for all {x,y \in [0,1]}. Fin the minimum of {\int_0^1 f} over all preferred functions.

Read more…

IMC 2016 – Day 2 – Problem 6

July 28, 2016 1 comment

Problem 6. Let {(x_1,x_2,...)} be a sequence of positive real numbers satisfying {\displaystyle \sum_{n=1}^\infty \frac{x_n}{2n-1}=1}. Prove that

\displaystyle \sum_{k=1}^\infty \sum_{n=1}^k \frac{x_n}{k^2} \leq 2.

Read more…

%d bloggers like this: