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