Finite elements, Lagrange's equation, and Numerical Python
The following picture shows some approximate eigenfunctions of Lagrange's equation on a disk clamped at the edge — that is, a simulation of some of the modes of vibration of an idealised drum. Green points lie in the plane of the paper, blue points above it and red points below.
The program that generated them is here; it's a short piece of Python 2.3 code using the numarray numeric library, and outputs a large (potentially very large) Postscript file showing the vibrational modes. The calculation proceeds in five stages:
To be more exact, the image consists of N triangles of points (N=5 for the title image, N=6 for the sample results), with the points within the triangles linked in a triangular mesh, and a "zipper" running down to link the sets of triangles together. See the two illustrations below, one of which has one of the five triangles from the title image highlighted, and the other of which, besides demonstrating that with a touch-pad everyone draws like a five-year-old, tries to show how the points within and between triangles are linked.
To do this, we run over the set of points. At each point P, we expand φ in a Taylor series around P, and specialise this Taylor series to each of the points P' that adjoin P. This gives an overdetermined series of linear equations in the Taylor-series coefficients, which we solve by least-squares; the least-squares solution to Ma=b is linear in b, so the best estimates for the d²/dx² and d²/dy² coefficients are linear functions of the values of φ at the P'. We approximate ∇²φ by the sum of these two coefficient estimates.
I'm not sure how I'd make the matrix generation very substantially better.
There are clearly lots of things to do in the way of mesh generation; it should be possible to define an area in the plane and have software fill it with nicely-distributed points (for example, by starting off with randomly-distributed points and then applying some physical law equivalent to having the points repel one another and slide along the edges of the area), then calculating the neighbours of each point using a Voronoi-like algorithm.
The matrices produced from this kind of finite-element method are very sparse, since ∇²φ(P) is calculated from the values of φ at only a few other points near P. There are certainly algorithms for computing eigenvectors of sparse matrices efficiently; whole regions of engineering and simulation are built around them. But they're not implemented in numarray yet, and I know nothing about them; I'd much appreciate some good references. If I could get sparse methods to work efficiently, it might well be possible to work with hundreds of thousands to millions of vertices on my laptop, rather than the two thousand or so that I'm currently constrained to.
There are obvious ways to improve the visualisation, given access to OpenGL or DirectX from within Python; working in two orthogonal directions, you could either visualise the drum-heads in three dimensions by drawing triangles with vertices x,y,φ, or you could produce much nicer two-dimensional drawings by plotting Gouraud-shaded triangles with vertices x,y and vertex colour dependent on phi. But access to OpenGL from within Python appears only to be implemented under the Irix operating system at the moment; and I'm not altogether sure that my laptop speaks OpenGL at all.
Please send comments on this page (things I've forgotten to mention, things I've mentioned that ought to be common knowledge to the merest imbecile, mathematical errors, issues of Python style, notes on Python packages to do interesting things, references to papers about sparse-matrix methods ...) to me at firstname.lastname@example.org.