Chebfun Logo
Oxford University
Mathematical Institute

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.


HTML Comment Box is loading comments...

Please contact us with any questions and comments.
Copyright © 2013, The University of Oxford & The Chebfun Team.