Example of optimization under constraints – using Lagrange multipliers

Consider the functional defined by

and the set defined by .

1. Prove that admits a minimizer on and find it.
2. Let be the set . Does admit minimizers on ? If yes determine all points in which the maximum is attained.

Proof: This is a problem which illustrates well concepts related to optimality conditions for constrained problems. It is possible to solve these problems by two methods: one using only basic ideas and a second one using results from optimization theory. One can also see by examining the basic method how to recover naturally the optimality condition given by the Lagrange multipliers.

Advertisements

Optimizing a quadratic functional under affine constraints

Consider the quadratic functional defined by

for a real matrix which is symmetric positive definite and a vector . Consider also the application defined by where is a real matrix and . Furthermore suppose that and that is of maximal rank. The objective is to study the following optimization problem under equality constraints:

where is defined by .

1. Show that has a unique minimizer on .
2. Prove that if is the minimizer of on then there exists a unique vector such that solves the linear system

   

3. Show that the matrix of the above system is invertible.

IMO 2018 Problems – Day 2

Problem 4. A site is any point in the plane such that and are both positive integers less than or equal to 20.

Initially, each of the 400 sites is unoccupied. Amy and Ben take turns placing stones with Amy going first. On her turn, Amy places a new red stone on an unoccupied site such that the distance between any two sites occupied by red stones is not equal to . On his turn, Ben places a new blue stone on any unoccupied site. (A site occupied by a blue stone is allowed to be at any distance from any other occupied site.) They stop as soon as a player cannot place a stone.

Find the greatest such that Amy can ensure that she places at least red stones, no matter how Ben places his blue stones.

Problem 5. Let be an infinite sequence of positive integers. Suppose that there is an integer such that, for each , the number

is an integer. Prove that there is a positive integer such that for all .

Problem 6. A convex quadrilateral satisfies . Point lies inside so that and . Prove that .

Source: AoPS

IMO 2018 Problems – Day 1

Problem 1. Let be the circumcircle of acute triangle . Points and are on segments and respectively such that . The perpendicular bisectors of and intersect minor arcs and of at points and respectively. Prove that lines and are either parallel or they are the same line.

Problem 2. Find all integers for which there exist real numbers satisfying , and

For .

Problem 3. An anti-Pascal triangle is an equilateral triangular array of numbers such that, except for the numbers in the bottom row, each number is the absolute value of the difference of the two numbers immediately below it. For example, the following is an anti-Pascal triangle with four rows which contains every integer from to

Does there exist an anti-Pascal triangle with rows which contains every integer from to ?

Source: AoPS.

Balkan Mathematical Olympiad 2018

Problem 1. A quadrilateral is inscribed in a circle , where and is not parallel to . Point is the intersection of the diagonals and and the perpendicular from to intersects the segment at the point . If bisects the angle , prove that is a diameter of the circle .

Problem 2. Let be a positive rational number. Two ants are initially at the same point in the plane. In the -th minute each of them chooses whether to walk due north, east, south or west and then walks the distance of meters. After a whole number of minutes, they are at the same point in the plane (non necessarily ), but have not taken exactly the same route within that time. Determine all the possible values of .

Problem 3. Alice and Bob play the following game: They start with two non-empty piles of coins. Taking turns, with Alice playing first, each player chooses a pile with an even number of coins and moves half of the coins of this pile to the other pile. The came ends if a player cannot move, in which case the other player wins.

Determine all pairs of positive integers such that if initially the two piles have and coins, respectively, then Bob has a winning strategy.

Problem 4. Find all primes and such that divides .

Source: https://bmo2018.dms.rs/wp-content/uploads/2018/05/BMOproblems2018_English.pdf

Sum of the Euler Totient function

Given a positive integer , the Euler totient function is defined as the number of positive integers less than which are co-prime with (i.e. they have no common factors with ). There are formulas for computing starting from the factorization of . One such formula is

where the product is made over all primes dividing .

If you have to compute for all numbers less than a threshold then another property could be useful: is arithmetic, that is, whenever . Therefore you could store all values computed until and for computing the value there are two possibilities: is a prime power and then or is composite and with . Then use the stored values to compute .

I now come to the main point of this post: computing the sum of all values of the totient function up to a certain :

One approach is to compute each and sum them. I will call this the brute-force approach. For all numerical purposes I will use Pari-GP in this post. On my computer it takes less than a second to compute and about seconds to compute . This is super linear in time, since the algorithm computes the factorization for each and then sums the values. Using the sieve approach could improve the timing a bit, but the algorithm is still super linear.

In some Project Euler problems it is not uncommon to have to compute something like or even larger. Therefore, there must be more efficient ways to compute out there, so let’s study some of the properties of . In another post I dealt with the acceleration of the computation of the sum of the divisor function.

We have which is the number of pairs with such that . It is not difficult to see that the total number of such pairs is . Moreover, the possible values of are . Now, if for we search instead for pairs satisfying then we have with and we get

There fore the number of pairs with gcd equal to is . Now we arrive at an interesting recursive formula:

At a first sight this looks more complicated, but there is a trick to keep in mind whenever you see a summation over of terms of the form : these quantities are constant on large intervals. Indeed,

Therefore we can change the index of summation from to . The range of for which the interval contains more than one integer is of order . Indeed, . Therefore for we should have at least one integer in the interval . The part where is larger than corresponds to smaller than . Therefore, we can split into two sums, each of order . and get that

where in the last sum we must make sure that in order to avoid duplicating terms in the sum.

Therefore we replaced a sum until to two sums with upper bound . The complexity is not , but something like since we have a recursive computation. Nevertheless, with this new formula and using memoization, to keep track of the values of already computed, we can compute very fast:

is computed instantly (vs second with brute force)

takes second (vs seconds with brute force)

takes seconds (vs over minutes with brute force)

takes seconds

takes about minutes

etc. Recall that these computations are done in Pari GP, which is not too fast. If you use C++ you can compute in seconds, in second and in seconds and in under a minute, if you manage to get past overflow errors.

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: