Archive

Archive for the ‘matlab’ Category

Project Euler Problem 285

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.

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…

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)



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