My Favorite ODE

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 ˙y directly, so it is implicit.

˙y2+y2=1

We don't need a computer to verify that the analytic solution is a combination of sint and cost. 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 ˙y we get two solutions, differing in sign.

˙y=±1y2

That ± 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, sint. 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 sint beyond π/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,

˙y=f(y)

where

f(y)=1y2

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

J=fy

In this case J is 1-by-1.

J=y1y2

When y approaches ±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)=cost 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.

ddt(˙y2+y2=1)

2˙y¨y+2y˙y=0

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

¨y+y=0

First order system.

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

˙y1=y2,  ˙y2=y1

Now there is no trouble computing sint and cost.

   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 π 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




Published with MATLAB® R2016a

|
  • print

댓글

댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.

Loading...
Go to top of page