%% Lee and Greengard ODE examples % Nick Trefethen, 12th June 2012 %% % (Chebfun Example ode/LeeGreengardODEs.m) % [Tags: #linearODE, #boundarylayer] %% % In 1997 Lee and Greengard published a paper called "A fast adaptive % numerical method for stiff two-point boundary value problems" [1]. The % algorithm described there, being adaptive, can handle far stiffer % problems than Chebfun. Nevertheless Chebfun does pretty well with Lee % and Greengard's interesting collection of examples. These problems are % linear. %% EXAMPLE 1. VISCOUS SHOCK % The first example is % $$\varepsilon u''(x) + 2 x u'(x) = 0, \quad u(-1) = 1, ~ u(1) = 1.$$ % The following anonymous function produces a Chebfun solution as a % function of ep: uep = @(ep) chebop(@(x,u) ep*diff(u,2) + 2*x.*diff(u),-1,1)\0; %% % Here we plot the solution for ep = 0.01, 0.0001. It works fine, but Lee % and Greengard can go down to 1e-14. FS = 'fontsize'; LW = 'linewidth'; for k = 1:2 ep = 10^(-2*k); subplot(2,1,k) tic, u = uep(ep); t = toc; plot(u,'m',LW,2), ylim(1.4*[-1 1]); grid on s = sprintf('Ep = %5.1e Length =%4d Time =%6.3f',ep,length(u),t); title(s,FS,10) end %% EXAMPLE 2. BESSEL EQUATION % The second example is % $$u''(x) + x^{-1} u'(x) + {x^2-\nu^2\over x^2} u(x) = 0, % \quad u(0) = 0, ~ u(600) = 1$$ % with $\nu = 100$. We follow the same pattern as before (multiplying % the equation through by $x^2$): unu = @(nu) chebop(@(x,u) x.^2.*diff(u,2)+x.*diff(u)+(x.^2-nu^2).*u,[0,600],0,1)\0; %% % The solution for nu = 100 looks just as in Lee and Greengard: nu = 100; tic, u = unu(nu); t = toc; clf, plot(u,LW,1.2), grid on s = sprintf('nu = %3d Length =%4d Time =%6.3f',nu,length(u),t); title(s,FS,10) %% EXAMPLE 3. TURNING POINT % This example is an Airy equation, % $$\varepsilon u''(x) - x u(x) = 0, \quad u(-1) = 1, ~ u(1) = 1.$$ % We proceed as usual: uep = @(ep) chebop(@(x,u) ep*diff(u,2)-x.*u,1,1)\0; clf for k = 1:2 ep = 10^(-3*k); subplot(2,1,k) tic, u = uep(ep); t = toc; plot(u,'r',LW,1.2), grid on s = sprintf('Ep = %5.1e Length =%4d Time =%6.3f',ep,length(u),t); title(s,FS,10) end %% EXAMPLE 4. POTENTIAL BARRIER % The fourth example is % $$\varepsilon u''(x) + (x^2-0.25)u(x) = 0, \quad u(-1) = 1, ~ u(1) = 2.$$ % Here we go: uep = @(ep) chebop(@(x,u) ep*diff(u,2)+(x.^2-0.25).*u,1,2)\0; for k = 1:2 ep = 10^(-3*k); subplot(2,1,k) tic, u = uep(ep); t = toc; plot(u,'color',[0 .7 0],LW,1.2), grid on s = sprintf('Ep = %5.1e Length =%4d Time =%6.3f',ep,length(u),t); title(s,FS,10) end %% EXAMPLE 5. CUSP % This time we have % $$\varepsilon u''(x) + x u'(x) - 0.5u(x) = 0, \quad u(-1) = 1, ~ u(1) = 2.$$ % With a global discretization and standard defaults, Chebfun can go down % to 1e-5 or so. With a breakpoint introduced at $x=0$ by specifying % the domain [-1 0 1], we get a little further, though not as far % as Lee and Greengard: uep = @(ep) chebop(@(x,u) ep*diff(u,2)+x.*diff(u)-0.5*u,[-1 0 1],1,2)\0; for k = 1:2 ep = 10^(-3*k); subplot(2,1,k) tic, u = uep(ep); t = toc; plot(u,'color',[1 .5 .5],LW,2), grid on s = sprintf('Ep = %5.1e Length =%4d Time =%6.3f',ep,length(u),t); title(s,FS,10) end %% EXAMPLE 6. EXPONENTIAL ILL-CONDITIONING % Finally we consider % $$\varepsilon u''(x) - x u'(x) + u(x) = 0, \quad u(-1) = 1, ~ u(1) = 2.$$ % The pictures look fine down to Lee and Greengard's value of ep = 1/70: uep = @(ep) chebop(@(x,u) ep*diff(u,2)-x.*diff(u)+u,[-1 0 1],1,2)\0; for k = 1:2 ep = (1/35)/k; subplot(2,1,k) tic, u = uep(ep); t = toc; plot(u,'color',[0 .8 .8],LW,2), grid on s = sprintf('Ep = %5.1e Length =%4d Time =%6.3f',ep,length(u),t); title(s,FS,10) end %% % However, this problem is highly ill-conditioned and I have not % investigated how accurate the solution really is. %% % REFERENCE % % [1] J.-Y. Lee and L. Greengard, "A fast adaptive numerical method for % stiff two-point boundary value problems", SIAM Journal on Scientific % Computing 18 (1997), 403-429.