MATLAB Examples

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