Code covered by the BSD License  

Highlights from
Chebfun V4

image thumbnail

Chebfun V4

by

 

30 Apr 2009 (Updated )

Numerical computation with functions instead of numbers.

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

A parameter dependent ODE with breakpoints

A parameter dependent ODE with breakpoints

Asgeir Birkisson, 25th January 2012

Contents

(Chebfun example ode/ParameterODE.m)

(This example is inspired by a discussion with Paul Constantine [1]).

Let the ODE

(a(x,s)u')' = 1,   u(0) = u(1) = 0,

be given, where

a(x,s) = 1+4s(x^2-x).

and ' denotes the derivative w.r.t. x.

The exact solution can be shown to be

u(x,s) = 1/(8s) log(1+4s(x^2-x)) = 1/(8s) log(a(x,s))

It is clear that for s = 1, the solution has a singularity at x = 1/2. Here, we explore what happens when we solve the problems for values of s getting closer and closer to the critical value s = 1.

Setting up the problem

We begin by rewriting the differential equation to the form

a(x,s)u'' + a'(x,s)u' = 1

as it will be simpler to work with. We now set up anonymous functions to represent a and a'

a =  @(x,s)  1 + 4*s*(x.^2-x);
ap = @(x,s)  4*s*(2*x-1);

as well as an anonymous function for the exact solution, and the chebfun x on the interval [0,1]:

uexact = @(x,s) log(a(x,s)) / (8*s);
chebx = chebfun('x',[0 1]);

We can now set up a chebop to represent the BVP operator. However, since we want to explore what the solution looks like for different values of s, we define the chebop as an anonymous function as well (which output will be a chebop). The two last argument correspond to imposing homogenous Dirichlet conditions on the solution

Ns = @(s) chebop(@(x,u) a(x,s).*diff(u,2) + ap(x,s).*diff(u),[0 1], 0, 0);

Since we want to be taking values of s closer and closer to 1, we rewrite s in the following form:

s = 1-10^{-\gamma}

where \gamma takes integer values (giving s = 0.9, 0.99, 0.999, ...). We thus define s as an anonymous function

s = @(gamma)  1-10^(-gamma);

We can then obtain the solution of the problem for different values of gamma. Again, we use anonymous functions to achieve the desired effect.

ugamma = @(gamma) solvebvp(Ns(s(gamma)),1);

Here, the solvebvp method is another way to call the chebop backslash method. The second arguments corresponds to the right-hand side of the differential equation.

Solutions for different values of gamma

We're now all set to solve the problem for different values of gamma. We store the solution in a quasimatrix, as well as the norm of the residuals and the true errors:

gammas = 1:3;
sols = chebfun;
res = []; trueerror = [];
for gamma = gammas
    solgamma = ugamma(gamma);
    plot(solgamma), title(['Soln for gamma = ', num2str(gamma)]), drawnow
    sols(:,gamma) = solgamma;
    res(gamma) = norm(feval(Ns(s(gamma)),solgamma)-1);
    trueerror(gamma) = norm(solgamma - uexact(chebx,s(gamma)));
end

Here we are required to use the feval method to evaluate the residual since Matlab doesn't allowing double indexing, i.e. we can't call Ns(s(gamma))(solgamma).

Values of gamma up to 3 work fine, however, for gamma = 4, the solver fails to converge using 1025 points. We can see the lengths of the columns of the quasimatrix by displaying it:

sols
sols = 
   chebfun column 1 (1 smooth piece)
       interval       length   endpoint values   
[       0,       1]       97 -6.5e-14  1.3e-17   
vertical scale = 0.35 
   chebfun column 2 (1 smooth piece)
       interval       length   endpoint values   
[       0,       1]      303  3.4e-12        0   
vertical scale = 1.8 
   chebfun column 3 (1 smooth piece)
       interval       length   endpoint values   
[       0,       1]      387 -2.5e-11  6.3e-17   
vertical scale =  16 

We plot the solutions together in a single plot

figure, plot(sols)
title('u(x,s(\gamma))')
legend('\gamma = 1', '\gamma = 2', '\gamma = 3','Location','SE')
xlabel('x')

Looking at the entries in the vector storing the values of the residual reveals that they grow extremely fast with the value of gamma

figure, semilogy(gammas,res,'-*')
title('Norm of residual'), xlabel('\gamma')

This is due to large contributions to the residual close to the boundaries of the interval, as the following plot reveals

figure,semilogy(abs(feval(Ns(s(3)),sols(:,3))-1))

However, the true error remains much better under control:

figure, semilogy(gammas,trueerror,'-*')
title('Norm of true error'), xlabel('\gamma')

Introducing a breakpoint

The plot above of the solutions for different values of gamma reveals that the solution gets more and more difficult to represent close to x = 1/2 as gamma increases (i.e. s gets closer to 1). That makes a good case for introducing a breakpoint in the solution at x = 1/2, so rather than the solution being represented by a global chebfun, it is represented by two pieces.

We introduce a breakpoint in the operator as follows (notice the second argument to the chebop constructor):

Nsbreak = @(s) chebop(@(x,u) a(x,s).*diff(u,2) + ap(x,s).*diff(u), ...
                      [0 .5 1], 0, 0);

We now redefine the anonymous function which gives the solution

ugammabreak = @(gamma) solvebvp(Nsbreak(s(gamma)),1);

We're now all set to solve the problem using breakpoints for different values of gamma. Here, values of gamma up to 6 work with the default chebop settings, but we can obtain solutions for even larger values of gamma, up to 7, by increasing the number of gridpoints the solver can use (beyond 7, we run into problems with matrices being close to singular)

cheboppref('maxdegree',2^12+1)
chebx = chebfun('x',[0 .5 1]);
gammas = 1:7;
sols = chebfun;
res = []; trueerror = []; legs = [];
for gamma = gammas
    solgamma = ugammabreak(gamma);
    plot(solgamma), title(['Soln for gamma = ', num2str(gamma)]), drawnow
    sols(:,gamma) = solgamma;
    res(gamma) = norm(feval(Nsbreak(s(gamma)),solgamma)-1);
    trueerror(gamma) = norm(solgamma - uexact(chebx,s(gamma)));
    legs = [legs; {['\gamma = ' num2str(gamma)]}]; % Used for legend later
end

The length of the chebfuns are more modest now than before without breakpoints:

sols
sols = 
   chebfun column 1 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]       46 -1.4e-12    -0.32   
[     0.5,       1]       39    -0.32        0   
Total length = 85   vertical scale = 0.32 
   chebfun column 2 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]       90 -1.2e-12    -0.58   
[     0.5,       1]       71    -0.58  9.6e-18   
Total length = 161   vertical scale = 0.58 
   chebfun column 3 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]      126  4.1e-12    -0.86   
[     0.5,       1]      126    -0.86 -1.1e-16   
Total length = 252   vertical scale = 0.94 
   chebfun column 4 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]      110 -9.4e-14     -1.2   
[     0.5,       1]      110     -1.2  1.1e-16   
Total length = 220   vertical scale = 1.5 
   chebfun column 5 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]      186 -1.4e-12     -1.4   
[     0.5,       1]      186     -1.4        0   
Total length = 372   vertical scale = 1.7 
   chebfun column 6 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]      253  2.4e-12     -1.7   
[     0.5,       1]      253     -1.7 -5.6e-17   
Total length = 506   vertical scale =   2 
   chebfun column 7 (2 smooth pieces)
       interval       length   endpoint values   
[       0,     0.5]      443  1.2e-11       -2   
[     0.5,       1]      443       -2        0   
Total length = 886   vertical scale = 2.3 

We plot all the solutions in a single plot

figure, plot(sols)
title('u(x,s(\gamma))')
legend(legs,'Location','SE')
xlabel('x')

Again, the true error of the solutions remain pretty well under control

figure, semilogy(gammas,trueerror,'-*')
title('Norm of true error'), xlabel('\gamma')

(But the residual norm gives infinity due to the jump in the residual at the breakpoint).

Reference [1] Paul Constantine's website: http://www.stanford.edu/~paulcon/

Contact us