# Lee and Greengard ODE examples

Nick Trefethen, 12th June 2012

## Contents

(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

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND =  1.612982e-20.


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.