%% 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.