## 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:

1. Generate the mesh: the mesh consists of a single point at the centre and then a number of points arranged in concentric circles. Except for some edge cases, each point is linked to six neighbours, being the two adjacent to it on its circle, two on the next circle in, and two on the next circle out.

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.

2. Generate the matrix: we want to find a matrix which, when multiplied by a vector giving the height φ of the drum-head at each of the points on the mesh, gives a vector whose coefficients are ∇²φ at those points.

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.

3. Find eigenvectors of the matrix: an eigenvector of the above matrix corresponds to a function φ with ∇²φ=λφ. Computing eigenvectors of matrices with inexact coefficients is a difficult task, but thankfully it's a routine difficult task, and numarray.linear_algebra.eigenvectors does it for us, no questions asked. It takes a few seconds for a few hundred vertices, and ten or so minutes for 2000 vertices, on my new laptop PC.
4. Visualise the result: we know the coefficients of the vertices of the mesh, so, for each real-valued eigenvector, I write out Postscript code to draw a circle at each vertex with colour depending on the value of φ at that vertex. The advantage of using Postscript is that, via Ghostscript, it can be arbitrarily scaled, and displayed with edge anti-aliasing.

### What next?

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.