### Archive

Posts Tagged ‘numerical analysis’

## FreeFem++ Tutorial – Part 1

February 24, 2016 2 comments

First of all, FreeFem is a numerical computing software which allows a fast and automatized treatment of a variety of problems related to partial differential equations. Its name, FreeFem, speaks for itself: it is free and it uses the finite element method. Here are a few reasons for which you may choose to use FreeFem for a certain task:

1. It allows the user to easily define 2D (and 3D) geometries and it does all the work regarding the construction of meshes on these domains.
2. The problems you want to solve can be easily written in the program once we know their weak forms.
3. Once we have variables defined on meshes or solutions to some PDE, we can easily compute all sorts of quantities like integral energies, etc.

Before showing a first example, you need to install FreeFem. If you are not familiar with command line work or you just want to get to work, like me, you can install the visual version of FreeFem which is available here. Of course, you can find example programs in the FreeFem manual or by making a brief search on the internet.

I’ll present some basic stuff, which will allow us in the end to solve the Laplace equation in a circular domain. Once we have the structure of the program, it is possible to change the shape of the domain in no time.

## Even small errors can be fatal

Machines and computers represent numbers as sequences of zeros and ones called bits. The reason for doing this is the simplicity of constructing circuits dealing with two states. This fact coupled with limited memory capacities means that from the start we cannot represent all numbers in machine code. It is true that we can get as close as we want to any number with great memory cost when precision is important, but in fact we always use a fixed precision. In applications this precision is fixed to a number of bits (16, 32, 64) which correspond to significant digits in computations. Doing math operations to numbers represented on bits may lead to loss of information. Consider the following addition using ${15}$ significant digits:

$\displaystyle 5.00000000000002+6.00000000000003 = 11.0000000000000,$

notice that using only the first ${15}$ significant digits we have made an error of ${5\times 10^{-14}}$. This may seem small, but if we do not pay attention and we let errors of this kind accumulate we may have a final error which is unacceptable. This is what happened in the case of the Patriot missile case which I’ll discuss briefly below.

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

October 14, 2014 3 comments

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 ${D}$ which contains the shape ${\Omega}$ and add a penalization on the exterior of your domain ${\Omega}$. The problem to solve becomes something like:

$\displaystyle (-\Delta +\mu I) u = 1 \text{ on }D$

where ${\mu}$ is defined by

$\displaystyle \mu(x) = \begin{cases} 0 & x \in \Omega \\ + \infty & x \notin \Omega\end{cases}. \ \ \ \ \ (1)$

Note that doing this we do not need to impose the boundary condition ${u=0}$ on ${\partial \Omega}$. This is already imposed by ${\mu}$, and the fact that ${u}$ is forced to be zero outside ${\Omega}$.

## Simple optimization problem regarding distances

Given ${n}$ points ${A_1,...,A_n}$ in the plane we can study the problem

$\displaystyle \min_M MA_1+...+MA_n.$

For ${n=3}$ the solution is known and for a triangle which has an angle less than ${120^\circ}$ is the Toricelli (or Fermat) point, which lies at the intersection of the circumcircles of the equilateral triangles built outside the triangle having the sides of the triangle as support.

For ${n=4}$, if we have a convex quadrilateral, then a simple triangle inequality argument proves that the optimal position for ${M}$ is the intersection of the diagonals.

In general, we cannot precisely say what is the position of the solution. Providing a numerical algorithm which computes the optimal position of ${M}$ is not a hard task, and I will do that in the end of this post.

## Numerical method – minimizing eigenvalues on polygons

December 23, 2013 1 comment

I will present here an algorithm to find numerically the polygon ${P_n}$ with ${n}$ sides which minimizes the ${k}$-th eigenvalue of the Dirichlet Laplacian with a volume constraint.

The first question is: how do we calculate the eigenvalues of a polygon? I adapted a variant of the method of fundamental solutions (see the general 2D case here) for the polygonal case. The method of fundamental solutions consists in finding a function which already satisfies the equation ${-\Delta u=\lambda u}$ on the whole plane, and see that it is zero on the boundary of the desired shape. We choose the fundamental solutions as being the radial functions ${\Phi_n^\omega(x)=\frac{i}{4}H_0(\omega|x-y_n|)}$ where ${y_1,..,y_m}$ are some well chosen source points and ${\omega^2=\lambda}$. We search our solution as a linear combination of the functions ${\Phi_n}$, so we will have to solve a system of the form

$\displaystyle \sum \alpha_i \Phi_i^\omega(x) =0 , \text{ on }\partial \Omega$

in order to find the desired eigenfunction. Since we cannot solve numerically this system for every ${x \in \partial \Omega}$ we choose a discretization ${(x_i)_{i=1..m}}$ on the boundary of ${\Omega}$ and we arrive at a system of equations like:

$\displaystyle \sum \alpha_i \Phi_i^\omega(x_j) = 0$

and this system has a nontrivial solution if and only if the matrix ${A^\omega = (\Phi_i^\omega(x_j))}$ is singular. The values ${\omega}$ for which ${A^\omega}$ is singular are exactly the square roots of the eigenvalues of our domain ${\Omega}$.

## Best approximation of a certain square root

Let ${\lambda}$ be a real number such that the inequality

$\displaystyle 0 < \sqrt{2002}-\frac{a}{b} < \frac{\lambda}{ab}$

holds for an infinity of pairs ${(a,b)}$ of natural numbers. Prove that ${\lambda\geq 5}$.

## Distance from a point to an ellipsoid

October 17, 2013 5 comments

Here I will discuss the problem of finding the point which realizes the minimal distance from a given ${x_0}$ to an ellipsoid. First I will look at the theoretical properties which will lead us to ${x_0}$ and then I will present a numerical algorithm to find the nearest point.

This problem is interesting in itself, but has some other applications. Sometimes in optimization problems we need to impose a quadratic constraint, so after a gradient descent we always project the result on the constraint which is an ellipsoid.

Let’s begin with the problem, which can be formulated in a general way like this: we are given ${y \in \Bbb{R}^n}$, and a matrix ${A}$, positive semidefinite, which gives an ellipsoid ${x^TAx=1}$. Without loss of generality assume that ${A}$ is diagonal (if not, we diagonalize it using an orthogonal matrix and work in a new coordinate system). The problem becomes

$\displaystyle \min_{x^TAx=1} \|x-y\|^2.$

(note that the square is not important here, but simplifies the computations later)