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.

Aeroponics

I've always wanted to build a little aeroponics setup ever since I read some articles about NASA developing aeroponic grow chambers to use on really long space missions. Since we've moved from the gulf coast up to the mid-west, it would be especially nice to have some healthy growing things to brighten up the house during the winter.

There are little kitchen appliances you can buy that grow tomatoes and herbs using an aeroponic method, but where's the fun in just buying one? Much better to build it.

First we need a utility tub as the root chamber where the nutrient ladden fog will be applied.

Now that blue or red LED Christmas lights are common and inexpensive, they make a good light source. They have really low power consumption compared to other types of lights. Why not a white light source? Plants don't use the green part of the spectrum (that's why they look green), so it's a waste of energy to generate white light when we could just generate the red and blue that they do use. To turn a string of ordinary Christmas lights into a plant light we can mount the lights on a sheet of white foam board.

Next we need an ultrasonic humidifier to generate the fog to apply to the roots. Rather than just water in our humidifier though, we need to add some water soluble fertilizer.

We also need to modify some plastic pots so that they become "net pots", basically just slice lots of holes in the sides so the roots can grow out into the fog chamber.

That's it. All told it looks like a home-made aeroponics setup would run about $100, the largest chunk of that being the ultrasonic humidifier ($30).

Poor Man's ADC

Can you use a commodity webcam as an el-cheapo analog to digital converter (ADC)? The folks at slashdot want to know so they can reduce their energy consumption.
Here's the example image:
Octave has plenty of image processing functions in the image package from Octave Forge. Here's the edges detected with the built-in edge() function:

Now we need to perform a smart correlation with the gauge needles so we can map the angle of the needle to a decimal value in [0,10). One way to do it would be to zoom in on the center of the gauge:

Then we can perform a simple least squares regression on the non-zero pixel locations to get a rough estimate of the "slope" of the needle.

needleI = pngread( "needle.png" );
[i,j] = find(needleI(1:35,1:35));
y = 40-j;
x = i;
X = [ ones(length(x),1), x ];
A = X'*X;
b = (X')*y;
beta = A\b;

For the first needle pictured above we get a slope of -0.19643. This seems a tad low, so we probably need to correct for the slant in the original image.

Octave script with all the calculations demonstrated above.

Also check out this book or this page for lots of image processing examples in MATLAB/Octave.

Stupid Google Calculator Tricks

How many fortnights would it take light to travel the circumference of the earth?

What about the speed of light in furlongs per fortnight?

They've got the mass of the proton and the mass of the electron so you can easily calculate their ratio.

I guess we've no need for that gigantic CRC Handbook any more...but what if you need to look up a constant and you're on a desert island with no internet connection, better still pack that CRC just to be safe.

SCons for Latex

Latex is really great for publishing papers with nicely typeset equations and graphics. It automates all of the little formatting stuff and lets you spend time developing the content of your paper rather than wasting effort fiddling with format. I've been using GNU Make for years to process my Latex sources into dvi's and pdf's, but it's always been pretty painful and clumsy (mostly because I'm not a Make guru).

The makefiles for a Latex document can get pretty ugly and complicated especially if your document depends on various data or graphics files. SCons is a new and improved build-tool that implements its configuration files as Python scripts, the file is called SConstruct now instead of Makefile. The great thing from my point of view is that it has built-in support for Latex. Making a pdf from my Latex source file is simple:

env = Environment()
env.PDF('main.tex')

The PDF builder is smart enough to know it needs to run pdflatex and then based on the output possibly run bibtex and then pdflatex again to track down undefined references (or not). All in two lines of Python.

The really neat thing is how easy it is to extend the built-in build environments.

env = Environment()
plot_bld = Builder(action = 'gnuplot < $SOURCE > $TARGET,
suffix = '.png'
src_suffix = '.gp')
env.Append(BUILDERS = {'Plot' : plot_bld})
env.PDF('main.tex')

Now my SCons can automatically update my png graphics files based on gnuplot commands in any files ending with .gp.

Practical Development Environments has a nice little write-up in Chapter 5 about SCons.

EM Flash Cards

I just discovered a great new plug-in (in time for mid-terms) for OpenOffice.org Impress called Open Cards. It is uses the titles of the slides as the question and the body of the slide as the answer to the flashcards. It keeps track of your learning state, has a couple of different learning modes, and seems to use some pretty sophisticated training algorithms.

I put together a set of flashcards for electrodynamics if you're interested in trying it out (of course you need Open Office installed to use the file).

Complex Step

Finite differences are cool, but you are limited by subtractive cancellation (section 3.4: subtractive cancellation exercises). What if you want to make a numerical estimate of the derivative of a function, like a 2D Gaussian:


A finite difference approximation would be (in Octave):

dfdx = ( f(x+h) - f(x) )/ h;

Complex step is even cooler, because you don't have any subtraction (there's no difference), so you can choose a very small step size without loosing accuracy due to subtractive cancelation:

complex_step = complex( 0, 1e-320 );

The derivative is approximated by just the imaginary part:

x = x + complex_step;
dfdx = imag( f(x) )/imag( complex_step );

So the derivative with respect to x looks like this:

And the estimated derivative with respect to y:


This approach is really useful for design sensitivity analysis, and since modern Fortran supports complex types we can even use this method for serious number crunching!

2008 Election Outcome

This is not a political blog, but the Economist has a pretty good article out on John McCain that I felt like commenting on. The last line sums it up: "...pity, there are few better men in American politics".

While he is a very different man than his fellow Naval Academy grad Jimmy Carter, both are men of high integrity, who took that academy training in "duty, honor, country" seriously, and lived it in their public lives after their military service. Though I wasn't alive for the Carter presidency, so I only know it was ineffective from reading about it, I have been alive for President Carter's active work with Habitat for Humanity and his attempts at peace brokering around the world. I think it's safe to say that he embodies Christian ideals of service and peace-seeking. While John McCain is not a paragon of Christian values, he is an exemplar of honesty and sacrifice in public service.

It's too bad that such good men are such unsuccessful politicians in our country. Hopefully Senator McCain's public service after his unsuccessful presidential bid will be as good for the world as President Carter's continues to be following his unsuccessful presidency. Either way, they are a couple of Academy grads I'm proud of (even if they did go to the wrong academy).

Go Air Force, Sink Navy, Beat Army!

A winning combination: Intel and Heineken

I just got a little Atom-based Intel board to play around with and it currently resides in a rubermaid container as its case. I've also been enjoying the new Heineken draught kegs, so I have a few empties sitting in my kitchen floor. This led to the obvious question:

Can a little Atom motherboard

fit inside a Heineken Draught keg?

These new Atom based boards are little mini-itx form factor, that's ~6.93x6.93". The little micro ATX power supplies are ~5x4x3" and a hardrive is ~3.5x5x1". So it might be just a tad to close on the diameter to be able to slide the motherboard right in.

That would be a pretty sweet case though...