Code covered by the BSD License

Chebfun V4

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

Exact solutions of some ODEs

Exact solutions of some ODEs

Nick Trefethen, December 2010

Contents

(Chebfun example ode/ExactSolns.m)

Chapter 1 of the textbook by Bender and Orszag contains an intense review of a number of methods for solving ODEs exactly. Here are some examples illustrating techniques presented in that chapter. In each case we solve an ODE with Chebfun and compare with the exact solution. For simplicity we pose all the equations on the domain [1,2]:

d = [1 2];
x = chebfun('x',d);
N = chebop(d);


EXAMPLE 1: SEPARATION OF VARIABLES

Consider first the problem

   xy' = y^2 - 2y + 1,    y(1) = 0.

We can separate variables to get

   dy/(1-y)^2 = dx/x,

which can be integrated to give the exact solution

   y = 1 - 1/(C+log(x))

for some constant C. For the given boundary condition the constant is C=1, so we have

exact = 1 - 1./(1+log(x));


A Chebfun solution goes like this:

N.op = @(y) x.*diff(y) - y.^2 + 2*y;
N.lbc = 0;
y = N\1;
err = norm(y-exact,inf);


Here is a plot of the solution, using '.-' to show the grid needed to resolve the solution.

LW = 'linewidth'; FS = 'fontsize'; MS = 'markersize';
plot(y,'.-',LW,1,MS,18), grid on
title(sprintf('xy'' = y^2-2y+1     Error = %6.2e',err),FS,14)


EXAMPLE 2: SEPARATION OF VARIABLES

As another example, consider

   y' = sin(y),   y(1) = pi/2.

Separating variables now gives

   dy/sin(y) = dx

which implies

   log(tan(y/2)) = x + C

   y = 2 atan(Cexp(x))

for some constant C. We can satisfy the boundary condition with C=exp(-1), so the exact solution is

exact = 2*atan(exp(x-1));


A Chebfun solution is as follows:

N.op = @(y) diff(y) - sin(y);
N.lbc = pi/2;
y = N\0;
err = norm(y-exact,inf);


Here is a plot:

plot(y,'.-',LW,1,MS,18), grid on
title(sprintf('y'' = sin(y)     Error = %6.2e',err),FS,14)


EXAMPLE 3: ORDER REDUCTION FOR AUTONOMOUS EQUATION

Consider the autonomous problem

   yy' = 2(y')^2,   y(1) = 1, y(2) = 2,

where we think of y as a function of x, i.e., y' = dy/dx. Now introduce the new variable u = y', which we think of as a function of y, i.e., u' = du/dy. The equation becomes a new equation in u and y of first order,

   yuu' = 2u^2.

Separating variables gives u'/u = 2/y, which leads to log(u) = 2log(y)+C for some constant C, i.e., u = Cy^2, i.e.

   y' = Cy^2.

Separating variables again and integrating gives

   y = 1/(Cx+D)

for constants C and D. We can satisfy the boundary conditions with D=3/2, C=-1/2, giving the exact solution

exact = 2./(3-x);


A Chebfun solution is as follows:

N.op = @(y) y.*diff(y,2) - 2*diff(y).^2;
N.lbc = 1; N.rbc = 2;
y = N\0;
err = norm(y-exact,inf);


Here is a plot:

plot(y,'.-',LW,1,MS,18), grid on
title(sprintf('yy'''' = 2(y'')^2     Error = %6.2e',err),FS,14)


EXAMPLE 4: BERNOULLI EQUATION

An equation of the form y' = a(x)y +b(x)y^p, known as a Bernoulli equation, can be made linear by the change of variables u=y^(1-p). Consider for example the problem

   y' = y/x + x/y,  y(1) = 1.

Setting u=y^2 reduces the problem to

   u'/2 = u/x + x.

Multiplying by the integrating factor x^(-2) gives

   u'/x^2 - 2u/x^3 = 2/x,

that is,

   (u/x^2)' = 2/x,

which can be integrated to give

   y = xsqrt(C+2log(x)).

For our boundary conditions the solution is accordingly

exact = x.*sqrt(1+2*log(x));


For a Chebfun solution, we get an error message if we specify no initial guess as the default guess of the zero function causes division by zero. Accordingly here is a solution in which an initial guess is specified.

N.op = @(y) diff(y) - y./x - x./y;
N.lbc = 1; N.rbc = [];
N.init = 1;
y = N\0;
err = norm(y-exact,inf);


And here's the plot:

plot(y,'.-',LW,1,MS,18), grid on
title(sprintf('y'' = y/x + x/y     Error = %6.2e',err),FS,14)


Reference:

C. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill, 1978.