I went to a magic show the other day. Nick Trefethen was giving a demo of Chebfun, a Matlab extension package he is building in collaboration with his Oxford students and colleagues. In the course of the talk, several mathematical rabbits were pulled out of numerical hats.
The key idea in Chebfun is to represent any function of a real variable by a polynomial approximation.
>> f = chebfun('sin(x) + sin(x.^2)', [0 10]);
>> plot(f)
That wiggly line looks like a graph of y = sin(x) + sin(x2), but that’s an illusion. What is being plotted here is a certain polynomial of degree 118 that happens to approximate sin(x) + sin(x2) with high precision.
As I understand it, the chebfun construction algorithm works something like this. First you select N+1 points in the interval where the function is defined, and construct the unique polynomial of degree N that passes through all the points. If the error of this approximation is below a threshold, you’ve found your chebfun. Otherwise, choose a larger sample of points and try again.
The sample points are not evenly spaced across the interval. They are Chebyshev points, whose distribution varies as a cosine function, denser at the extremes and sparser in the middle. In this case, the process converged with 119 Chebyshev points:
In one respect the example above is an easy one: The function is quite smooth. Here’s something more challenging:
>> hat = 1-abs(x-5)/5; >> h = max(f, hat); >> plot(h)
This is where we pull the rabbit out of the hat—or at least several pairs of rabbit ears. To deal with the discontinuities sharp corners in this curve, the Chebfun system assembles 25 polynomial segments, each defined on a different interval. Some are linear, some of higher degree. But the entire structure is still treated as a single function, which can be operated on by other functions. For example, sum(h) calculates the integral over [0, 10], returning the result 8.598303617326401. And here’s the square root of those rabbit ears:
These are neat tricks, but why would one want to work with polynomial approximations to a function, rather than with the function itself? I’m too new to all this to answer that question with confidence, so I’ll quote the Chebfun Guide:
The aim of Chebfun is to “feel symbolic but run at the speed of numerics”. More precisely our vision is to achieve for functions what floating-point arithmetic achieves for numbers: rapid computation in which each successive operation is carried out exactly apart from a rounding error that is very small in relative terms.
For those who want to know more, I offer a few pointers:
The first published paper on Chebfun:
Battles, Zachary, and Lloyd N. Trefethen. 2004. An extension of MATLAB to continuous functions and operators. SIAM Journal on Scientific Computing 25:1743–1770. (PDF)
Trefethen’s argument favoring floating-point arithmetic over symbolic computation or exact rational arithmetic:
Trefethen, Lloyd N. 2007. Computing numerically with functions instead of numbers. Mathematics in Computer Science 1:9–19. (PDF)
A provocative account of why polynomial approximation is not as wonky as you may think:
Trefethen, Lloyd N. 2011. Six myths of polynomial interpolation and quadrature. Mathematics Today. (PDF)
Finally, Trefethen has a forthcoming book on Chebfun and related matters (which I have only just begun to read):
Trefethen, Lloyd N. To appear. Approximation Theory and Approximation Practice. (PDF)
Chebfun runs inside Matlab, the numerical computing environment from Mathworks. Chebfun itself has recently become open-source software (under a BSD license), but Matlab is proprietary. As far as I can tell, Chebfun does not not (yet?) run under Octave, the open-source alternative to Matlab.