Is it really O(h**4)?

Over the course of a couple posts (1 and 2) we've shown that you can extend the idea of finite difference approximations to include differences in the imaginary direction. That is you can take complex steps. I presented a formula for the second derivative on a compact stencil (using only nearest neighbours) which I claimed was fourth order accurate. But is it really?

Here's a simple convergence study using cos(), whose second derivative is -cos(). We'll use python to loop over an array with the number of points per wavelength we want to use in our approximation and then calculate a residual, or error, between the analytical solution and our approximation.

for ni in n:
x = 2*pi*arange(0,ni)/(ni-1)
h = x[1]-x[0]
f = cos(x)
analytical_d2f = -cos(x)
# take a positive imaginary step
x_step = x + complex(0,h)
f_pc = cos(x_step)
# take a negative imaginary step
x_step = x - complex(0,h)
f_mc = cos(x_step)
# take a positive real step
x_step = x + h
f_ph = cos(x_step)
# take a negative real step
x_step = x - h
f_mh = cos(x_step)

approx_d2f_c = ( f_ph + f_mh - f_pc.real - f_mc.real)/(2*h*h)
approx_d2f_h = ( f_ph + f_mh - 2*f )/(h*h)

resid_c[i] = sqrt( sum( (analytical_d2f - approx_d2f_c)**2 ) )
resid_h[i] = sqrt( sum( (analytical_d2f - approx_d2f_h)**2 ) )
i += 1

The approximated and analytical second derivatives are shown below for 10 points per wavelength using both a standard central difference with only a real step as well as the complex step formula.

The convergence behaviour for the two approximations is shown below.

Obviously we're getting much better convergence with our complex step formula.

More fun with the PIL

That's the Python Imaging Library. I demonstrated a couple of edge detection techniques (gradient and laplacian) using the PIL to simply read in the image files and make the pixel values conveniently available for twiddling. There's lots of neat stuff in the library other than just image reading/writing utilities. The blend method for instance.
I created a gray-scale image of the magnitude of the gradient (the edges in the image), and then blended it with the original color image. The python to do that takes just two lines:

cblended = Image.blend(outimg.convert("RGB"), img, 0.4)
cblended.save("kuai_laying_color_blend.jpg")

Image.blend() returns an image that is a linear interpolation of the first two arguments (images). The third argument controls if the emphasis is placed on the first image (number closer to zero) or the second image (number closer to one).

Makes it look kind of like a water color, no?

Complex Step Laplace

We can use the results from the previous post about complex step finite differences to solve the 2D Laplace's equation. Instead of solving the equation on a single rectangular grid as we would with normal finite differences, we now solve on three parallel grids (shown below).

One of the grids lies on the real x-y plane, one of the grids is offset slightly in the positive imaginary direction (blue arrow in the figure above), and one of the grids is offset slightly in the negative imaginary direction. We care about the real part of those two imaginary solutions for our second derivative approximation:


The central difference takes care of our solution on the real grid, but what about the two imaginary grids? For those we need asymmetric differences, and I'll put those up in another post.

You might recognize the look of the example mesh if you've ever played around with the 3D modeling/animation program Blender.

Laplacian Edge Detection

In Poor Man's ADC Revisit, I demonstrated gradient based edge detection using Python's Imaging Library to read in a jpg. You can also use the second derivative (a 2D Laplacian operator) to detect edges. Very similar to the gradient based scheme, just change the finite differences to give an approximation to the second derivative rather than the first.

# calculate second derivative for all three channels, sum of squares for mag
for y in xrange(2,h-1):
for x in xrange(2,w-1):
rdx = inpix[x+1,y][0] + inpix[x-1,y][0] - 2*inpix[x,y][0]
gdx = inpix[x+1,y][1] + inpix[x-1,y][1] - 2*inpix[x,y][1]
bdx = inpix[x+1,y][2] + inpix[x-1,y][2] - 2*inpix[x,y][2]
rdy = inpix[x,y+1][0] + inpix[x,y-1][0] - 2*inpix[x,y][0]
gdy = inpix[x,y+1][1] + inpix[x,y-1][1] - 2*inpix[x,y][1]
bdy = inpix[x,y+1][2] + inpix[x,y-1][2] - 2*inpix[x,y][2]
outpix2[x,y] = math.sqrt( rdx*rdx + gdx*gdx + bdx*bdx + rdy*rdy + gdy*gdy + bdy*bdy )

To perform image segmentation it would probably be smart to use the pixel value, the gradient (mag and direction) and the second derivative, no sense in throwing away useful information. Gray-scale image of the magnitude of the second derivative (edges are really where the second derivative changes sign, but an edge tends to have extrema in the second derivative on either side of it) shown below.


Using Psyco or F2Py has been shown on similar problems to give several orders of magnitude speed up over a straight Python implementation as shown above.

Poor Man's ADC Revisit

I wrote a bit ago about using a webcam as an analog-to-digital converter. The really quick hack I put together used Octave's built-in edge() function, which outputs a black and white image of the edges.

A little bit better result can be had with Python's Imaging Library. Here's the basic script to read in the image, calculate the x and y derivatives for the red, blue and green channels, and then set the gray level in the output image to the magnitude of the derivative.

# takes two arguments, input image file, output image file
import sys
import math
import Image

img = Image.open( sys.argv[1] )

w, h = img.size

outimg = Image.new("L",img.size)
outpix = outimg.load()

inpix = img.load()
x = 2;
y = 2;
for y in xrange(2,h-1):
for x in xrange(2,w-1):
rdx = ( inpix[x+1,y][0] - inpix[x-1,y][0] )/2
gdx = ( inpix[x+1,y][1] - inpix[x-1,y][1] )/2
bdx = ( inpix[x+1,y][2] - inpix[x-1,y][2] )/2
rdy = ( inpix[x,y+1][0] - inpix[x,y-1][0] )/2
gdy = ( inpix[x,y+1][1] - inpix[x,y-1][1] )/2
bdy = ( inpix[x,y+1][2] - inpix[x,y-1][2] )/2
outpix[x,y] = math.sqrt( rdx*rdx + gdx*gdx + bdx*bdx + rdy*rdy + gdy*gdy + bdy*bdy );

outimg.save(sys.argv[2],"JPEG")

The resulting gray image looks a lot better than the black and white version generated with Octave (because we aren't throwing away so much information from the original image).

Complex Step**2

I wrote a while back about using a complex step to approximate the first derivative of a function. The neat thing being that there was no subtractive cancellation in the approximation formula as there would be in a normal finite difference. This meant that we could choose very small step sizes without suffering the effects of magnified round-off error.

What about higher order derivatives? The simple result for the first derivative was a second order approximation that only required one function evaluation. Can we achieve the same benefits for a second derivative? The answer is not quite, but it's still pretty good.

The usual way to approximate a second derivative with finite differences is to take a central difference. This uses Taylor series expansions about the location of interest in two directions:


Adding the idea of using complex steps gives two additional points to include in the approximation formula:


Now we just have to solve a linear system to figure out what the coefficients are in our approximation formula:

Based on the Taylor series expansions above, we want the real part of the answer (rather than the imaginary part as was the case with the first derivative). The coefficients are 1/2 for the functions with the real step and -1/2 for the functions with the imaginary step (and zero for the function at x). So we weren't able to avoid subtractive cancellation, but we got an explicit 4th order accurate approximation for the second derivative on a compact stencil!

Sparkle Plenty

I just had a pretty neat problem on my Electrodynamics (course text) mid-term. The problem is to solve for the electric field for two charged conducting spheres of different radii connected by a conducting filament so they reach the same potential.

The result (pdf of full solution) is that the electric field is greater near the surface of the smaller sphere even though it carries less of the total charge. The physical implications of this is that the electric field will be high near the parts of a charged conductor that have a small radius of curvature (ie corners). That's why you get arcing and sparking on forks and crumpled aluminum foil that you put in the microwave.

New Maxima Wiki

My favorite computer algebra system has a new wiki, from the Maxima mailing list (with slight edits to make the links work):

I've set up a new Maxima wiki at http://apps.sourceforge.net/mediawiki/maxima/

You should have a SF account to edit the wiki, see

There are some pages recovered from the old wiki. Please, check the old wiki dump before recreating pages (do not write them from scratch). Eventually I'll copy it to the apps.sourceforge.net/mediawiki/maxima.

Note that the MediaWiki markup differs from the PHPWiki markup. Consult the MediaWiki manual

--
Setting Orange, Aftermath 28 YOLD 3174
Alexey Beshenov http://beshenov.ru/

If you haven't used Maxima before, give it a try. It can do much of what the commercial packages (Mathematica, Maple, MathCad) can do, and it's free software.

As I wrote a while back, it can be really handy for developing numerical software in Fortran or Python (with F2Py).

LED Light Board

The first part of my aeroponics project (overview) is turning those blue and red LED Christmas lights into a grow lamp.

The basic design is a 2 dimensional array of alternating red and blue lights mounted to a piece of white foam board.

Here's a plot of the layout (generated in Octave of course):

And here is the light board as built:

You'll notice two empty holes on the right hand side of the board, the wire between the lights was a really tight fit to meet the spacing I chose (should have checked that a little more closely before drilling all those holes).

Those mini-LEDs look pretty bright, but I'm not sure if it is producing enough light to actually grow a plant, might need to add a string or two. The packaging claims that you can safely plug 43 (!!) of the 60 light strings into a single outlet.