CHEBFUN GUIDE 1: INTRODUCTION

Lloyd N. Trefethen, April 2008

Contents

1.1 What is a chebfun?

A chebfun is a function of one variable defined on an interval [a,b]. The syntax for chebfuns is almost exactly the same as the usual Matlab syntax for vectors, with the familiar Matlab commands for vectors overloaded in natural ways. Thus, for example, whereas sum(f) returns the sum of the entries when f is a vector, when f is a chebfun it returns a definite integral.

The aim of the chebfun system 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 [Trefethen 2007].

The implementation of chebfuns is based on the mathematical fact that smooth functions can be represented very efficiently by polynomial interpolation in Chebyshev points, or equivalently, thanks to the Fast Fourier Transform, by expansions in Chebyshev polynomials. For a simple function, 20 or 30 points often suffice, but the process is stable and effective even for functions complicated enough to require 1000 or 1,000,000 points. The chebfun system makes use of adaptive procedures that aim to find the right number of points automatically so as to represent each function to roughly machine precision (about 15 digits of relative accuracy).

The mathematical foundations of the chebfun system are for the most part well established by results scattered throughout the 20th century. A key early figure, for example, was Bernstein in the 1910s. Nevertheless it is hard to find the relevant material collected in one place. A new reference on this subject will be [Trefethen 2009]. A printed book adapted from this Guide is also being prepared [Trefethen et al. 2009].

The chebfun system was originally created by Zachary Battles and Nick Trefethen at Oxford during 2002-2005 [Battles & Trefethen 2004]. Battles left the project in 2005, and meanwhile three new members have been added to the team: Ricardo Pachon (from 2006), Rodrigo Platte (from 2007), and Toby Driscoll (from 2008).

1.2 Constructing simple chebfuns

The "chebfun" command constructs a chebfun from a specification such as a string or an anonymous function. If you don't specify an interval, then the default interval [-1,1] is used. For example, the following command makes a chebfun corresponding to cos(20x) on [-1,1] and plots it.

  f = chebfun('cos(20*x)');
  plot(f)

From this little experiment, you cannot see that f is represented by a polynomial. One way to see this is to find the length of f:

  length(f)
ans =
    51

Another is to remove the semicolon that suppresses output:

  f
f = 
   chebfun column (1 smooth piece)
          interval          length   values at Chebyshev points
(       -1,        1)        51         0.41       0.44       0.55 ...       0.41
 

These results tell us that f is represented by a polynomial interpolant through 51 Chebyshev points, i.e., a polynomial of degree 50. These numbers have been determined by an adaptive process. We can see the data points by plotting f with the '.-' option:

  plot(f,'.-')

The formula for N+1 Chebyshev points in [-1,1] is

         x(j) = -cos(j pi/N) ,     j = 0:N,

and in the figure we can see that the points are clustered accordingly near 1 and -1. Note that in the middle of the grid, there are about 5 points per wavelength, which is evidently what it takes to represent this cosine to 15 digits of accuracy. For intervals other than [-1,1], appropriate Chebyshev points are obtained by a linear scaling.

The curve between the data points is the polynomial interpolant, which is evaluated by the barycentric formula introduced by Salzer [Berrut & Trefethen 2004, Salzer 1972]. This method of evaluating polynomial interpolants is stable and efficient even if the degree is in the millions [Higham 2004].

What is the integral of f from -1 to 1? Here it is:

  sum(f)
ans =
   0.091294525072762

This number was computed by integrating the polynomial (Clenshaw-Curtis quadrature -- see Section 2.1), and it is interesting to compare it to the exact answer from calculus:

  exact = sin(20)/10
exact =
   0.091294525072763

Here is another example, now with the chebfun defined by an anonymous function instead of a string. (The construction process is often faster if you use anonymous functions.) In this case the interval is specified as [0,100].

  g = chebfun(@(t) besselj(0,t),[0,100]);
  plot(g), ylim([-.5 1])

The function looks complicated, but it is actually a polynomial of surprisingly small degree:

  length(g)
ans =
    89

Is it accurate? Well, here are three random points in [0,100]:

  x = 100*rand(3,1)
x =
  60.640308996562851
  44.279610499164477
  17.326317266735106

Let's compare the chebfun to the true Bessel function at these points:

  exact = besselj(0,x);
  error = g(x) - exact;
  [g(x) exact error]
ans =
  -0.101108148158334  -0.101108148158334  -0.000000000000000
   0.105741737901417   0.105741737901417   0.000000000000000
  -0.129935754450782  -0.129935754450782   0.000000000000000

If you want to know the first 5 zeros of the Bessel function, here they are:

  r = roots(g); r = r(1:5)
r =
   2.404825557695895
   5.520078110286203
   8.653727912911036
  11.791534439014139
  14.930917708487751

Notice that we have just done something nontrivial and potentially useful. How else would you find zeros of the Bessel function so readily? As always with numerical computation, we cannot expect the answers to be exactly correct, but they will usually be very close. In fact, these computed zeros are accurate to about 13 digits:

  besselj(0,r)
ans =
   1.0e-13 *
  -0.630671148896423
  -0.365160041970356
  -0.064984650697585
  -0.330000027195769
   0.073482599805332

Most often we get a chebfun by operating on other chebfuns. For example, here is a sequence that uses plus, times, divide, and power operations on an initial chebfun "x" to produce a famous function of Runge:

  x = chebfun('x');
  f = 1./(1+25*x.^2);
  length(f)
  clf, plot(f)
ans =
   110

1.3 Operations on chebfuns

There are more than 100 commands that can be applied to a chebfun. For a complete list you can type "methods":

  methods chebfun
Methods for class chebfun:

abs           conv_inf      flipud        norm          sech          
acos          conv_tmp      floor         not           set           
acosh         cos           get           null          sign          
acot          cosh          gmres         ode113        sin           
acoth         cot           horzcat       ode15s        sinh          
acsc          coth          imag          ode45         size          
acsch         csc           interp1       orth          spy           
asec          csch          inv           pinv          sqrt          
asech         ctranspose    isempty       plot          std           
asin          cumprod       isequal       plot3         subsasgn      
asinh         cumsum        isreal        plus          subspace      
atan          define        ldivide       poly          subsref       
atanh         diag          legpoly       polyfit       sum           
bvp4c         diff          length        power         svd           
bvp5c         display       log           prod          tan           
ceil          domain        log10         qr            tanh          
cf            eq            log2          range         times         
chebfun       erf           max           rank          transpose     
chebpade      erfc          mean          rdivide       uminus        
chebpoly      erfcx         min           real          uplus         
chebpolyplot  erfinv        minandmax     remez         vander        
comet         exp           minus         repmat        var           
complex       feval         mldivide      restrict      vertcat       
cond          find          mrdivide      roots         
conj          fix           mtimes        round         
conv          fliplr        ne            sec           

To find out what a command does, you can use "help".

  help chebfun/domain
  DOMAIN   Domain of definition.
  I = DOMAIN(F) returns the domain of definition of the chebfun F. This
  includes breakpoint information if F is a single chebfun, but not if F is
  a quasimatrix.
  
  [A,B] = DOMAIN(F) returns the endpoints of the domain as scalars.
 
  See also domain/domain.
 
  See http://www.comlab.ox.ac.uk/chebfun for chebfun information.

Every command in the list exists in ordinary Matlab, with the exceptions of "domain", "restrict", "chebpoly", "define", and the chebfun constructor "chebfun". We have already seen "length" and "sum" in action. In fact we have already seen "subsref" too, since that is the Matlab command for (among other things) evaluating arguments in parentheses. Here is another example of its use:

  f(0.5)
ans =
   0.137931034482759

Here for comparison is the true result:

  1/(1+25/4)
ans =
   0.137931034482759

In this Runge function example, we have also implicitly seen "times", "plus", "power", and "rdivide", all of which have been overloaded from their usual Matlab uses to apply to chebfuns.

In the next part of this tour we shall explore many of these commands systematically. First, however, we should see that chebfuns are not restricted to smooth functions.

1.4 Piecewise smooth chebfuns

Many functions of interest are not smooth but piecewise smooth. In this case a chebfun may consist of a concatenation of smooth pieces, each with its own polynomial representation. Each of the smooth pieces is called a "fun", and funs are implemented as a subclass of chebfuns. This enhancement of the chebfun system was developed initially by Ricardo Pachon during 2006-2007, then also by Rodrigo Platte starting in 2007. Essentially funs consist of the "classic chebfuns" for smooth functions on [-1,1] originally implemented by Zachary Battles.

Later we shall describe the options in greater detail, but for the moment let us see some examples. One way to get a piecewise smooth function is directly from the constructor, taking advantage of its capability of automatic edge detection. For example,

  splitting on
  f = chebfun('abs(x-.3)'); plot(f,'.-')

The "length" command reveals that f is defined by four data points, namely two for each linear interval:

  length(f)
ans =
     4

We can can the definition of f in detail by typing f without a semicolon:

  f
f = 
   chebfun column (2 smooth pieces)
          interval          length   values at Chebyshev points
(       -1,      0.3)         2          1.3          0
(      0.3,        1)         2            0        0.7
Total length = 4 
 

This output confirms that f consists of two funs, each defined by two points and two corresponding function values. We can see its structure from another angle with "struct", Matlab's command for seeing the various fields within an object:

  struct(f)
ans = 
     funs: [1x2 fun]
    nfuns: 2
      scl: 1.300000000000000
     ends: [-1 0.300000000000000 1]
     imps: [1.300000000000000 0 0.700000000000000]
    trans: 0

This output again shows that f consists of two funs with breakpoints at -1, 1, and a number very close to 0.3. The "imps" field refers to "impulses", which relate to values at breakpoints, including possible information related to delta functions, discussed in Section 2.4.

Another way to make a piecewise smooth chebfun is to construct it explicitly from various pieces. For example, the following command specifies three functions x^2, 1, and 4-x, together with a vector of endpoints indicating that the first function applies on [-1,1], the second on [1,2], and the third on [2,4]:

  f = chebfun('x.^2',1,'4-x',[-1 1 2 4]);
  plot(f)

We expect f to consist of three pieces of lengths 3, 1, and 2, and this is indeed the case:

  f
f = 
   chebfun column (3 smooth pieces)
          interval          length   values at Chebyshev points
(       -1,        1)         3            1          0          1 
(        1,        2)         1            1 
(        2,        4)         2            2          0
Total length = 6 
 

Our eyes see pieces, but to the chebfun system, f is just another function. For example, here is its integral.

  sum(f)
ans =
   3.666666666666667

Here is an algebraic transformation of f, which we plot in another color for variety.

  plot(1./(1+f),'r')

Some chebfun commands naturally introduce breakpoints in a chebfun. For example, the "abs" command first finds zeros of a function and introduces breakpoints there. Here is a chebfun consisting of 6 funs:

  f = abs(exp(x).*sin(8*x));
  plot(f)

And here is an example where breakpoints are introduced by the "max" command, leading to a chebfun with 13 pieces:

  f = sin(20*x);
  g = exp(x-1);
  h = max(f,g);
  plot(h)

As always, h may look complicated to a human, but to the chebfun system it is just a function. Here are its mean, standard deviation, minimum, and maximum:

  mean(h)
ans =
   0.578242020778010
  std(h)
ans =
   0.280937455806246
  min(h)
ans =
   0.135335283236613
  max(h)
ans =
   1.000000000000000

A final note about piecewise smooth chebfuns is that the automatic edge detection or "splitting" feature can be turned on and off. (Earlier we used "splitting on" to make sure we were in the right mode.) For example, if we compute sin(x) over [0,1000] with automatic edge detection, we end up with a chebfun with many pieces:

  tic, f = chebfun('sin(x)',[0 1000*pi]); toc
  struct(f)
Elapsed time is 0.462871 seconds.
ans = 
     funs: [1x32 fun]
    nfuns: 32
      scl: 1
     ends: [1x33 double]
     imps: [1x33 double]
    trans: 0

In such a case it is more efficient -- and more interesting mathematically -- to turn off the splitting and construct one global chebfun:

  splitting off
  tic, f2 = chebfun('sin(x)',[0 1000*pi]); toc
  struct(f2)
Elapsed time is 0.020192 seconds.
ans = 
     funs: [1x1 fun]
    nfuns: 1
      scl: 0.999999998433811
     ends: [0 3.141592653589793e+03]
     imps: [0 -3.214166459275634e-13]
    trans: 0

On the other hand with splitting off there is no hope for a function like the absolute value:

  warning on
  tic, f = chebfun('abs(x-.3)'); toc
Warning: Function not resolved, using 65537 pts. Have you tried 'splitting on'? 
Elapsed time is 0.124180 seconds.

So we had better turn it on again. In a chebfun computation, splitting can be turned on and off freely to handle different functions appropriately. The default or "factory" value is splitting off.

  splitting on
  tic, f = chebfun('abs(x-.3)'); toc
Elapsed time is 0.036519 seconds.

1.5 Rows, columns, and quasimatrices

Matlab doesn't only deal with column vectors: there are also row vectors and matrices. The same is true of chebfuns. The chebfuns shown so far have all been in column orientation, which is the default for the chebfun system, but one can also take the transpose, compute inner products, and so on:

  x = chebfun('x')
x = 
   chebfun column (1 smooth piece)
          interval          length   values at Chebyshev points
(       -1,        1)         2           -1          1
 
  x'
ans = 
   chebfun row (1 smooth piece)
          interval          length   values at Chebyshev points
(       -1,        1)         2           -1          1
 
  x'*x
ans =
   0.666666666666667

One can also make matrices whose columns are chebfuns or whose rows are chebfuns, like this:

  A = [1 x x.^2]
A = 
   chebfun column 1 (1 smooth piece)
          interval          length   values at Chebyshev points
(       -1,        1)         1            1 
 
   chebfun column 2 (1 smooth piece)
          interval          length   values at Chebyshev points
(       -1,        1)         2           -1          1
 
   chebfun column 3 (1 smooth piece)
          interval          length   values at Chebyshev points
(       -1,        1)         3            1          0          1 
 
  A'*A
ans =
   2.000000000000000                   0   0.666666666666667
                   0   0.666666666666667                   0
   0.666666666666667                   0   0.400000000000000

These are called "quasimatrices", and will be discussed in a later section.

1.6 References

[Battles & Trefethen 2004] Z. Battles and L. N. Trefethen, "An extension of Matlab to continuous functions and operators", SIAM Journal on Scientific Computing 25 (2004), 1743-1770.

[Berrut & Trefethen 2005] J.-P. Berrut and L. N. Trefethen, "Barycentric Lagrange interpolation", SIAM Review 46 (2004), 501-517.

[Higham 2004] N. J. Higham, "The numerical stability of barycentric Lagrange interpolation", IMA Journal of Numerical Analysis 24 (2004), 547-556.

[Salzer 1972] H. E. Salzer, "Lagrangian interpolation at the Chebyshev points cos(nu pi/n), nu = 0(1)n; some unnoted advantages", Computer Journal 15 (1972), 156-159.

[Trefethen 2007] L. N. Trefethen, "Computing numerically with functions instead of numbers", Mathematics in Computer Science 1 (2007), 9-19.

[Trefethen 2009] L. N. Trefethen, Neoclassical Numerics, book in preparation.

[Trefethen et al. 2009] L. N. Trefethen et al., Chebfun, Chebop, and Chebyshev Technology, book in preparation.