%% CHEBFUN GUIDE 9: INFINITE INTERVALS, INFINITE FUNCTION VALUES, AND SINGULARITIES
% Lloyd N. Trefethen, November 2009, revised February 2011

%%
% This chapter presents some features of
% Chebfun that are less robust than what is
% described in the first eight chapters.  With classic bounded
% chebfuns on a bounded interval [a,b], you can do amazingly
% complicated things often without encountering any difficulties.
% Now we are going to let the intervals and the functions
% diverge to infinity -- but please
% lower your expectations!  Partly because the software is relatively new,
% and partly for good mathematical reasons,
% one cannot expect the same reliability with these features.

%% 9.1 Infinite intervals
% If a function converges reasonably rapidly to a constant at infinity,
% you can define a corresponding chebfun.  Here are a couple of
% examples on [0,inf].  First we plot a function and find its maximum:
f = chebfun('0.75 + sin(10*x)./exp(x)',[0 inf]);
plot(f)
maxf = max(f)

%%
% Next we plot another function and integrate it from 0 to inf:
g = chebfun('1./(gamma(x+1))',[0 inf]);
sumg = sum(g)
plot(g,'r')

%%
% Where do f and g intersect?  We can find out using roots:
plot(f), hold on, plot(g,'r')
r = roots(f-g)
plot(r,f(r),'.k')

%%
% Here's an example on [-inf,inf] with a calculation of
% the location and value of the minimum:
g = chebfun(@(x) tanh(x-1),[-inf inf]);
g = abs(g-1/3);
clf, plot(g)
[minval,minpos] = min(g)

%%
% Notice that a function on an infinite domain is by default plotted
% on an interval like [0,10] or [-10,10].  You can use an extra
% 'interval' flag to plot
% on other intervals, as shown by this example of a function
% of small norm whose largest values are near x=30:
hh = @(x) cos(x)./(1e5+(x-30).^6);
h = chebfun(hh,[0 inf]);
plot(h,'interval',[0 100])
normh = norm(h)
%%
% Chebfun provides a convenient tool for the
% numerical evaluation of integrals over infinite domains:
g = chebfun('(2/sqrt(pi))*exp(-x.^2)',[0 inf])
sumg = sum(g)

%%
% The cumsum operator applied to this integrand gives us the error function,
% which matches Matlab's erf function well:
errorfun = cumsum(g)
disp('          erf               errorfun')
for n = 1:6, disp([erf(n) errorfun(n)]), end

%%
% One should be cautious in evaluating integrals over infinite intervals,
% however, for as mentioned in Section 1.5, the accuracy
% is sometimes disappointing,
% especially for functions that do not decay very quickly:
sum(chebfun('(1/pi)./(1+s.^2)',[-inf inf]))

%%
% Here's an example of a function that is too wiggly to be
% fully resolved:
sinc = chebfun('sin(pi*x)./(pi*x)',[-inf inf]);
plot(sinc,'m','interval',[-10 10])

%%
% Chebfun's capability of handling infinite intervals was introduced by
% Rodrigo Platte in 2008-09.  The basis of these computations is a change
% of variables, or mapping, which simplifies the infinite interval to
% [-1,1].  Let's take a look at what is going on in the case
% of the function g just constructed.  We'll do this by digging inside
% Chebfun a bit -- with a warning that the details
% of these inner workings are not fixed, but may change
% in future releases.

%%
% First we look at the different fields that make up a chebfun:
struct(g)

%%
% The first field, funs,
% indicates in this case that g consists of a single fun, that is, it
% is not split into pieces.  We can look at that fun like this:
f = g.funs;
struct(f)

%%
% The crucial item here is map, which in turn has a number of fields:
m = f.map

%%
% What's going on here is the use of a rational mapping from the
% y variable in [-1,1] to the x variable in [0,inf].
% The forward map is a Matlab anonymous function,
m.for

%%
% The inverse map is another anonymous function,
m.inv

%%
% The derivative of the forward map is used in the calculation
% of integrals and derivatives:
m.der

%%
% Matlab's anonymous functions don't make the values of the
% parameters a and s visible, though as it happens these numbers can be found
% as the first and last entries of the par field,
m.par

%%
% The use of mappings to simplify an unbounded domain to a bounded one
% is an idea that has been employed many times
% over the years.  One of the references we have benefitted especially from, which
% also contains pointers to other works in this area, is the book [Boyd 2001].

%% 9.2 Poles
% Chebfun can handle certain "vertical" as well as
% "horizontal" infinities -- especially, functions that blow up according
% to an integer power, i.e., with a pole.  If you know the nature of the blowup,
% it is a good idea to specify it using the 'exps' flag.
% For example, here's a function with a pole at 0.  We can use
% 'exps' to tell the constructor that the function looks like x^(-1)
% at the left endpoint and x^0 (i.e., smooth) at the right
% endpoint.
f = chebfun('sin(50*x) + 1./x',[0 4],'exps',[-1,0]);
plot(f), ylim([-5 30])

%%
% Here's the same function but over a domain that contains the
% singularity in the middle.  We tell the constructor where the
% pole is and what the singularity is like:
f = chebfun('sin(50*x) + 1./x',[-2 0 4],'exps',[0,-1,0]);
plot(f), ylim([-30 30])

%%
% Here's the tangent function:
f = chebfun('tan(x)', pi*((-5/2):(5/2)), 'exps', -ones(1,6));
plot(f), ylim([-5 5])

%%
% Rootfinding works as expected:
x2 = chebfun('x/2',pi*(5/2)*[-1 1]);
hold on, plot(x2,'k')
r = roots(f-x2); plot(r,x2(r),'or','markersize',8)

%%
% And we can manipulate the function in various other familiar ways:
g = sin(2*x2)+min(abs(f+2),6);
hold off, plot(g)

%%
% If you don't know what singularities your function may
% have, Chebfun has some ability to find them if the flags
% "blowup" and "splitting" are on:
gam = chebfun('gamma(x)',[-4 4],'splitting','on','blowup',1);
plot(gam), ylim([-10 10])

%%
% But it's always better to specify the breakpoints and powers if you know them:
gam = chebfun('gamma(x)',[-4:0 4],'exps',[-1 -1 -1 -1 -1 0]);

%%
% If you know the breakpoints but not the strengths of the poles, you
% can specify NaN at locations of unknown strength:
%
f = chebfun(@(t) exp(t)./(exp(t)-1),[0 1],'exps',[NaN 0])

%%
% It's also possible to have poles of different strengths
% on two sides of a singularity.  In this case, you specify two
% exponents at each internal breakpoint rather than one:
f = chebfun(@(x) cos(100*x)+sin(x).^(-2+sign(x)),[-1 0 1],'exps',[0 -3 -1 0]);
plot(f), ylim([-30 30])


%% 9.3  Singularities other than poles
% Less reliable but also sometimes useful is the possibility of working
% with functions with algebraic singularities that are not poles.
% Here's a function with inverse square root singularities at each end:
w = chebfun('(2/pi)./(sqrt(1-x.^2))','exps',[-.5 -.5]);
plot(w,'m'), ylim([0 10])

%%
% The integral is 2:
sum(w)

%%
% We pick this example because
% Chebyshev polynomials are the orthogonal polynomials with respect
% to this weight function, and Chebyshev coefficients are defined by
% inner products against Chebyshev polynomials with respect to this
% weight.  For example, here we compute inner products of x^4 + x^5
% against the Chebyshev polynomials T0,...,T5.  (The integrals
% in these inner products
% are calculated by Gauss-Jacobi quadrature using methods implemented by
% Nick Hale; for more on this subject see the command jacpts.)
x = chebfun('x');
T = chebpoly(0:5)';
f = x.^4 + x.^5;
chebcoeffs1 = T*(w.*f)

%%
% Here for comparison are the Chebyshev coefficients as obtained
% from chebpoly:
chebcoeffs2 = flipud(chebpoly(f)')

%%
% Notice the excellent agreement except for coefficient a0.  As mentioned
% in Section 4.1, in this special case
% the result from the inner product must be multiplied by 1/2.

%% 
% You can specify singularities for functions that don't blow up, too.
% For example, suppose we want to work with sqrt(x*exp(x)) on the interval
% [0,2].  A first try fails completely:
ff = @(x) sqrt(x.*exp(x));
d = domain(0,2);
f = chebfun(ff,d)

%%
% We could turn splitting on and resolve the function by many
% pieces, as illustrated in Section 8.3:
f = chebfun(ff,d,'splitting','on')

%%
% A much better representation, however, is constructed if we
% tell Chebfun about the singularity at x=0:
f = chebfun(ff,d,'exps',[.5 0])
plot(f)

%%
% Under certain circumstances Chebfun will introduce singularities
% like this of its own accord.  For example, just as abs(f) introduces
% breakpoints at roots of f, sqrt(abs(f)) introduces breakpoints and
% also singularities at such roots:
theta = chebfun('t',[0,4*pi]);
f = sqrt(abs(sin(theta)))
plot(f)
sumf = sum(f)

%%
% If you have a function that blows up but you don't know the
% nature of the singularities, even whether they are poles or
% not, Chebfun will try to figure them
% out automatically if you run in 'blowup 2' mode.  Here's an example

%%
f = chebfun('x.*(1+x).^(-exp(1)).*(1-x).^(-pi)','blowup',2)

%%
% Notice that the 'exps' field shows values close
% to -e and -pi, as is confirmed by looking
% at the numbers to higher precision:
f.funs.exps

%%
% The treatment of blowups in Chebfun
% was initiated by Mark Richardson in an MSc thesis at
% Oxford [Richardson 2009], then further developed by
% Richardson in collaboration with Rodrigo Platte and Nick Hale.  

%% 9.4 Another approach to singularities
% We have just seen how certain algebraic singularities can be
% captured in Chebfun via the 'exps' flag, which represents the
% function in question by a polynomial times certain algebraic terms.
% One example looked like this:
ff = @(x) sqrt(x.*exp(x));
d = domain(0,2);
f = chebfun(ff,d,'exps',[.5 0])

%%
% It is easy to run into trouble with such calculations, however, for
% example if we define
gg = @(x) 1+sqrt(x.*exp(x));

%%
% and then try g = chebfun(gg,d,'exps',[.5 0]).  
% Better results may sometimes be achieved -- for bounded functions
% only -- by using a different
% experimental Chebfun facility based on the same kind of changes of
% variable that make infinite intervals possible as described in
% Section 9.1.  For example, we may write
g = chebfun(gg,d,'singmap',[.5 0])

%%
% Here's another example:
ff = @(x) cos(22*x)./(1+x.^2)  + gamma(1.2+x).*real(airy(10*x)).*(1-x).^.3;
f1 = chebfun(ff,'singmap',[0 .3])
plot(min(f1,1))

%% 9.5 References
%
% [Boyd 2001] J. P. Boyd, Chebyshev and Fourier Spectral Methods,
% 2nd ed., Dover, 2001.
%
% [Richardson 2009] M. Richardson, Approximating Divergent
% Functions in the Chebfun System, thesis, MSc
% in Mathematical Modelling and Scientific Computing,
% Oxford University, 2009.
%