Thursday, December 15, 2011

2D root finding algorithm using a recursive bisection method for convex quadrilaterals wich allows to find more than 1 zero

Lately I've been occupied in finding critical points in a 2D vector field. These are locations on a X, Y domain where both components of a vector field go to zero, i.e., (U, V) = (0,0). This is the same has having two functions, f(x,y) and g(x,y), and I want to find the set positions (xp,yp) where simultaneously f(xp,yp)=0 and g(xp,yp)=0.
This is not a topic I'm studying and I end up finding these points by mapping the domain with streamlines and verifying for patterns which indicate that a critical point is present (either saddles, nodes or focus). By verifying I meant visual inspection, or creating streamlines and looking at those, which is very boring.

Fortunately I've found the work of A. Globus, C. Levit, and T. Lasinski, where for the 2D case they do this process analytically by fitting bilinear interpolations and verifying the several particular cases (no solution inside the quadrilateral, 1 solution, 2 solutions, imaginary roots or infinite solutions).

Throughout this process I thought of a method of refining a root inside a quadrilateral, akin to the bisection method, but different. Instead of dividing to cut my quadrilateral in order to refine my solution, I create a tree where I further divide my element, letting me search for more roots in other locations. The best way to implement this was with a recursive function, as the quadrilateral divisions and root finding tests are the same to the children of the first element and subsequently. The test to see if a quadrilateral should be refined further is very basic: if the U values at the vertexes are all either >0 or <0, and the same happens to the V values at the vertexes, then this quadrilateral does not have any zero. Of course, there are functions where the topology would destroy this basic technique, such as one that has 4 roots inside the quadrilateral, located in a way that U and V is positive at all vertexes.

However for some functions this would allow to find more than one root. The vector field functions, U and V, may be discrete (leading to the bilinear interpolation case) or have an analytic representation. This would be achieved by changing the ufun and vfun.
Unfortunately it is not robust enough to sort out between repeated zeros (for example, if 2 or 4 children of a quadrilateral have the same vertex which coincides with a root, neither to work well under infinite solutions.

Either way it may be a useful piece of code for the future.

Friday, October 21, 2011

Colorbar titles in matplotlib

In matplotlib sometimes you want to put some title in a colorbar, when producing a filled contour plot or similar. One way is to simply:

Another way is to use the set_title or the set_xlabel of your colorbar axes. For the first:

However in this way I don't know how to set the title at a distance from the axes. For the second, this can be done using the labelpad:

Tuesday, October 11, 2011

Avoid possible division by zero by limiting the divisor

The friendly "division by zero" warning has appeared to almost anyone that had a pocket calculator. Although this can be a nice greeting from our little digital friends that is bidding us to give it some more though on our problem and the required computations, it is generally a pain if one has a large code and has to wait several hours, if not days, to finish a computation. It is sometimes natural that in the progress of the numerical steps needed to solve a PDE where the initial and some boundary conditions are dubious, that the solution needs to wander a bit before giving up and returning those anxiously awaited results.

So, for an operation "a/b" where "b" is may be 0 at some point, a limiter can be applied such that a minimum threshold, dictated by parameter "SMALL", forces the value to be higher than zero and also to avoid underflows.

Welcome

I am human, thus I forget things. Well, sometimes these are quite useful things that in the instances when these are really needed, coincidently it is when these are harder to remember. I have used notebooks for a long time but when these are also needed, they are never around, which is a beautiful fact adding to the arguments of Murphy's law. As such, tired of searching for these multiple pieces of code when a difficulty comes by, I've decided to create this blog to put the stuff that I learn, from other people or by trial and error, so that it is available through some simple mouse clicks.

My daily activity concerns programming to solve problems related with numerical things.
The main programming I work with is FORTRAN 90 and the interpreted languages I use are Python (mainly numpy, scipy, matplotlib), Matlab, Maxima. My preferred shell is Bash and my text editor is good old Vim. Regarding text processing, Latex is a must.