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