Ganesh Swami

One of the many things I learnt from my nonlinear physics project is the pseudospectral method. I’m going to try explaining this method the way I understand it.

Any function can be expressed as an infinite sum of the set of basis functions:

For the purpose of an example, if we take the polynomial basis of order 2, then the function can be approximated as

The coefficients can be solved for by building a system of equations. The three unknowns require three points on the interval where we know the exact value of the function. Let these three points be x0, x1 and x2 and the value of the function at these points be b0, b1 and b3.

The system of equations is now,

written compactly as

The system can be solved by inverting the matrix A, which is of complexity . The function p(x) is an interpolating polynomial in the interval of interest. As you can see, using a higher order polynomial quickly becomes prohibitively expensive.

Instead of picking the polynomial basis, you can pick any other basis. The Chebyshev basis is commonly used. A friend of mine picked the Fourier-Bessel basis. Depends on the application.

Fourier basis

The Fourier basis is fundamental to linear systems theory. Linear operators become diagonal in the Fourier basis. This fact comes handy in numerical schemes such as exponential integrators, where you’ll have to take the exponential of a matrix.

The second advantage falls out of the uncertainty principle. Increasing the number of modes in the Fourier domain is equivalent to shrinking the spacing between steps in the space domain. This leads to an exponential convergence in the number of points required as opposed to the usual polynomial convergence with finite differences.

The third key point in the use of the Fourier basis is that the FFT lets us do the equivalent of a matrix inversion in . This is very fast.

Summary

In summary, three random facts come together to form a beautiful theory:

• The FFT is much faster than a matrix inversion.
• Linear operators are diagonal in the Fourier basis.
• Exponential convergence of spectral methods.

Hope you enjoyed reading this post as much as I enjoyed writing it.

Suppose you had a process that generated data every second, but only data on the timescale of weeks or months is useful. Some of my friends recently used RRDtool for this. This tool is extremely useful for any kind of time-series and is available under a liberal license (GPL.)

Quoting from the tutorial,

RRDtool refers to Round Robin Database tool. Round robin is a technique that works with a fixed amount of data, and a pointer to the current element. Think of a circle with some dots plotted on the edge — these dots are the places where data can be stored. Draw an arrow from the center of the circle to one of the dots — this is the pointer. When the current data is read or written, the pointer moves to the next element. As we are on a circle there is neither a beginning nor an end, you can go on and on and on. After a while, all the available places will be used and the process automatically reuses old locations. This way, the dataset will not grow in size and therefore requires no maintenance. RRDtool works with with Round Robin Databases (RRDs). It stores and retrieves data from them.

I recently came across an article on re-evaluating classic algorithms in the modern age. The article makes a strong point: the relative cost of computer-vs-programmer time has dramatically dropped in recent years. It isn’t capitally efficient to hand tune certain classes of algorithms anymore. It’s a good read and highly recommended.

On a similar note, I’m having doubts about a PDE solver that I’m writing for my course project. I’m implementing the Pseudo-Spectral Method for the space dynamics. The technique is quite simply a projection onto another basis that requires fewer coefficients in the series expansion. At each time step, I take the forward and backward Fourier transform of the space variables. This might have made sense a few years back when the algorithmic complexity of the FFT beats a space convolution hands down, but the added code complexity is probably not worth it. The only reason I’m doing it now is because I’ve never done anything like this before and it’s pretty cool.