# My Favorite ODE

Posted by **Cleve Moler**,

My favorite ordinary differential equation provides a good test of ODE solvers, both numeric and symbolic. It also provides a nice illustration of the underlying existence theory and error analysis.

### Contents

#### Implicit equation

Here is one of my favorite ordinary differential equations. The dot means differentiation with respect to $t$. The equation does not express $\dot{y}$ directly, so it is implicit.

$$ \dot{y}^2 + y^2 = 1 $$

We don't need a computer to verify that the analytic solution is a combination of $\sin{t}$ and $\cos{t}$. But the constants $y=1$ and $y=-1$ are also solutions.

#### Explicit equation

Our ODE solvers require an explicit equation, but when we solve this equation for $\dot{y}$ we get two solutions, differing in sign.

$$ \dot{y} = \pm \sqrt{1-y^2} $$

That $\pm$ sign makes this interesting.

#### First attempt

Let's plunge ahead with the $+$ sign. Specify the equation, an interval for $t$ and an initial condition.

```
f = @(t,y) sqrt(1-y^2)
format short
tspan = [0 pi]
y0 = 0
```

f = @(t,y)sqrt(1-y^2) tspan = 0 3.1416 y0 = 0

An attempt to compute a numeric solution with `ode45` runs into trouble when $y$ gets slightly larger than $1$ and the square root of a negative number generates an imaginary result. So temporarily use a shorter $t$ interval.

try ode45(f,tspan,y0); catch ode45(f,[0 1.47],y0) end pi_axis([-0.1 1.1])

Error using odeplot (line 63) Error updating the ODEPLOT window. Solution data may have been corrupted. Argument Y cannot be complex. Error in ode45 (line 435)

#### Keep in bounds

Change the equation so that it provides a guard against the square root of negative numbers.

f = @(t,y) sqrt(1-min(y,1)^2) ode45(f,tspan,y0); pi_axis([-0.1 1.1])

f = @(t,y)sqrt(1-min(y,1)^2)

This is the correct solution, $\sin{t}$. For a while. But when $y$ reaches $1$ it has to stay there because its derivative is always nonnegative.

#### Switch sign

If we want to continue with $\sin{t}$ beyond $\pi/2$ we have to switch to the minus sign.

f = @(t,y) (1-2*double(t>pi/2))*real(sqrt(1-y.^2)) ode45(f,tspan,y0) pi_axis([-0.1 1.1])

f = @(t,y)(1-2*double(t>pi/2))*real(sqrt(1-y.^2))

#### Jacobian

We're solving the autonomous differential equation,

$$ \dot{y} = f(y) $$

where

$$ f(y) = \sqrt{1-y^2} $$

An important quantity in both the existence theory for ODEs and the error analysis for numeric methods is the *Jacobian*.

$$ J = \frac{\partial{f}}{\partial{y}} $$

In this case $J$ is 1-by-1.

$$ J = -\frac{y}{\sqrt{1 - y^2}} $$

When $y$ approaches $\pm 1$ the Jacobian approaches infinity.

Existence and uniqueness theorems for ODEs rely on a bounded Jacobian (or a bounded Lipschitz constant). The derivations of numerical methods rely on a bounded Jacobian. Error control in modern ODE software relies on a bounded Jacobian. All of these are up in the air with this example.

#### Nonuniqueness

We can see the lack of uniqueness happening when we change the initial condition to $y(0) = -1$. The Jacobian is infinite at the start and both the constant $y(t) = -1$ and the trig function $y(t) = -\cos{t}$ are solutions. A small perturbation in the initial condition switches the numerical solution from one to the other.

f = @(t,y) sqrt(1-max(min(y,1),-1)^2) y0 = -1 ode45(f,tspan,y0); pi_axis([-1.2 0]) hold on format long y0 = -.99999999 ode45(f,tspan,y0); pi_axis([-1.1 1.1]) hold off

f = @(t,y)sqrt(1-max(min(y,1),-1)^2) y0 = -1 y0 = -0.999999990000000

#### Symbolic solutions

I was curious to see how the Symbolic Toolbox handles this equation with this initial condition. There is good news and surprising news.

```
syms y(t)
ode = diff(y)^2+y^2==1
init = y(0)==-1
ysym = dsolve(ode,init)
```

ode(t) = diff(y(t), t)^2 + y(t)^2 == 1 init = y(0) == -1 ysym = -1 -cosh(t*1i)

I am pleased to see the two solutions confirmed, but I am surprised to see the second one expressed as a hyperbolic trig function of an imaginary argument.

#### Simplify

Simplification generates a more familiar form.

simplify(ysym)

ans = -1 -cos(t)

#### Second order equation

Let's return to the original implicit equation and differentiate the entire equation with respect to $t$.

$$ \frac{d}{dt} \left( \dot{y}^2 + y^2 = 1 \right) $$

$$ 2 \dot{y} \ddot{y} + 2 y \dot{y} = 0 $$

We can divide through this equation by $2 \dot{y}$, provided that ${\dot{y}}$ is nonzero. This eliminates the constant solutions and produces an old friend, the harmonic oscillator.

$$ \ddot{y} + y = 0 $$

#### First order system.

Our ODE solvers need the second order equation to be expressed as a first order system.

$$ \dot{y_1} = y_2, \ \ \dot{y_2} = -y_1 $$

Now there is no trouble computing $\sin{t}$ and $-\cos{t}$.

f = @(t,y)[y(2); -y(1)] y0 = [0; -1] ode23(f,tspan,y0) pi_axis([-1 1])

f = @(t,y)[y(2);-y(1)] y0 = 0 -1

#### pi_axis

I have been using this little function to label the x-axis with multiples of $\pi$ and scale the y-axis.

```
type pi_axis
```

function pi_axis(yaxis) axis([0 pi yaxis]) set(gca,'xtick',0:pi/4:pi, ... 'xticklabels',{0 '\pi/4' '\pi/2' '3\pi/4' '\pi'}) end

Get the MATLAB code

Published with MATLAB® R2016a

### Note

Comments are closed.

## Recent Comments