Spectral solution of 2-point boundary value problems using Chebfun¶

Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington.

See http://faculty.washington.edu/rjl/classes/am590a2013/codes.html for more IPython Notebook examples.

This notebook demonstrates using the Matlab Chebfun package to solve 2-point boundary value problems via a Chebyshev spectral method, following the examples in Chapter 21 of

[ATAP] L.N. Trefethen, "Approximation Theory and Approximation Practice", http://www.ec-securehost.com/SIAM/OT128.html

For more about how to do this, see

http://www.mathworks.com/matlabcentral/fileexchange/23972-chebfun/content/chebfun/guide/html/guide7.html

Using Matlab and Chebfun from IPython...¶

First we need to set things up to use Matlab and set the path to find chebfun.

You must have pymatbridge installed and chebfun available for this to work.

In [111]:

import pymatbridge
ip = get_ipython()
pymatbridge.load_ipython_extension(ip)

In [112]:

%%matlab
path(path,'/Users/rjl/chebfun/chebfun_v4.2.2889/chebfun')

Setting up the problem¶

Define the chebfun x and the linear operator $\partial_x^2 + I$ with Dirichlet boundary conditions:

In [113]:

%%matlab
x = chebfun('x');
L = chebop(@(u) diff(u,2) +  u)

L = chebop
   Linear operator operating on chebfuns defined on:          
     interval [-1,1]                                          
   representing the operator:                                 
     @(u)diff(u,2)+u                                          
   with n = 6 realization:                                    
     42.6000  -68.3607   40.8276  -23.6393   17.5724   -8.0000
     21.2859  -30.5331   12.6833   -3.6944    2.2111   -0.9528
     -1.8472    7.3167   -9.0669    5.7889   -1.9056    0.7141
      0.7141   -1.9056    5.7889   -9.0669    7.3167   -1.8472
     -0.9528    2.2111   -3.6944   12.6833  -30.5331   21.2859
     -8.0000   17.5724  -23.6393   40.8276  -68.3607   42.6000

Now let's specify some boundary conditions:

In [114]:

%%matlab
uleft = 3.;
uright = 4;
L.lbc = uleft;
L.rbc = uright;

Note that this changes the $6\times 6$ matrix shown: the first 4 rows now map function values at 6 Chebyshev points to values of $Lu$ at 4 Chebyshev points...

In [115]:

%%matlab
L

L = chebop
   Linear operator operating on chebfuns defined on:          
     interval [-1,1]                                          
   representing the operator:                                 
     @(u)diff(u,2)+u                                          
   left boundary condition:                                   
     3                                                        
   right boundary condition:                                  
     4                                                        
   with n = 6 realization:                                    
     33.0353  -51.2657   27.8620  -14.2070   10.1514   -4.5760
     -0.6075    5.7509   -9.1443    6.9332   -3.2805    1.3482
      1.3482   -3.2805    6.9332   -9.1443    5.7509   -0.6075
     -4.5760   10.1514  -14.2070   27.8620  -51.2657   33.0353
      1.0000         0         0         0         0         0
           0         0         0         0         0    1.0000

Specify the right hand side in the ODE $u''(x) + u(x) = f(x)$ :

In [116]:

%%matlab
f = 0*x

f = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[      -1,       1]        1        0        0   
vertical scale =   0

Solve the equation using Chebfun:¶

In [117]:

%%matlab
g = L\f
plot(g,'r')
title('The Chebfun solution')

g = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[      -1,       1]       29        3        4   
vertical scale = 6.5

Compare to the exact solution:¶

In this case we know the exact solution is a linear combination of $\sin(x)$ and $\cos(x)$ and the boundary conditions give us a $2\times 2$ linear system to solve for the coefficients:

In [118]:

%%matlab
A = [cos(-1) sin(-1); cos(1) sin(1)];
b = [uleft; uright];
c = A\b;
u = c(1)*cos(x) + c(2)*sin(x);
error = max(abs(u-g))
plot(u)
title('The exact solution')

error =
     1.232273452396662e-11

Airy equation¶

Next we solve the problem in (21.4) of [ATAP]:

$u'' - xu = 0, \quad u(-30) = 1, ~u(30) = 0.$

In [119]:

%%matlab
L = chebop(@(x,u) diff(u,2) - x.*u, [-30,30]);
L.lbc = 1; L.rbc = 0; u = L\0;
plot(u)
u

u = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[     -30,      30]      184        1  2.2e-16   
vertical scale = 6.1

The nonlinear pendulum¶

Finally we solve (21.5a) and (21.5b), the nonlinear pendulum problem

$\theta'' + \sin(\theta) = 0, \quad \theta\in [0,6]$ with boundary conditions

$u(0) = -\pi/2, \quad u(6) = \pi/2.$

In [120]:

%%matlab
N = chebop(0,6);
N.op = @(theta) diff(theta,2) + sin(theta);

In [121]:

%%matlab
N.lbc = -pi/2;
N.rbc = pi/2;
theta = N\0

theta = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[       0,       6]       70     -1.6      1.6   
vertical scale = 5.4

In [122]:

%%matlab
plot(-cos(theta)), grid on, ylim([-1 1])
title('Nonlinear pendulum (21.5)','fontsize',20)
xlabel('t','fontsize',20), ylabel('height -cos(\theta)','fontsize',20)

With different boundary conditions:

In [123]:

%%matlab
N.lbc = -pi/2; N.rbc = 5*pi/2; theta = N\0;
plot(-cos(theta)), grid on, ylim([-1 1])
title('Nonlinear pendulum (21.5), another solution','fontsize',20)
xlabel('t','fontsize',20), ylabel('height -cos(\theta)','fontsize',20)

Note that theta increases past $\pi$ , indicating the pendulum goes past the vertical and then back down:

In [124]:

%%matlab
plot(theta)

Imposing boundary conditions¶

Solve

$u''(x) + u(x) = exp(x)$ with boundary conditions

$u(0)=3, u'(2) = 1$

In [125]:

%%matlab
L = chebop(@(u) diff(u,2) + u,[0,2],@(u) u-3,@(u) diff(u)-1)

L = chebop
   Linear operator operating on chebfuns defined on:          
     interval [0,2]                                           
   representing the operator:                                 
     @(u)diff(u,2)+u                                          
   left boundary condition:                                   
     @(u)u-3 = 0                                              
   right boundary condition:                                  
     @(u)diff(u)-1 = 0                                        
   with n = 6 realization:                                    
     33.0353  -51.2657   27.8620  -14.2070   10.1514   -4.5760
     -0.6075    5.7509   -9.1443    6.9332   -3.2805    1.3482
      1.3482   -3.2805    6.9332   -9.1443    5.7509   -0.6075
     -4.5760   10.1514  -14.2070   27.8620  -51.2657   33.0353
      1.0000         0         0         0         0         0
     -0.5000    1.1056   -1.5279    2.8944  -10.4721    8.5000

In [126]:

%%matlab
f = chebfun('exp(x)', [0,2]);
u = L\f

u = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[       0,       2]       23        3      3.6   
vertical scale = 3.6

In [127]:

%%matlab
plot(u)

Clearly $u(0) = 3$ . Check that $u'(2) = 1$ :

In [128]:

%%matlab
uprime = diff(u);
uprime(2)

ans =
   1.000000000001889

Plot residual $u''(x) - u(x) - e^x$ :

In [129]:

%%matlab
residual = diff(uprime) + u - f;
max_residual = max(abs(residual))
plot(residual)

max_residual =
     4.465903202799382e-10

In [129]: