Graham’s Biggest Little Hexagon

Think of the following problem: What is the largest area of a Hexagon with diameter equal to 1?

As is the case with many questions similar to the one above, called polygonal isoperimetric problems, the first guess is the regular hexagon. For example, the largest area of a Hexagon with fixed perimeter is obtained for the regular hexagon. However, for the initial question, the regular hexagon is not the best one. Graham proved in his paper “The largest small hexagon” that there exists a better competitor and he showed precisely which hexagon is optimal. More details on the history of the problem and more references can be found in Graham’s paper, the Wikipedia page or the Mathworld page.

I recently wanted to use this hexagon in some computations and I was surprised I could not find explicitly the coordinates of such a hexagon. The paper “Isodiametric Problems for Polygons” by Mossinghoff was as close as possible to what I was looking for, although the construction is not explicit. Therefore, below I present a strategy to find what is the optimal hexagon and I will give a precise (although approximate) variant for the coordinates of Graham’s hexagon.

Matrices with entries of the form f(i,j) in Matlab

It is often necessary to define a matrix with entries where are line and column indices. Using loops is possible, but is not the most elegant solution. Also, it may not be the fastest. Here is a way to construct such matrices fast and easy: simply use meshgrid to construct the arrays of indices . Take a look at the example below and construct your own:

N = 10;
[I,J] = meshgrid(1:N,1:N);
A = N+1-max(I,J);

Mesh a hollow sphere

This is an interesting experiment I’ve done today and I’d like to share it. There are nice software out there which allow you to mesh regular surfaces, like the sphere, torus or other examples.

I needed something different: I wanted to mesh a sphere with a hole and an inner surface. This kind of surfaces may be useful for those who want to design objects for 3D printing. When using 3D printers, the cost is usually proportional to the volume of material used, since this also consumes impression time. If you manage to make a hollow object, the cost goes down immediately (as the resistance of the object…). However, to make a hole, you need to consider an inner parallel surface. I’ll not handle this here, since it’s more complicated. I’ll just mesh a hollow sphere.

I did this in Matlab using the nice library Distmesh, which is free and really easy to use. The whole difficulty was making a function whose zero level set is the hollow sphere. Recall that it is possible to find the union of two shapes using the minimum of two level set function and the intersection using the maximum. This is all that you need to understand what I did below. You can see the details in the Matlab code below.

function [p,t] = Hole_Spher_Mesh(mh)

  R = 0.5;
  r = 0.4;
  rh = 0.1;

  fh=@(p) sqrt(min(p(:,1).^2+p(:,2).^2,ones(size(p(:,1)))));
  fh=@(p) sqrt((p(:,1).^2+p(:,2).^2));
  [p,t]=distmeshsurface(@(x)fd(x,R,r,rh),@huniform,mh,[-1.1,-1.1,-1.1;1.1,1.1,1.1]);
  points = p';
  npt = prod(size(points))/3

clf
patch('Faces',t,'Vertices',p,'FaceColor',[0,0.7,0.9],'EdgeColor','none','FaceAlpha',0.5);

function res = fd(p,R,r,rh)

res1 = sqrt(sum(p.^2,2))-R;
res2 = (rh-sqrt(p(:,1).^2+p(:,2).^2));
res3 = p(:,3);
res2 = max(-res2,res3);
res4 = sqrt(sum(p.^2,2))-r;
res2 = min(res2,res4);
res = max(res1,-res2);

And here is the resulting surface for the parameter mh = 0.01 (if you put a higher precision computations will take longer).

Using parfor in Matlab

We all know that loops don’t behave well in Matlab. Whenever it is possible to vectorize the code (i.e. use vectors and matrices to do simultaneous operations, instead of one at a time) significant speed-up is possible. However, there are complex tasks which cannot be vectorized and loops cannot be avoided. I recently needed to compute eigenvalues for some 10 million domains. Since the computations are independent, they could be run in parallel. Fortunately Matlab offers a simple way to do this, using parfor.

There are some basic rules one need to respect to use parfor efficiently:

1. Don’t use parfor if vectorization is possible. If the task is not vectorizable and computations are independent, then parfor is the way to go.
2. Variables used in each computation should not overlap between processors. This is for obvious reasons: if two processors try to change the same variable using different values, the computations will be meaningless in the end.
3. You can use an array or cell to store the results given by each processor, with the restriction that processors should work on disjoint parts of the array, so there is no overlap.

The most restrictive requirement is the fact that one cannot use the same variables in the computations for different processors. In order to do this, the simplest way I found was to use a function for the body of the loop. When using a matlab function, all variables are local, so when running the same function in parallel, the variables won’t overlap, since they are local to each function.

So instead of doing something like

parfor i = 1:N
   commands ...
   array(i) = result
end

you can do the following:

parfor i=1:N
   array(i) = func(i);
end

function res = func(i)
   commands...

This should work very well and no conflict between variables will appear. Make sure to initialize the array before running the parfor, a classical Matlab speedup trick: array = zeros(1,N). Of course, you could have multiple outputs and the output array could be a matrix.

There is another trick to remember if the parpool cannot initialize. It seems that the parallel cluster doesn’t like all the things present in the path sometimes. Before running parfor try the commands

c = parcluster('local');
c.parpool

If you recieve an error, then run

restoredefaultpath
c = parcluster('local');
c.parpool

and add to path just the right folders for your code to work.

Project Euler 607

If you like solving Project Euler problems you should try Problem number 607. It’s not very hard, as it can be reduced to a small optimization problem. The idea is to find a path which minimizes time, knowing that certain regions correspond to different speeds. A precise statement of the result can be found on the official page. Here’s an image of the path which realizes the shortest time:

 

Fireworks – Project Euler 317

You can read the text of the problem here. The idea is that we have an explosion at a given height in a uniform gravitational field (no friction/wind). Supposing that all particles go outwards from the explosion point with a constant initial speed, what is the shape of the body formed by these particles? A hint is given in the picture below with a nice Matlab simulation.

Now what happens if we add a little bit of wind? Here’s the perturbation obtained when adding a 10m/s uniform wind speed vs the initial configuration. Something like this would be a lot more challenging 🙂

Project Euler Problem 144 – Laser light escaping an ellipse

Project Euler – problem 144 – visualization of the solution in Matlab

FreeFem to Matlab – fast mesh import

I recently wrote a brief introduction to FreeFem++ in this post. FreeFem is a software designed for the numerical study of partial differential equations. It has the advantage of being able to easily define the geometry of the domain, construct and modify meshes, finite element spaces and solve problems on these meshes.

I use Matlab very often for numerical computations. Most of the numerical stuff I’ve done (take a look here if you want) was based on finite differences methods, fundamental solutions and other classical techniques different from finite elements methods. Once I started using finite elements I quickly realized that Matlab is not that easy to work with if you want some automated quality meshing. PDEtool is good, but defining the geometry is not easy. There is also a simple tool: distmesh which performs a simple mesh construction for simple to state geometries. Nevertheless, once you have to define finite element spaces and solve problems things are not easy…

This brings us to the topic of this post: is it possible to interface Matlab and FreeFem? First, why would someone like to do this? Matlab is easier to code and use than FreeFem (for one who’s not a C expert…), but FreeFem deals better with meshes and solving PDE with finite elements. Since FreeFem can be called using system commands, it is possible to call a static program from Matlab. FreeFem can save meshes and variables to files. Let’s see how can we recover them in Matlab.

There is a tool called “FreeFem to Matlab” developed by Julien Dambrine (link on Mathworks). There’s also a longer explanation in this blog post. I recently tried to use the tool and I quickly found that it is not appropriate for large meshes. It probably scans the mesh file line by line which makes the loading process lengthy for high quality meshes. Fortunately there’s a way to speed up things and I present it below. I will not cover the import of the data (other than meshes) since the function importdata from the FreeFem to Matlab tool is fast enough for this.

Print your Matlab models in 3D

The emergence of 3D printers opens a whole new level of creation possibilities. Any computer generated model could be materialized as soon as it can be transformed in a language that the 3D printer can use. This is also the case with objects and structures which emerge from various mathematical research topics. Since I’m working on shape optimization problems I have lots of structures that would look nice printed in 3D. Below you can see an example of a 3D model and its physical realization by a 3D printer.

I want to show below how can  you can turn a Matlab coloured patch into a file which can be used by a 3D printer. The first step is to export the Matlab information regarding the position of the points, the face structure and the colours into an obj file format. This is not at all complicated. Vertex information is stored on a line of the form

where is exactly the character , give the coordinates of the points and give the colour associated to the point in the RGB format. The face information can be entered in a similar fashion:

where are the indices of the points in the corresponding face. Once such an obj file is created, it can be imported in MeshLab (a free mesh editing software). Once you’re in MeshLab you should be able to export the structure into any file format you want, which can be understood by a 3D printer (like STL). Once you have the stl file, you can go on a 3D printing website like Sculpteo and just order your 3D object.

SIAM 100 digit challenge no 4 – Matlab solution

In 2002 the SIAM News magazine published a list of problems in numerical analysis. The answer to each of the problems should be given with an accuracy of at least digits. And so there came the Hundred-Dollar, Hundred-Digit Challenge. I present below my attempt to solving the th problem in this series, as it deals to finding the global minimum of a highly oscillating function. I confess that my approach is rather simplistic, but solves the problem. First let’s state the problem below:

Problem 4. What is the global minimum of the function

ICPC 2015 World Final Problem D

This is a solution to problem D from the International Collegiate Programming Contest. The list of the problems can be found here. This is one of the problems which was really appealing for me, from the start. The idea is the following: suppose you have a Swiss cheese cube of (the units are centimeters). This cheese cube has multiple holes in it, and for the sake of simplicity, we assume they are spheres. Suppose also that their centers and their radii are known. The problem is how to cut this piece of cheese into horizontal slices having equal weights (that is equal volumes). The holes are supposed to be contained in the cube, and non-intersecting.

The problem is really simple if you don’t have any holes. Just cut at evenly distributed points along the height of the cheese cube. Things change if you add holes, because cutting at evenly distributed points may lead to slices of different weights, since the holes may be not evenly distributed in the cube.

The approach I had in mind was to forget about the slices for a moment. Working with slices means that you may need to cut two caps off a sphere, and that leads to complications. Instead of finding the slices, we could search for the planes dividing the cheese cube into two pieces of given volume ratio. That is, given a proportion of the volume, find the height of the cut, such that that proportion is below the cut. To do this, we only need to be capable to find the volume of the pierced cube which lies under a plane of height .

In order to find this volume, it suffices to subtract from the volumes of the portions of the spheres which lie under the plane of height . Thus, we only need the formula for the volume of a spherical cap, which is , where is the radius of the sphere, and is the height of the cap. Such a function, which computes the volume under a plane of height can be easily vectorized, and implemented in Matlab.

Once we have the function which computes the volume under a certain height, it suffices to do a binary search, in order to find the heights which give volumes . Once we have the heights, we can find the widths of the slices. I implemented the algorithm in Matlab below. I don’t have some complicated test cases, but in the ones I have, the division into slices takes seconds.

function Prob3_2015

n = 2; %number of spheres
s = 100; %number of slices

r = [10000 40000]; %radii of the spheres

%coordinates of the centers
%the box is parametrized as 100000^3
%and the real dimensions are 100^3 mm^3

ct = [10000 20000 20000;
      40000 50000 60000];
%total volume of the spheres
totalv = 4*pi/3*sum(r.^3);
%actual volume of the cheese cube
actualv = 100000^3-totalv;
slices = zeros(1,s);

for i = 1:s-1
target  = i*actualv/s;

m = 0;
M = 100000;

while M-m>1e-6
   mid  = (m+M)/2;
   amid = slice_vol(mid,ct,r);
   if amid>target
      M = mid;
   else
      m = mid;
   end
end
fprintf('Cumulative slice %d : %9.6f mm\n',i,mid/1000);
slices(i) = mid;
end
slices(end)= 100000;
slices = [0 slices];
disp(' ');
for i = 1:s
  fprintf('Width slice %d : %9.6f mm\n',i,(slices(i+1)-slices(i))/1000);
end

% function calculating volume under a plane
function res = slice_vol(h,ct,r)

if isempty(ct)
   pars = 0;
else 
   pars = min(max(h-(ct(:,3)-r(:)),0),2*r(:));
end;

res  = h*10^10-sum(pi*pars.^2.*r(:)-pi*pars.^3/3);

ICPC 2015 World Final Problem B

This the the solution to Problem B from the International Collegiate Programming Contest. The list of problems can be found here. The idea is to find the maximal area of the intersection of two moving polygons. The inputs give the initial positions of two convex polygons and the vectors giving their speed in the plane.

One idea of a solution is as follows: take a maximal time and a discretization of the interval by dividing it into parts. Iterate over the translations at times given by this discretization, compute at each step the area of the intersection of the two polygons (if there is any intersection at all), and in the end find the time for which this area is maximized. Now, even if the discretization step is large (greater than the demanded precision of ), we can conclude that is an approximation of the final time with an error smaller than . This is due to the fact that the function representing the area of the intersection has two monotonicity intervals, as a consequence of the fact that the polygons are convex. On the first interval, the intersection area is increasing, on the second one it is decreasing. Thus, once we have a discrete maximum, we are close to the real maximum.

Now, all we are left to do in order to achieve any desired precision is to refine the search near this discrete maximum, find a new, better approximation, and refine again, until we have enough precision. Of course, one first problem with this algorithm is the initial search using a maximal time . It is possible that if or are not large enough, then we do not detect any intersection of the two polygons. Thus, an initial guess, based on a mathematical argument is needed in order to reduce the search of the optimal time to an interval which is small enough to have enough initial precision to detect the discrete maximum.

The algorithm presented below, uses Matlab’s predefined function polybool which can compute the intersection of two polygons. Of course, this is an overkill, since dealing with intersections of convex polygons is not that complicated (but still, I didn’t have enough time to play with the problem, in order to provide a more optimized version). I do not treat the search for the initial time interval. As I think about it, I guess a n argument based on finding some line intersections should give us a narrow enough time interval (with some care for the case when the two speed directions are collinear). The algorithm presented below solves the sample cases, but could fail in more general situations.

function Prob1_2015(p1,p2)

% choose initial polygons; see the text of the problem
p1 = [6 3 2 2 4 3 6 6 6 7 4 6 2 2 2];
p2 = [4 18 5 22 9 26 5 22 1 -2 1];

%p1 = [4 0 0 0 2 2 2 2 0 1 1];
%p2 = [4 10 0 10 2 12 2 12 0 -1 1];

np1 = p1(1);
np2 = p2(1);
cp1 = p1(2:1+2*np1);
cp2 = p2(2:1+2*np2);

vp1 = p1(end-1:end);
vp2 = p2(end-1:end);

cp1 = cp1(:);
cp1 = reshape(cp1,2,np1);

xp1 = cp1(1,:);
yp1 = cp1(2,:);
cp2 = cp2(:);
cp2 = reshape(cp2,2,np2);

xp2 = cp2(1,:);
yp2 = cp2(2,:);

%set precision
prec = 1;

%here there should be a clever initial choice of 
%the starting time and stopping time
%this choice works well in this case
start = 0;
stop  = 5;

while prec>1e-6
n = 100;

timex = linspace(start,stop,n);
areas = zeros(size(timex));

for i = 1:n
   t = timex(i);
   xt1 = xp1+t*vp1(1);
   yt1 = yp1+t*vp1(2);
   xt2 = xp2+t*vp2(1);
   yt2 = yp2+t*vp2(2);
   
   [x,y] = polybool('intersection',xt1,yt1,xt2,yt2);
   areas(i) = polyarea(x,y);
end
[m,I] = max(areas);
if m<1e-6
   % if no polygonal intersection is detected
   % there is nothing further to do
   display('never');
   break
end

tapp = timex(I);
fprintf('Precision %f | Time %f\n',prec,tapp);
start = tapp-prec;
stop  = tapp+prec;
prec  = prec/10;
end
clf
%this draws the final position of the polygon
%as well as their intersection
fill(xt1,yt1,'blue')
axis equal
hold on
fill(xt2,yt2,'red')
fill(x,y,'green')
axis equal
hold off

Here’s what you get for the first sample test provided in the questions:

 >> Prob1_2015
    res = 4.193548

ICPC 2015 World Final Problem A

This is the solution to Problem A from the International Collegiate Programming Contest. The list of problems can be found here. This first problem consists simply of reading the parameters of the function defined below, and computing its values on the set . Then, you need to find the maximum decrease in the function values.

The inputs are parameters , the function is

and the values to be considered are . Below is a Matlab code which works, at least for the given Sample Inputs. This should be optimized by removing the for loop. It is instant for but it should be modified to work until .

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;

vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);

mat = zeros(n,1);
for i = 1:n
mat(i) = vals(i)-min(vals(i:end));
end

res = max(mat)

The sample outputs are:

 >> Prob1_2015(42,1,23,4,8,10)  
      104.855110477394

 >> Prob1_2015(100,7,615,998,801,3)
      0

 >> Prob1_2015(100,432,406,867,60,1000)
      399.303812592112

New version which works for large . It suffices to eliminate the computation of the min at each of the phases of the iteration. This solves the problem in 0.2 seconds for

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;
vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);
mat       = zeros(n,1);
lastmax   = vals(1);
for i = 1:n
  lastmax   = max(lastmax,vals(i));
  mat(i)    = lastmax-vals(i);
end
res = max(mat)

Identifying edges and boundary points – 2D Mesh – Matlab

A triangulation algorithm often gives as output a list of points, and a list of triangle. Each triangle is represented by the indexes of the points which form it. Sometimes we need extra information on the structure of the triangulation, like the list of edges, or the list of boundary points. I will present below two fast algorithms for doing this.

Finding the list of edges is not hard. The idea is to go through each triangle, and extract all edges. The algorithm proposed below creates the adjacency matrix in the following way: for each triangle we set the elements (and their symmetric correspondents) to be equal to .

In order to find the points on the boundary (in two dimensions), it is enough to look for the edges which are sides to only one triangle. We can do this using the adjacency matrix. Note that if is the adjacency matrix, then stores the number of paths of length (two sides) between two points of the triangulation. Note that any edge which is not on the boundary will contain the starting and ending point of two paths of length . If is a triangle such that points are on the boundary, then (there is one path of length going through . We also have . Conversely, if and then are connected, and there is a unique path of length going from to . Thus, is an edge on the boundary. Therefore, we just need to identify the indexes such that there exists with .

Below are two short Matlab codes doing these two algorithms. I guess they are close to being optimal, since only sparse and vectorized operations are used.

%p is the list of points
%T is the list of triangles, ap is the number of points
%this computes the adjacency matrix
A = min(sparse(T(:,1),T(:,2),1,ap,ap)+sparse(T(:,2),T(:,3),1,ap,ap)+sparse(T(:,3),T(:,1),1,ap,ap),1);
A = min(A+A',1);
% this finds the boundary points, whose indexes are stored in Ibord
B = A^2.*A==1;
Ibord = find(sum(B,2)>0);

Simple triangle mesh – Matlab code

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 , and define the (only) initial triangle as . Then, at each step, take a triangle from the list , add the midpoints of to the list of points and replace the triangle with the four smaller triangles determined by the vertices and midpoints of . 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).

Solving Poisson’s equation on a general shape using finite differences

One of the questions I received in the comments to my old post on solving Laplace equation (in fact this is Poisson’s equation) using finite differences was how to apply this procedure on arbitrary domains. It is not possible to use this method directly on general domains. The main problem is the fact that, unlike in the case of a square of rectangular domain, when we have a general shape, the boudary can have any orientation, not only the orientation of the coordinate axes. One way to avoid approach this problem would be using the Finite Element Method. Briefly, you discretize the boundary, you consider a triangulation of the domain with respect to this discretization, then you consider functions which are polynomial and have support in a few number of triangles. Thus the problem is reduced to a finite dimensional one, which can be written as a matrix problem. The implementation is not straightforward, since you need to conceive algorithms for doing the discretization and triangulation of your domain.

One other approach is to consider a rectangle which contains the shape and add a penalization on the exterior of your domain . The problem to solve becomes something like:

where is defined by

Note that doing this we do not need to impose the boundary condition on . This is already imposed by , and the fact that is forced to be zero outside .

Two more mind games solved using Matlab

After doing this, solving the Sudoku puzzle with Matlab, I wondered what else can we do using integer programming? Many of the mind games you can find on the Internet can be modelled as matrix problems with linear constraints. I will present two of them, as well as their solutions, below. The game ideas are taken from www.brainbashers.com. Note that in order to use the pieces of code I wrote you need to install the YALMIP toolbox.

1. Three in a row

You have a table and on each row and column you need to have squares coloured with color and squares coloured with color such that no three squares taken vertically or horizontally have the same color. In the begining you are given the colors of some of the square such that it leads to a unique solution satisfying the given properties. The goal is to find the colors corresponding to each small square starting from the ones given such that the end configuration satisfies the given properties.

Solve a Sudoku Puzzle using Matlab

There was a time when I loved to solve Sudoku puzzles. It happened sometimes to find some really difficult puzzles, which I was unable to solve. Having a program which solves Sudoku puzzles is cool, but knowing how it works is cooler.

There are multiple types of solvers. One can imagine to take each free cell in a Sudoku puzzle and associate to it the set of integers in which can lie in there (see the picture below). If the set of possible values contains only one element, then you have found the element corresponding to that cell. If not, then the program needs to make one guess, and see if it gets him to a contradiction. If it arrives at a contradiction, return to the choice step and choose another option, and this needs to be done recursively, which can take huge amounts of time.

One other method which can be implemented in Matlab is to use the matrix structure of the Sudoku table and write the constraints of non repetition on lines, columns and blocks as linear constraints. Once all these constraints are well defined, you can call a semidefinite programming solver to find a solution. A semidefinite programming solver looks to find the vector which optimizes a given cost functions under certain constraints. A good sudoku puzzle has only one solution, and therfore every admissible input which satisfies the constraints is a solution. This means that the cost function doesn’t matter here.

I will present the implementation using YALMIP. You can find another equivalent Matlab implementation here.

Finite difference Method for 1D Laplace Equation

I will present here how to solve the Laplace equation using finite differences. First let’s look at the one dimensional case:

In this case it is easy to solve explicitly the differential equation and get that . However, sometimes it is not so easy or not possible at all to find an explicit solution, and there are some numerical methods which can be used to approximate our solution. In the one dimensional case, or in the case of a rectangular domain or a box in superior dimensions, the method of finite differences can be sucessfully used.