%% CHEBFUN GUIDE 8: CHEBFUN PREFERENCES
% Lloyd N. Trefethen, November 2009, revised February 2011
%% 8.1 Introduction
% Like any software package, Chebfun is based on certain design
% decisions. Some of these are fixed, like the principle of representing
% functions by Chebyshev expansions. Others are adjustable, like the
% maximum number of points at which a function will be sampled before Chebfun
% gives up trying to resolve it. A starting point in exploring
% these matters is to type the command help chebfunpref. (For chebops,
% there is help cheboppref.) Or just to see the list of preferences, you
% can simply execute chebfunpref. Here we execute it with the argument
% 'factory' to ensure that all preferences are set to their factory values:
chebfunpref('factory')
%%
% In this chapter of the Chebfun Guide we explore some of these adjustable
% preferences, showing how special effects can be achieved by modifying
% them. Besides showing off some useful techniques, this review will also
% serve to deepen the user's understanding of Chebfun by poking about a
% bit at its edges.
%%
% A general point to be emphasized is the distinction between creating a
% chebfun directly from the constructor and creating one by operating on
% previous chebfuns. In the former case we can include preferences
% directly in the constructor command, and we recommend this as good
% practice:
f = chebfun('x.^x',[0,1],'splitting','on');
%%
% In the latter case, however, one must (say) turn the preference on and
% off again.
x = chebfun('x',[0,1]);
splitting on
f = x.^x;
splitting off
%% 8.2 domain: the default domain
% Like Chebyshev polynomials themselves, chebfuns are defined by default on
% the domain [-1,1] if no other domain is specified. However, this default
% choice of the default domain can be modified. For example, we can work
% with trigonometric functions on [0,2pi] conveniently like this:
chebfunpref('domain',[0 2*pi])
f = chebfun(@(t) sin(19*t));
g = chebfun(@(t) cos(20*t));
plot(f,g), axis equal, axis off
%% 8.3 splitting: splitting on/off
% Perhaps the preference that users wish to control most often is the
% choice of splitting off or on. Splitting off is the factory default,
% though splitting on was the default in Chebfun Version 2.
%%
% In both splitting off and splitting on modes, a chebfun may consist of a
% number of pieces, called funs. For example, even in splitting off mode,
% the following sequence makes a chebfun with four funs:
chebfunpref('factory');
x = chebfun(@(x) x);
f = min(abs(x),exp(x)/6);
format short, f.ends
plot(f)
%%
% One breakpoint is introduced at x=0, where the constructor recognizes
% that abs(x) has a zero, and two more breakpoints are introduced at
% -0.1443 and at 0.2045, where it recognizes that abs(x) and exp(x)/6 will
% intersect.
%%
% The difference between splitting off and splitting on pertains to
% additional breakpoints that may be introduced in the more basic chebfun
% construction process, when the constructor makes a chebfun solely by
% sampling point values. For example, suppose we try to make the same
% chebfun as above from scratch, by sampling an anonymous function, in
% splitting off mode. We get a warning message:
ff = @(x) min(abs(x),exp(x)/6);
f = chebfun(ff);
%%
% With splitting on, Chebfun's built-in edge detector quickly finds the
% singular points and introduces breakpoints there:
f = chebfun(ff,'splitting','on');
f.ends
%%
% This example involves specific points of singularity, which the
% constructor has duly located. In addition to this, in splitting on mode
% the constructor will subdivide intervals recursively at non-singular
% points when convergence is not taking place fast enough. For example,
% with splitting off we cannot successfully construct a chebfun for the
% square root function on [0,1] (unless we use exponentials or mappings as
% described in the next chapter):
f = chebfun(@(x) sqrt(x),[0 1]);
%%
% With splitting on, however, all is well:
f = chebfun(@(x) sqrt(x),[0 1],'splitting','on');
length(f)
format long, f((.1:.1:.5)'.^2)
%%
% Inspection reveals that Chebfun has broken the interval into a
% succession of pieces, each 100 times smaller than the next:
f.ends
%%
% In this example all but one of the subdivisions have occurred near an
% endpoint, for the edge detector has estimated that the difficulty of
% resolution lies there. For other functions, however, splitting will take
% place at midpoints. For example, here is a function that is complicated
% throughout [-1,1], especially for larger values of x.
ff = @(x) sin(x).*tanh(3*exp(x).*sin(15*x));
%%
% With splitting off, it gets resolved by a global polynomial of rather
% high degree.
f = chebfun(ff);
length(f)
plot(f)
%%
% With splitting on, the function is broken up into pieces, and there is
% some reduction in the overall length:
f = chebfun(ff,'splitting','on');
length(f)
format short, f.ends
%%
% When should one use splitting off, and when splitting on? If the goal is
% simply to represent complicated functions, especially when they are more
% complicated in some regions than others, splitting on sometimes has
% advantages. An example is given by the function above posed on [-3,3]
% instead of [-1,1]. With splitting off, the global polynomial has a degree
% in the tens of thousands:
f3 = chebfun(ff,[-3 3]);
length(f3)
plot(f3)
%%
% With splitting on the representation is much more compact:
f3 = chebfun(ff,[-3 3],'splitting','on');
length(f3)
%%
% On the other hand, splitting off mode has advantages of robustness. In
% particular, operations involving derivatives generally work better when
% functions are represented by global polynomials, and chebops
% for the most part requires this. Also, for educational purposes, it is
% very convenient that Chebfun can be used so easily to study the
% properties of pure polynomial representations.
%% 8.4 splitdegree: degree limit in splitting on mode
% When intervals are subdivided in splitting on mode, as just illustrated,
% the parameter splitdegree determines where this will happen. With the
% factory value splitdegree=128, splitting will take place if a polynomial
% of degree 128 proves insufficient to resolve a fun. Let us confirm for
% the chebfun f constructed a moment ago that the degrees of the individual
% funs are all less or equal than 128:
f.funs
%%
% Alternatively, suppose we wish to allow individual funs to have degree up
% to 512. We can do that like this:
f = chebfun(ff,'splitting','on','splitdegree',512);
length(f)
format short, f.ends
f.funs
%% 8.5 maxdegree: maximum degree
% As just mentioned, in splitting off mode, the constructor tries to make a
% global chebfun from the given string or anonymous function. For a
% function like abs(x) or sign(x), this will typically not be possible and
% we must give up somewhere. The parameter maxdegree, set to 2^16 in the
% factory, determines this giving-up point.
%%
% For example, here's what happens normally if we try to make a chebfun for
% sign(x).
f = chebfun('sign(x)');
%%
% Suppose we wish to examine the interpolant to this function through 50
% points instead of 65537. One way is like this:
f = chebfun('sign(x)',50);
plot(f)
%%
% Notice that no warning message is produced since we have asked explicitly
% for exactly 50 points. On the other hand we could also change the
% default maximum to this number (or more precisely the default degree to
% one less than this number), and then there would be a warning message:
f = chebfun('sign(x)','maxdegree',49);
%%
% Perhaps more often one might wish to adjust this preference to enable use
% of especially high degrees. On the machines of 2011, Chebfun is
% perfectly capable of working with polynomials of degrees in the millions.
% The function abs(x)^(3/2) on [-1,1] provides a nice example, for it is
% smooth enough to be resolved by a global polynomial, provided it is of
% rather high degree:
tic
f = chebfun('abs(x).^1.5','maxdegree',1e6);
lengthf = length(f)
format long, sumf = sum(f)
plot(f)
toc
%%
% (Much more efficient ways of resolving this function, by eliminating the
% singularity, are described in Chapter 9.)
%% 8.6 minsamples: minimum number of sample points
% At the other end of the spectrum, the parameter minsamples determines the
% minimum number of points at which a function is sampled during the
% chebfun construction process, and the factory value of this parameter is
% 9. This does not mean that all chebfuns have length at least 9. For
% example, if f is a cubic, then it will be sampled at 9 points, Chebyshev
% expansion coefficients will be computed, and 5 of these will be found to
% be of negligible size and discarded. So the resulting chebfun is a
% cubic, even though the constructor never sampled at fewer than 9 points.
chebfunpref('factory');
f = chebfun('x.^3');
lengthf = length(f)
%%
% More generally a function is sampled at 9, 17, 33,... points until a set
% of Chebyshev coefficients are obtained with a tail judged to be
% negligible.
%%
% Like any process based on sampling, this one can fail. For example, here
% is a success:
splitting off
f = chebfun('-x -x.^2 + exp(-(30*(x-.5)).^2)');
length(f)
plot(f)
%%
% But if we change the exponent to 4, we get a failure:
f = chebfun('-x -x.^2 + exp(-(30*(x-.5)).^4)');
length(f)
plot(f)
%%
% What has happened can be explained as follows. The function being sampled
% has a narrow spike near x = 0.5, and the closest grid points lie near
% 0.383 and 0.707. In the case of the exponent 2, we note that at x=0.383,
% exp(-(30(x-.5)^2))=4.5e-6, which is large enough to be noticed by the
% Chebfun constructor. On the other hand in the case of exponent 4, we have
% exp(-(30(x-.5)^4))=1.2e-66, which is far below machine precision. So in
% the latter case the constructor thinks it has a quadratic and does not
% try a finer grid.
%%
% If we increase minsamples, the correct chebfun is found:
chebfunpref('minsamples',17)
f = chebfun('-x -x.^2 + exp(-(30*(x-.5)).^4)');
length(f)
plot(f)
%%
% Incidentally, if the value of minsamples specified is not one greater
% than a power of 2, it is rounded up to the next such value. So
% chebfunpref('minsamples',10) would give the same result as
% chebfunpref('minsamples',17).
%%
% The default minsamples=9 was chosen as a good compromise between
% efficiency and reliability. In practice it rarely seems to fail, but
% perhaps it is most vulnerable when applied in splitting on mode to
% functions with narrow spikes. For example, the following chebfun is
% missing most of the spikes that should be there:
chebfunpref('factory')
ff = @(x) max(.8,sin(x+x.^2)) - x/20;
f = chebfun(ff,[0,10],'splitting','on');
plot(f)
%%
% Increasing minsamples fills them in:
chebfunpref('minsamples',17)
f = chebfun(ff,[0,10],'splitting','on');
plot(f)
%% 8.7 resampling: resampling on/off
% We now turn to a particularly interesting preference for Chebfun geeks,
% relating to the very idea of what it means to sample a function.
%%
% When a chebfun is constructed, a function is normally sampled at 9, 17,
% 33,... Chebyshev points until convergence is achieved. Now Chebyshev
% grids are nested, so the 17-point grid, for example, only contains 8
% points that are not in the 9-point grid. By default, the Chebfun
% constructor takes advantage of this property so as not to recompute
% values that have already been computed. (The default went the other way
% until 2009.)
%%
% For example, here is a chebfun constructed in the usual factory mode:
chebfunpref('factory');
ff = @(x) besselj(x,exp(x))
tic, f = chebfun(ff,[0 8]); toc
length(f)
%%
% Let us see what happens if we set 'resampling on', so that previously
% computed values are not reused:
tic, f = chebfun(ff,[0 8],'resampling','on'); toc
length(f)
%%
% One might wonder why 'resampling on' is an option at all, but in fact, it
% introduces some very interesting possibilities. What if the "function"
% being sampled is not actually a fixed function, but depends on the grid?
% For example, consider this prescription:
ff = @(x) length(x)*sin(2*x);
%%
% The values of f at any particular point will depend on the length of the
% vector in which it is embedded! What will happen if we try to make a
% chebfun, disabling the "sampletest" feature that is usually applied by
% the constructor as a safety test? The constructor tries the 9-point
% Chebyshev grid, then the 17-point grid, then the 33-point grid. On the
% last of these it finds the Chebyshev coefficients are sufficiently small,
% and proceeds to truncate to length 18. We end up with a chebfun of length
% 18 that precisely matches the function 33sin(2x).
f = chebfun(ff,'sampletest',0,'resampling','on');
length(f)
max(f)
plot(f,'.-')
%%
% This rather bizarre example encourages us to play further. What if we
% change length(x)*sin(2*x) to sin(length(x)*x)? Now there is no
% convergence, for no matter how fine the grid is, the function is
% underresolved.
hh = @(x) sin(length(x)*x);
h = chebfun(hh,'sampletest',0,'resampling','on');
%%
% Here is an in-between case where convergence is achieved on the grid of
% length 65, and the resulting chebfun then trimmed to length 44.
kk = @(x) sin(length(x).^(2/3)*x);
k = chebfun(kk,'sampletest',0,'resampling','on');
length(k)
plot(k,'.-')
%%
% Are such curious effects of any use? Yes indeed, they are at the heart
% of the chebop system. When the chebop system solves a differential
% equation by a command like u = L\f, for example, the chebfun u is
% determined by a "sampling" process in which a matrix problem obtained by
% Chebyshev spectral discretization is solved on grids of size 9, 17, and
% so on. The matrices change with the grids, so the sample values for u
% are crucially grid-dependent. Without resampling, chebops would not
% work.
%%
% Besides chebops, are there other practical uses of the Chebfun resampling
% feature? We do not currently know the answer and would be pleased to hear
% from users who may have ideas.
%% 8.8 eps: Chebfun constructor tolerance
% One of the controllable preferences is all too tempting: you can weaken
% the tolerance used in constructing a chebfun. The chebfunpref parameter
% eps is set by default to machine precision:
chebfunpref('eps')
%%
% However, one can change this with a command like chebfunpref('eps',1e-6).
%%
% There are cases where weakening the tolerance makes a big
% difference. For example, this happens in certain applications in 2D and
% in certain applications involving differential equations. (Indeed, the
% Chebfun differential equations commands have their own tolerance control
% strategies.) However, Chebfun does such a good job at resolving many
% functions that the eps-adjustment feature is not as useful as you might
% imagine, and we recommend that users not change eps unless they are
% having real problems with standard precision.
%% 8.9 Additional preferences
% Information about additional Chebfun preferences can be found by typing
% help chebfunpref. Here is a quick summary. In most cases various
% keywords are permitted such as 'on' or 1 or true, 'off' or 0 or false.
%%
% 'sampletest' controls whether a function is evaluated at an extra point
% as a safety check of convergence. With the default 'on' value, this test
% is indeed carried out.
%%
% 'blowup' relates to the construction of chebfuns that diverge to
% infinity, as described in Chapter 9. blowup=0 is used for no
% singularities, blowup=1 if for functions with poles (blowups with a
% negative integer power) and blowup=2 for functions with branch points
% (blowups with an arbitrary power).
%%
% 'chebkind' is set by default to 2, but can be changed to 1 to force
% Chebfun to use Chebyshev points of the first kind (i.e., roots of
% Chebyshev polynomials) rather than the second kind (extrema). This
% feature is experimental and may not do what you expect.
%%
% 'extrapolate" can be turned on to avoid evaluating functions at endpoints
% if this may lead to trouble with NaN or inf. For example, here we get a
% warning message because the value is 0 at the x=0 and 1 elsewhere,
f = chebfun('sign(x)',[0 1]);
%%
% whereas here there is no difficulty:
f = chebfun('sign(x)',[0 1],'extrapolate','on');
%%
% 'plot_numpts', set by default to 2001, controls how many points a
% function is sampled at for plotting.
%%
% 'ADdepth' stands for 'automatic differentiation depth', and limits how
% deep an automatic differentiation stack can go in order to limit memory
% use.