Numerical algorithms are traditionally applied to functions discretized in space, but one may ask what their analogues would be if we could compute with continuous functions directly. This question is no longer just theoretical, thanks to the development of the chebfun system in object-oriented MATLAB. We will show the system in action and ask whether it can live up to the vision of “computing with the feel of symbolics but the speed of numerics”.
This is a very interesting story with some potential implications for the way we teach calculus (see previous discussions, Donald Knuth: Calculus via O notation and Calculus without limits), even if the project appears to be concerned with a technical enhancement of MATLAB.
Conventional MATLAB is built on the most advanced algorithms for vectors and matrices; this is a source of its power, since so much of scientific computing in the end comes down to numerical linear algebra. The idea of chebfun computation is to create a MATLAB class that behaves syntactically like the usual MATLAB vectors. The difference is that with chebfuns, the usual vector commands in MATLAB are overloaded with analogues for continuous functions.
For example, the command
>> f = chebfun(‘real(airy(10*x))’);
calls the chebfun constructor to produce a chebfun object that will behave, up to the usual IEEE double machine precision of about 16 digits, like the Airy function .
How is all this done? The mathematical basis of the chebfun system is a pair of closely related ideas:
- Polynomial interpolation in Chebyshev points implemented by Salzer’s barycentric formula
- Chebyshev expansions implemented by the Fast Fourier Transform
These are combined together with state-of-the-art numerical algorithms for integration, zero finding and other operations. In particular, two features that make the system fast and accurate are
- Adaptive determination of the number of points needed to represent each function to about 16 digits
- Zero finding via eigenvalues of Chebyshev companion matrices with divide-and-conquer recursion.
The result is a system with, we have found, rather astonishing capabilities. For example, the Airy function above may seem complicated, but its chebfun representation to 16-digit precision is a polynomial of degree just 59. The more irregular function
>> f = chebfun(‘sin(3*x)+.5*exp(x).*sin(100*x./(1+3.*x.^2))’)
plotted below only needs a polynomial of degree 237. Even a function represented by a polynomial of degree 100,000 is integrated by sum to full accuracy on our workstation in 0:2 secs.
Some of the mathematics underpinning the extraordinary speed, accuracy, and stability of high order polynomial interpolants is old, and some of it is new, but it is safe to say that until the arrival of the chebfun system few people were aware of these possibilities. The investigations proposed here, despite their classical foundations, will explore poorly charted territory even from an algorithmic point of view. From a software point of view they are completely new.
Basically, a function on a segment is replaced by its appropriate Chebyshef polynomial approximation, but, if it has name, it retains that name.
Now, let us look at related issues in methodology of mathematical teaching. Assume that chebfun approach became an industry standard for (semi-)numeric computations, and you have to teach a calculus course for engineering students. Should you pay attention tot he following two metamathematical factors?
- The way they will be used by your students, functions become intensional objects with names, not extensional entities.
- Uniform convergence, not a limit at a point, starts to dominate the field.