### Archive

Posts Tagged ‘MATLAB’

## 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

Categories: Geometry, matlab, Uncategorized

## 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

$\displaystyle v\ x\ y\ z\ R\ G\ B$

where ${v}$ is exactly the character ${v}$, ${x,y,z}$ give the coordinates of the points and ${R,G,B}$ give the colour associated to the point in the RGB format. The face information can be entered in a similar fashion:

$\displaystyle f\ t_1\ t_2\ t_3$

where ${t_1, t_2, t_3}$ 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.

Categories: matlab, Uncategorized

## SIAM 100 digit challenge no 4 – Matlab solution

In 2002 the SIAM News magazine published a list of ${10}$ problems in numerical analysis. The answer to each of the problems should be given with an accuracy of at least ${10}$ digits. And so there came the Hundred-Dollar, Hundred-Digit Challenge. I present below my attempt to solving the ${4}$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

$\displaystyle e^{\sin(50x)}+\sin(60e^y)+\sin(70\sin x)+\sin(\sin(80y))-\sin(10(x+y))+\frac{1}{4}(x^2+y^2).$

Categories: matlab, Optimization

## 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 ${10\times 10\times 10}$ (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 ${h}$.

In order to find this volume, it suffices to subtract from ${h\times \text{base area}}$ the volumes of the portions of the spheres which lie under the plane of height ${h}$. Thus, we only need the formula for the volume of a spherical cap, which is ${V_{cap} = \pi (h^2R - h^3/3 )}$, where ${R}$ is the radius of the sphere, and ${h}$ is the height of the cap. Such a function, which computes the volume under a plane of height ${h}$ 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 ${V/s,2V/s,...,(s-1)V/s}$. 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 ${100}$ slices takes ${0.1}$ 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 ${T}$ and a discretization of the interval ${[0,T]}$ by dividing it into ${N}$ 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 ${t_0}$ for which this area is maximized. Now, even if the discretization step ${T/N}$ is large (greater than the demanded precision of ${10^{-3}}$), we can conclude that ${t_0}$ is an approximation of the final time with an error smaller than ${T/N}$. 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 ${T}$. It is possible that if ${T}$ or ${N}$ 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