%% CHEBFUN GUIDE 3: ROOTFINDING AND MINIMA AND MAXIMA
% Lloyd N. Trefethen, October 2009, revised February 2011
%% 3.1 roots
% Chebfun comes with a global rootfinding capability -- the
% ability to find all the zeros of a function in its region
% of definition. For example, here is a polynomial with two
% roots in [-1,1]:
x = chebfun('x');
p = x.^3 + x.^2 - x;
r = roots(p)
%%
% We can plot p and its roots like this:
plot(p)
hold on, plot(r,p(r),'.r')
%%
% Of course one does not need Chebfun to find roots of a polynomial:
roots([1 1 -1 0])
%%
% A more substantial example of rootfinding involving a Bessel function
% was considered in Sections 1.2 and 2.4. Here is a similar calculation for
% the Airy functions Ai and Bi, modeled after the page on "Airy functions"
% at WolframMathWorld. (The reason for the "real" command is a
% bug in Matlab: Matlab's "airy" command gives spurious nonzero imaginary
% parts, at the level of rounding errors, when applied to some real
% arguments, such as x=-2.)
Ai = chebfun('real(airy(0,x))',[-10,3]);
Bi = chebfun('real(airy(2,x))',[-10,3]);
hold off, plot(Ai,'r')
hold on, plot(Bi,'b')
rA = roots(Ai); plot(rA,Ai(rA),'.r')
rB = roots(Bi); plot(rB,Bi(rB),'.b')
axis([-10 3 -.6 1.5]), grid on
%%
% Here for example are the three roots of Ai and Bi closest to 0:
[rA(end-2:end) rB(end-2:end)]
%%
% Chebfun finds roots by a method due to Boyd and Battles
% [Boyd 2002, Battles 2006]. If the chebfun is of degree greater than 100,
% it is broken into smaller pieces recursively. On each small piece
% zeros are then found as eigenvalues of a "colleague matrix", the analogue
% for Chebyshev polynomials of a companion matrix for monomials [Specht 1960, Good 1961].
% This method can be startlingly fast and accurate. For example,
% here is a sine function with 11 zeros:
f = chebfun('sin(pi*x)',[0 10]);
lengthf = length(f)
tic, r = roots(f); toc
r
%%
% A similar computation with 101 zeros comes out equally well:
f = chebfun('sin(pi*x)',[0 100]);
lengthf = length(f)
tic, r = roots(f); toc
fprintf('%22.14f\n',r(end-4:end))
%%
% And here is the same on an interval with 1001 zeros.
f = chebfun('sin(pi*x)',[0 1000]);
lengthf = length(f)
tic, r = roots(f); toc
fprintf('%22.13f\n',r(end-4:end))
%%
% With the ability to find zeros, we can solve a variety of
% nonlinear problems. For example, where do the curves
% x and cos(x) intersect? Here is the answer.
x = chebfun('x',[-2 2]);
hold off, plot(x)
f = cos(x);
hold on, plot(f,'k')
r = roots(f-x)
plot(r,f(r),'or')
%%
% All of the examples above concern chebfuns consisting of
% a single fun. If there are several funs, then roots are
% included at jumps as necessary. The following
% example shows a chebfun with 26 pieces. It has
% nine zeros: one within a fun, eight at jumps between funs.
x = chebfun('x',[-2 2]);
f = x.^3 - 3*x - 2 + sign(sin(20*x));
hold off, plot(f), grid on
r = roots(f);
hold on, plot(r,0*r,'.r')
%% 3.2 min, max, abs, sign, round, floor, ceil
% Rootfinding is more central to Chebfun than
% one might at first imagine, because a number
% of commands, when applied to smooth chebfuns, must
% produce non-smooth results, and it is rootfinding that
% tells us where to put the discontinuities.
% For example, the "abs" command introduces breakpoints
% wherever the argument goes through zero. Here we
% see that x consists of a single piece, whereas abs(x) consists
% of two pieces.
x = chebfun('x')
absx = abs(x)
subplot(1,2,1), plot(x,'.-')
subplot(1,2,2), plot(absx,'.-')
%%
% We saw this effect already in Section 1.4.
% Another similar effect shown in that section occurs with
% min(f,g) or max(f,g). Here, breakpoints are introduced
% at points where f-g is zero:
f = min(x,-x/2), subplot(1,2,1), plot(f,'.-')
g = max(.6,1-x.^2), subplot(1,2,2), plot(g,'.-'), ylim([.5,1])
%%
% The function "sign" also introduces breaks, as illustrated
% in the last section.
% The commands "round", "floor", and "ceil" behave like this too.
% For example, here is exp(x) rounded to nearest multiples of 0.5:
g = exp(x);
clf, plot(g)
gh = 0.5*round(2*g);
hold on, plot(gh,'r');
grid on
%% 3.3 Local extrema: roots(diff(f))
% Local extrema of smooth functions can be located by finding
% zeros of the derivative. For example, here is a variant of
% the Airy function again, with all its extrema in its range of
% definition located and plotted.
f = chebfun('exp(real(airy(x)))',[-15,0]);
clf, plot(f)
r = roots(diff(f));
hold on, plot(r,f(r),'.r'), grid on
%%
% This method will find non-smooth extrema as well as smooth ones,
% since these correspond to "zeros" of the derivative where
% the derivative jumps from one sign to the other. Here is an
% example.
x = chebfun('x');
f = exp(x).*sin(30*x);
g = 2-6*x.^2;
h = max(f,g);
clf, plot(h)
%%
% Here are all the local extrema, smooth and nonsmooth:
extrema = roots(diff(h));
hold on, plot(extrema,h(extrema),'.r')
%%
% Suppose we want to pick out the local extrema that are actually
% local minima. We can do that by checking for the second
% derivative to be positive:
h2 = diff(h,2);
maxima = extrema(h2(extrema)>0);
plot(maxima,h(maxima),'ok','markersize',12)
%% 3.4 Global extrema: max and min
% If "min" or "max" is applied to a single chebfun, it returns
% its global minimum or maximum. For example:
f = chebfun('1-x.^2/2');
[min(f) max(f)]
%%
% Chebfun computes such a result by checking the values of
% f at all endpoints and at zeros of the derivative.
%%
% As with standard Matlab, one can find the location of the
% extreme point by supplying two output arguments:
[minval,minpos] = min(f)
%%
% Note that just one position is returned even though the
% minimum is attained at two points. This is consistent
% with the behavior of standard Matlab.
%%
% This ability to do global 1D optimization in Chebfun
% is rather remarkable. Here is a nontrivial example.
f = chebfun('sin(x)+sin(x.^2)',[0,15]);
hold off, plot(f,'k')
%%
% The length of this chebfun is not as great as one might
% imagine:
length(f)
%%
% Here are its global minimum and maximum:
[minval,minpos] = min(f)
[maxval,maxpos] = max(f)
hold on
plot(minpos,minval,'.b','markersize',30)
plot(maxpos,maxval,'.r','markersize',30)
%%
% For larger chebfuns, it is inefficient to compute
% the global minimum and maximum separately like this -- each
% one must compute the derivative and find all its zeros.
% An alternative code is offered that computes both at once:
[extemevalues,extremepositions] = minandmax(f)
%% 3.5 norm(f,1) and norm(f,inf)
% The default, 2-norm form of the "norm" command was considered in Section 2.2.
% In standard Matlab one can also compute 1-, infinity-, and Frobenius norms
% with norm(f,1), norm(f,inf), and norm(f,'fro'). These
% have been overloaded in Chebfun, and in the first two cases, rootfinding
% is part of the implementation. (The 2- and Frobenius norms are
% equal for a single chebfun but differ for quasimatrices; see Chapter 6.)
% The 1-norm norm(f,1) is the integral of
% the absolute value, and Chebfun computes this by adding up segments
% between zeros, at which abs(f) will typically have a discontinuous slope.
% The infinity-norm is computed from the formula
% norm(f,inf) = max(max(f),-min(f));
%%
% For example:
f = chebfun('sin(x)',[103 103+4*pi]);
norm(f,inf)
norm(f,1)
%% 3.6 Roots in the complex plane
% Chebfuns live on real intervals, and the funs from which
% they are made live on real subintervals. But a polynomial
% representing a fun may have roots outside the interval of definition,
% which may be complex. Sometimes we may want to get our hands on these
% roots, and the "roots" command makes this possible in various
% ways through the flags 'all', 'complex', and 'norecurse'.
%%
% The simplest example is a chebfun that is truly intended to
% correspond to a polynomial. For example, the chebfun
f = 1+16*x.^2;
%%
% has no roots in [-1,1]:
roots(f)
%%
% We can extract its complex roots with the command
roots(f,'all')
%%
% The situation for more general chebfuns is more complicated. For example,
% the chebfun
g = exp(x).*f(x);
%%
% also has no roots in [-1,1],
roots(g)
%%
% but it has plenty of roots in the complex plane:
roots(g,'all')
%%
% Most of these are spurious.
% What has happened is that g is represented by a polynomial chosen
% for its approximation properties on [-1,1].
% Inevitably that polynomial will have roots
% in the complex plane, even if they have little to do with g.
%%
% One cannot expect Chebfun to solve this problem perfectly -- after all,
% it is working on a real interval, not in the complex plane, and analytic continuation
% from the one to the other is
% well known to be an ill-posed problem. Nevertheless, Chebfun may do a pretty good job of
% selecting genuine complex (and real)
% roots near the interval of definition if you use the 'complex' flag:
roots(g,'complex')
%%
% We will not go into detail here of how this is done, but the idea is that associated
% with any fun is a family of "Chebfun ellipses" in the complex plane, with foci at
% the endpoints, inside which
% one may expect reasonably good accuracy of the fun. Assuming the interval is
% [-1,1] and the fun has length L, the Chebfun ellipse associated with a parameter
% delta<<1 is the region in the complex plane bounded by the image under the map
% (z+1/z)/2 of the circle |z|=r, where r is defined by r^(-L)=delta. The
% command roots(g,'complex') first effectively does roots(g,'all'), then returns only
% those roots lying inside a particular Chebfun ellipse -- we take the one corresponding
% to delta equal to the square root of the Chebfun tolerance parameter eps.
%%
% One must expect complex roots of chebfuns to lose accuracy as one moves away from
% the interval of definition. Here's an example:
f = exp(exp(x)).*(x-0.1i).*(x-.3i).*(x-.6i).*(x-1i);
roots(f,'complex')
%%
% Notice that the first three imaginary roots are located with
% about 10, 8, and 6 digits of accuracy, while the fourth does not appear
% in the list at all.
%%
% Here is a more complicated example:
F = @(x) 4+sin(x)+sin(sqrt(2)*x)+sin(pi*x);
f = chebfun(F,[-100,100]);
%%
% This function has a lot of complex roots lying in strips on either side
% of the real axis.
r = roots(f,'complex');
hold off, plot(r,'.','markersize',6)
%%
% If you are dealing with complex roots of complicated chebfuns like this,
% it may be safer to add a flag 'nonrecurse' to the roots call, at the
% cost of slowing down the computation.
% This turns off the Boyd-Battles recursion mentioned
% above, lowering the chance of missing a few roots near interfaces of the
% recursion. If we try that here we find that the results look almost the
% same as before in a plot:
r2 = roots(f,'complex','norecurse');
hold on, plot(r,'om','markersize',8)
%%
% However, the accuracy has improved:
norm(F(r))
norm(F(r2))
%%
% To find poles in the complex plane as opposed to zeros, see Section 4.8.
%% 3.7 References
%
% [Battles 2006] Z. Battles, Numerical Linear Algebra for
% Continuous Functions, DPhil thesis, Oxford University
% Computing Laboratory, 2006.
%
% [Boyd 2002] J. A. Boyd, "Computing zeros on a real interval
% through Chebyshev expansion and polynomial rootfinding",
% SIAM Journal on Numerical Analysis 40 (2002), 1666-1682.
%
% [Good 1961] I. J. Good, "The colleague matrix, a Chebyshev
% analogue of the companion matrix", Quarterly Journal of
% Mathematics 12 (1961), 61-68.
%
% [Specht 1960] W. Specht, "Die Lage der Nullstellen eines Polynoms.
% IV", Mathematische Nachrichten 21 (1960), 201-222.