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