Hack the derivative

The article comments: "In a lot of respects, it’s quite amazing how accurate many calculations can be made with floating numbers, like orbital mechanics and heat equations"

I'm not sure it is so surprising. In Nick Trefethen's Numerical Analysis entry in The Princeton Companion to Mathematics, he notes:

"Thus, on a computer, the interval [1 , 2] , for example, is approximated by about 10^16 numbers. It is interesting to compare the fineness of this discretization with that of the discretizations of physics. In a handful of solid or liquid or a balloonful of gas, the number of atoms or molecules in a line from one point to another is on the order of 10^8 (the cube root of Avogadro’s number). Such a system behaves enough like a continuum to justify our definitions of physical quantities such as density, pressure, stress, strain, and temperature. Computer arithmetic, however, is more than a million times finer than this. Another comparison with physics concerns the precision to which fundamental constants are known, such as (roughly) 4 digits for the gravitational constant G, 7 digits for Planck’s constant h and the elementary charge e, and 12 digits for the ratio μ_e/μ_B of the magnetic moment of the electron to the Bohr magneton. At present, almost nothing in physics is known to more than 12 or 13 digits of accuracy. Thus IEEE numbers are orders of magnitude more precise than any number in science. (Of course, purely mathematical quantities like π are another matter.)" http://people.maths.ox.ac.uk/trefethen/NAessay.pdf

Had the following thought just to entertain myself, I don't expect it to have any utility in practice, but here goes:

Take a small interval around your smooth function that includes the point where you want to compute the derivative. Joint the right and left ends by just reflecting itself, so that you have a differentiable periodic function. Now sample that at a suitable rate, take its FFT, derive the FFT of the derivative by multiplying it with 2\pi I. Take the inverse FFT and evaluate it at the point of interest. May be useful if you want to compute the derivative at many points, but seems wasteful otherwise.

One of course has to dot the i's and cross the t's (of aliasing effects, Nyquist limits etc). Perhaps some other fast integral transform would work even better. I will be surprised if there isn't an old industry around it.

This is called a spectral method and is widely used in numerical codes.

I highly recommend Spectral Methods in MATLAB by Trefethen (who someone mentions above) for a very good tutorial. You can freely ignore the MATLAB part and use whatever programming language you want, as long as you have an FFT routine.

When you use this method you are implicitly making the function periodic. I can give you any function on some interval (sufficiently well behaved) and you can compute the Fourier series of it. Even though it's only defined on the interval, if you plotted the Fourier series you would still find it to be periodic. The same idea carries over to the spectral derivative. Even if the function isn't periodic, the method still works, it's just that the numerics near the discontinuity are going to be crap and can affect things further away. But if you are sufficiently far away, things are kosher. Note there are ways of dealing with this using Chebyshev interpolation (see http://math.mit.edu/~stevenj/fft-deriv.pdf section 6) which Trefethen's book also discusses.

And no FFTs are the only transformations I am aware of. I've never heard of anyone using other transformations in a general numerical context, outside of specialized problems.