Devnull Functions

    Starting in April 2008 I began posting a periodical article on the STWing Devnull listserve. In general the posts are not daily, but they are universally functional. In each edition I explain a particular technique in numerical methods and provide sample implementation code (generally MATLAB or Maple). The posts are archived here to improve their availability.

1. 4-11-2008, Calculating Pi
Today's function calculates pi. There is an infinite quotient law that goes like:

pi + 3 = 6 + 1^2/( 6 + 3^2 / (6 + 5^2 / (6 + 7^2 /... ) ) )

This is harder to implement than a summed series (which would probably be a smarter way to find pi). The book-keeping can be a little tricky, but I took care of the hard part for you:
Maple Worksheet
Command Line Maple Text

I'll ruin the suspense for you guys, after 2 million iterations, I got:
3.1415926535897932384313933520294853060721694077979578717
457499511945592313174818062610473589013222815416639175869
294742947798883313542744127200580870808734684903576358895
914262334879400086540973263037

This is off by:
(10^-19)*0.31250031250017578124999991577147949229194641113
257174968727192366987466440794786440484168926352832352313
956278196307819005301321047607648959926466812347094872071
7046142948576040783
so 19 decimal places is pretty crappy for 2 million iterations, this is why people don't use infinite quotients for calculations.

2. 4-11-2008 Calculating 2
In the theme of calculating important numbers (yesterday it was pi) today we will calculate the square root of 4.

Let us first form this problem as:
4^(1/2) = x
exp( (1/2) * ln(4) ) = x

Now to calculate ln(4) we can use a McLauren series for ln(1-x). However, this only converges for abs(x) < 1. Thus we divide by e until 4_shrunk is less than one (this is i times). Thus:

ln(4) = ln (4_shrunk) + i

The taylor series can then be calculated as
ln(4_shrunk) = ln(1 - (1 -4_shrunk)) = sum_n ( 4_shrunk^n / n)

This is multiplied by 1/2 and then the series solution for e^x can be used.

Anyway, shrinking by a factor of e^2, and using 10 elements in the natural log and exponential series I got:

2.0000297346062226

This is probably sufficient for most of your engineering needs. (Hint, this algorithm can be used generally if you don't have a pow function)

MATLAB Code

3. 4-13-2008 Solving Universal Gravitation Problems
Tonight we will algebraically solve the n-body problem in gravitation. Actually, it turns out that we can't solve anything higher than the 2-body problem algebraically. Even if I figured out a way to do so, I certainly wouldn't be publishing it on devnull. Instead, we will be solving the n-body problem numerically, using MATLAB.

The n-body problem can be expressed as a system of ordinary differential equations. If we cleverly apply reduction of order, we can characterize each body using 6 first order differential equations. For each body:

x, y, z - The position variables
u, v, w - The corresponding velocity terms

Now, the first order differential equations for the positions take the form:
dx/dt = u, dy/dt = v, dz/dt = w

The tricky parts are from the accelerations. For body 2 acting on body 1:

r = (x_2 - x_1)i + (y_2 - y_1)j + (z_2 - z_1)k

F_12 = m_1*a_12 = (G*m_1*m_2 / (r (dot) r))  * (r/norm(r))
a = (Gm_2 / (r (dot) r))  * (r/norm(r))

These accelerations are summed for every other body which acts on each body.

Below is this algorithm implemented in a MATLAB script. To integrate it, use some integrator like ode45 in the following form:
[T, Y] = ode45(@NBodyODE, [Tstart Tend], ICs);

T will be the column vector of the time steps for the integrated solution
Y will be a matrix with columns: m_1 x_1 y_1 z_1 u_1 v_1 w_1 m_2 x_2 y_2 z_2 u_2 v_2 w_2 ...
Tstart and Tend are the initial and final times for the integrator (how long you want to track for)
ICs are a column vector of the initial conditions of the form [m_1i; x_1i; y_1i; z_1i; u_1i; v_1i; w_1i; m_2i; x_2i; y_2i; z_2i; u_2i; v_2i; w_2i;...]

MATLAB Code

4. 4-14-2008 Planar Center of Mass
I have to run soon, so this will be a quick edition. Let's say that you have a matrix of values (like a picture for example), and you want to find the center of mass. Well, if you are in MATLAB, you are in luck, because the first moment can be done very compactly, Observe:

-------------------------------------------------
[height, width] = size(Matrix);
mass = sum(sum(Matrix));
CMx = round( sum( Matrix * [1:width]' ) / mass );
CMy = round( sum( [1:height] * Matrix ) / mass );
-------------------------------------------------

height and width are the matrix sizes
mass is the total value sum for the matrix
CMx and CMy are the indices of the center of mass for the grid

Pretty Cool, no?

5. 6-9-2008 Calculating Pi, again
A friend of mine once suggested that random numbers could be used to calculate
pi. I'm sitting in the fishbowl alone, and don't feel like walking all the way
back to  my apartment on 43rd street yet. Instead I have decided to try out his
theory.

So, first you center a circle with diameter 1 on the origin. Then you
circumscribe a unit square around it. Then you can randomly generate points
inside of the square {[-0.5,0.5], [-0.5,0.5]}. The ratio of points that lie in
the circle (sqrt(x^2 + y^2) < 0.5) to total points generated (f) gives you pi. where
pi = 4*f.

So I tested 1 billion points and got 3.141661

In principle this works, but it's a pretty bad way to calculate PI, or use CPU
cycles. However, it is a great way to break wave equations if you are so inclined.

On an entirely unrelated note, does anyone need to use force in the next few weeks?
MATLAB Code

6. 6-12-2008 Image Compression
Hi all, today's function is a ghetto image compression scheme. The basic principle of most of the lossy image compression formats (ex: jpeg) is the Fourier Series. Basically, all nicely behaved functions can be approximated as a linear combination of sine/cosine functions (really any orthogonal family, but we won't go there for now). So here's the basic idea: use calculus to find these coefficients, and store them instead of the entire image. Then, rebuild the image from the coefficients when you want to look at your pretty picture. JPEG uses this scheme in 8x8 pixel blocks, and does everything in full color. Because I'm lazy, I'm just going to use one block (the whole image) and black/white. I found that this works pretty well for periodic images, and only fair for more general cases.

EX: Grid, this works very well, the 315x315 image (left) is approximated as a 9x9 coefficient set
Original Grid             Compressed Grid

EX: Face, this doesn't work as well, the 368x491 image (left) is approximated as a 20x20 coefficient set (but it takes up 1/450 of the space)
Original Face        Compressed Face

MATLAB Code

7. 7-14-2008 Is it inside?
This edition of the DevNull function will address an important problem that arises in many contexts: finding whether points lie inside of a domain. If you are working in a CAD application for example you can try to generate invalid models where two contours intersect each other. The CAD program must be able to identify this problem so that it doesn't crash from trying to build bad geometries.
    In general this problem is very tricky, and may be On^2 for two curves(you have to check every point on one curve against the points on the other curve. However, an On solution exists for 2D star shaped domains. In this case, one of the contours must have the property that no ray drawn from its center intersects the boundary more than once.
    In the star shaped geometry a transform can be defined that maps the star shaped domain to a circle while keeping all interior points inside and all exterior points outside (this is a theorem). If this transform is applied to the second curve, then intersections can be identified with the circle equation (x^2 + y^2 = 1). Check it out (in MATLAB)
MATLAB Code

8. 7-16-2008 How big is it?
I was continuing my work on 2D contours and I came across a new problem that is closely related to the previous DevNull function. If you are given a series of points in a plane that define a contour, you may need to know the area inside the contour. The On^2 solution would be to mesh the area into rectangles or triangles and sum their area. However, there exists an On method that uses a little trick from vector calculus.
    Green's Theorem in the Plane states for two functions(P,Q) defined over a planar area bounded by a contour:
int_boundary(Pdx + Qdy) = int_area( (dQ/dx - dP/dy)dA)

Now if we are exceedingly clever, we can pick functions like Q=2x and P=y that are defined everywhere, and simplify to:
int_area(dA) = A = int_boundary(ydx + 2xdy)

The path integral is easy to compute, and obviously linear in the number of points on the boundary.

And just to prove that it works for all of you (from the MATLAB command line):
n = 1024; dt = 2*pi/n; t = (dt:dt:2*pi)';
X = [cos(t) sin(t)];
disp(sprintf('%1.10d', getArea(X)));

gives

3.1415729404e+000

which ought to be the correct answer (slightly less than Pi)

MATLAB Code