%% CHEBFUN GUIDE 1: GETTING STARTED WITH CHEBFUN % Lloyd N. Trefethen, October 2009, revised February 2011 %% 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, it returns a definite integral when f is a chebfun. %% % Chebfun with a capital C is the name of the software system. %% % 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 [Trefethen 2007]. %% % The implementation of Chebfun 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. Chebfun 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 Chebfun 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 the Chebfun-based book [Trefethen 2013]. %% % Chebfun was originally created by Zachary Battles % and Nick Trefethen at Oxford during 2002-2005 % [Battles & Trefethen 2004]. Battles left the project % in 2005, and soon four new members were added to the team: Ricardo % Pachon (from 2006), Rodrigo Platte (from 2007), % and Toby Driscoll and Nick Hale (from 2008). % Beginning in 2009, Asgeir Birkisson and Mark Richardson % also became involved. Additional contributors from Oxford % and elsewhere include Phil Assheton, Folkmar % Bornemann, Pedro Gonnet, Tom Maerz, Sheehan Olver, Simon Scheuring, % Alex Townsend, and Joris Van Deun. %% % This Guide is based on Chebfun Version 4, released % in 2010. Chebfun is available at % http://www.maths.ox.ac.uk/chebfun/. %% 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) %% % Another is to remove the semicolon that suppresses output: f %% % These results tell us that f is represented by a polynomial % interpolant through 49 Chebyshev points, i.e., a polynomial of % degree 48. 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) %% % 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 %% % Here is another example, now with the chebfun defined by % an anonymous function instead of a string. % 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) %% % Is it accurate? Well, here are three random points % in [0,100]: x = 100*rand(3,1) %% % 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] %% % If you want to know the first 5 zeros of the Bessel function, % here they are: r = roots(g); r = r(1:5) %% % 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 close to machine precision: besselj(0,r) %% % 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) %% 1.3 Operations on chebfuns % There are more than 200 commands that can be applied to % a chebfun. For a list of many of them you can type "methods": methods chebfun %% % (Futher commands can be unearthed with methods linop, % methods chebop, and methods domain.) % To find out what a command does, you can use "help". help chebfun/chebpoly %% % Most of the commands in the list exist in ordinary Matlab; % some exceptions are "domain", "restrict", % "chebpoly", "define", and "remez". % 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) %% % Here for comparison is the true result: 1/(1+25/4) %% % 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 Chebfun % was developed initially by Ricardo Pachon during 2006-2007, then % also by Rodrigo Platte starting in 2007 [Pachon, Platte and % Trefethen 2009]. % Essentially funs are the "classic chebfuns" for smooth % functions on [-1,1] originally implemented by Zachary Battles in % Chebfun Version 1. %% % 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, % in the default "splitting off" mode a function with a jump in % its derivative produces a warning message, f = chebfun('abs(x-.3)'); %% % The same function can be successfully captured with splitting on: f = chebfun('abs(x-.3)','splitting','on'); %% % The "length" command reveals that f is defined by four data points, % namely two for each linear interval: length(f) %% % We can see the structure of f in more detail by typing % f without a semicolon: f %% % This output confirms that f consists of two funs, each defined % by two points and two corresponding function values. We can see % the structure from another angle with "struct", Matlab's command % for seeing the various fields within an object: struct(f) %% % 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. The "trans" field is 0 for a column chebfun % and 1 for a row chebfun (Section 1.6 and Chapter 6). The "jacobian" and "ID" fields % are used for automatic differentiation (Chapter 10). %% % 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 %% % Our eyes see pieces, but to Chebfun, f is just % another function. For example, here is its integral. sum(f) %% % 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 Chebfun % it is just a function. Here are its mean, standard deviation, % minimum, and maximum: mean(h) %% std(h) %% min(h) %% max(h) %% % A final note about piecewise smooth chebfuns is that % the automatic edge detection or "splitting" feature, when it % is turned on, may subdivide functions even though they do not % have clean point singularities, and this may be desirable or % undesirable depending on the application. For example, % considering sin(x) over [0,1000] with splitting on, % we end up with a chebfun with many pieces: tic, f = chebfun('sin(x)',[0 1000*pi],'splitting','on'); toc struct(f) %% % In this case it is more efficient -- and more interesting % mathematically -- to omit the splitting and construct one % global chebfun: tic, f2 = chebfun('sin(x)',[0 1000*pi]); toc struct(f2) %% % 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; % see Chapter 8. %% 1.5 Infinite intervals and infinite function values % A major change from Chebfun Version 2 to Version 3 was the generalization of % chebfuns to allow certain functions on infinite intervals or which % diverge to infinity: the credit for these innovations belongs to % Nick Hale, Rodrigo Platte, and Mark Richardson. % For example, here is a function on the whole real axis, f = chebfun('exp(-x.^2/16).*(1+.2*cos(10*x))',[-inf,inf]); plot(f) %% % and here is its integral: sum(f) %% % Here's the integral of a function on [1,inf]: sum(chebfun('1./x.^4',[1 inf])) %% % Notice that several digits of accuracy have been lost here. Be careful! -- % operations involving infinities in Chebfun are not always as accurate % and robust as their finite counterparts. %% % Here is an example of a function that diverges to infinity, % which we can capture by including the flag 'blowup 2' (try help blowup % for details): h = chebfun('(1/pi)./sqrt(1-x.^2)','blowup',2); plot(h) %% % In this case the integral comes out just right: sum(h) %% % For more on the treatment of infinities in Chebfun, see Chapter 9. %% 1.6 Rows, columns, and quasimatrices % Matlab doesn't only deal with column vectors: there are % also row vectors and matrices. The same is true of Chebfun. % The chebfuns shown so far have all been in column orientation, which % is the default, but one can also take % the transpose, compute inner products, and so on: %% x = chebfun('x') %% x' %% x'*x %% % One can also make matrices whose columns are chebfuns % or whose rows are chebfuns, like this: A = [1 x x.^2] %% A'*A %% % These are called "quasimatrices", and they are discussed in % Chapter 6. %% 1.7 How this Guide is produced % This guide is produced in Matlab using the "publish" command. % The formatting is rather simple, not relying on TeX features or % other fine points of typesetting. To publish a chapter for % yourself, make sure the chebfun guide directory % is in your path and then type, for example, "open(publish('guide1'))". % Before publishing, we recommend executing "guidedefaults". %% 1.8 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. % % [Pachon, Platte & Trefethen 2009] R. Pachon, R. B. Platte % and L. N. Trefethen, "Piecewise smooth chebfuns", % IMA J. Numer. Anal., 2009. % % [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 201£] L. N. Trefethen, Approximation Theory and % Approximation Practice, SIAM, publication expected in 2013.