Ordinary Differential Equation Suite

MATLAB and Simulink have a powerful suite of routines for the numerical solution of ordinary differential equations. Today's post offers an introduction. Subsequent posts will examine several of the routines in more detail.

Contents

Simulation

MATLAB started its life as a "Matrix Laboratory." The capability for the numerical solution of ordinary differential equations was soon added. This, together with the block diagram programming environment, provides the multidomain simulation companion product Simulink. The linear equation solutions from the matrix library underlie the stiff solvers in the ode suite.

Larry Shampine

Larry Shampine is a highly regarded mathematician and computational scientist. He is a Professor, now Emeritus, at Southern Methodist University in Dallas. Before SMU, Larry was with Sandia National Laboratory in Albuquerque. He has been a consultant to the MathWorks for many years. He has written several books and many research and expository papers. His software, originally in Fortran and more recently in MATLAB, has been the basis for our ODE codes, as well as others in the industry. Two of Larry's PhD students, Jacek Kierzenka and Michael Hosea, are on the staff of the MathWorks.

The Suite

There are seven routines in the MATLAB suite. Their names all begin with " ode ". This is followed by two or three digits, and then perhaps one or two letters. The digits indicate the order. Roughly, higher order routines work harder and deliver higher accuracy per step. A suffix "s", "t", or "tb" designates a method for stiff equations. Here is the list of the functions in the suite.

  • ode45 --- The first choice for most nonstiff problems.
  • ode23 --- Less stringent accuracy requirements than ode45.
  • ode113 --- More stringent accuracy requirements than ode45.
  • ode15s --- The first choice for most stiff problems.
  • ode23s --- Less stringent accuracy requirements than ode15s.
  • ode23t --- Moderately stiff problems without numerical damping.
  • ode23tb --- Less stringent accuracy requirements than ode15s.

ode2

Here is a simple code, ode2, which I will use to discuss order and the digits in the suite naming convention. It uses a fixed step size and evaluates the function defining the differential equation twice per step. The first function evaluation at the start of the step provides the slope for a half-step to the midpoint of the interval. The second function evaluation at the midpoint provides the slope for the step all the way across the interval.

   type ode2
   function ode2(F,t0,h,tfinal,y0)
   % ODE2  A simple ODE solver.
   %   ODE2(F,t0,h,tfinal,y0) uses a midpoint rule
   %   with fixed step size h on the interval
   %      t0 <= t <= tfinal
   %   to solve
   %      dy/dt = F(t,y)
   %   with y(t0) = y0.
   
      y = y0;
      for t = t0 : h : tfinal-h
         s1 = F(t,y);
         ymid = y + h*s1/2;
         s2 = F(t+h/2, ymid);
         y = y + h*s2;
         disp(y)
      end

At first glance you might think that the "2" in the name ode2 just comes from the fact that the function F(t,y) is evaluated twice each step. But the reason for the 2 is deeper than that. It is actually because this is a second order solver. Let's see what second order means. Try solving one of the world's easiest differential equations,

$$ dy/dt = 1 $$

With initial condition $y(0) = 0$, the solution is

$$ y(t) = t $$

I'm going to solve this on the interval $0 \le t \le 10$ with a step size $h = 1$, so the result displays as integers.

   format short
   F = @(t,y) 1;
   ode2(F,0,1,10,0)
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10

Well, we got that right. Now try

$$ dy/dt = 2t $$

With initial condition $y(0) = 0$, the solution is

$$ y(t) = t^2 $$

   F = @(t,y) 2*t;
   ode2(F,0,1,10,0)
     1
     4
     9
    16
    25
    36
    49
    64
    81
   100

Again, we got that right. So try

$$ dy/dt = 3t^2 $$

Do we generate

$$ y(t) = t^3 $$

   F = @(t,y) 3*t.^2;
   ode2(F,0,1,10,0)
    0.7500
    7.5000
   26.2500
    63
  123.7500
  214.5000
  341.2500
   510
  726.7500
  997.5000

No, we do not get $t^3$ exactly.

So, our simple function ode2 computes the exact answer for an ODE whose solution is a polynomial in $t$ of degree 2, but not 3. This is what we mean by a solver of second order.

Here's a homework exercise. Write the other simple ode2, the one based on the trapezoid rule. Evaluate F(t,y) at the start of the interval. Then take a step all the way across. Evaluate F(t,y) a second time at the other end of the interval. Then use the average of the two slopes to take the actual step. You should find that the trapezoid version of ode2 has the save behavior as the midpoint version. It solves polynomials of degree 1 and 2 exactly, but not polynomials of degree 3.

In both of these cases, it turns out that "2" is both the number of function evaluations per step and the order. This is not always the case. In later posts, we will see that as we strive for higher accuracy and more efficient methods, the number of function evaluations required increases faster than the order.

ode4, Classical Runge-Kutta

The classical Runge-Kutta method was developed independently by a famous German applied mathematician, Carl Runge, in 1895, and another German applied mathematician, who was not quite as famous, Wilhelm Kutta, in 1901. It was used extensively for hand computations and is still probably in widespread use on digital computers today. You can see from this MATLAB version, a function called ode4, that it evaluates F(t,y) four times per step.

   type ode4
   function ode4(F,t0,h,tfinal,y0)
   % ODE4  Classical Runge-Kutta ODE solver.
   %   ODE4(F,t0,h,tfinal,y0) uses the classsical Runge-Kutta
   %   method with fixed step size h on the interval
   %      t0 <= t <= tfinal
   %   to solve
   %      dy/dt = F(t,y)
   %   with y(t0) = y0.
   
      y = y0;
      for t = t0 : h : tfinal-h
         s1 = F(t,y);
         s2 = F(t+h/2, y+h*s1/2);
         s3 = F(t+h/2, y+h*s2/2);
         s4 = F(t+h, y+h*s3);
         y = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
         disp(y)
      end

You will find that this function integrates $1$, $t$, $t^2$ and $t^3$ exactly. But how about $t^4$?

   format long e
   F = @(t,y) 5*t^4;
   ode4(F,0,1,10,0)
     1.041666666666667e+00
     3.208333333333334e+01
     2.431250000000000e+02
     1.024166666666667e+03
     3.125208333333333e+03
     7.776250000000000e+03
     1.680729166666666e+04
     3.276833333333333e+04
     5.904937500000000e+04
     1.000004166666667e+05

The exact solution would be the integers, $t^5$, but we didn't quite get them. So, classical Runge-Kutta computes the exact answer for an ODE whose solution is a polynomial in $t$ of degree 4, but not 5. This is what we mean by a solver of fourth order, and is the reason why this function is called ode4.

The name of ode4 has only one digit, not two, and classical Runge-Kutta does not qualify for a spot in the MATLAB suite because the method has a fixed step size. There is no error estimate. Modern numerical methods, and modern mathematical software, include error estimates and automatic step size control. This is the subject of my next blogs.

Lorenz graphic

I want to include some graphics in today's blog, so here let's use ode23 to plot the three components of the Lorenz chaotic differential equation described in my previous blog post.

   type lorenzplot

   lorenzplot
   function lorenzplot
   
      % LORENZPLOT  Plot the three components of the solution
      % to the Lorenz equations.
   
      sigma = 10;
      beta = 8/3;
      rho = 28;
      eta = sqrt(beta*(rho-1));
      A = [ -beta    0     eta
               0  -sigma   sigma 
            -eta   rho    -1  ];
         
      v0 = [rho-1 eta eta]';
      y0 = v0 + [3 2 -4]';
      tspan = [0 30];
      
      [t,y] = ode23(@lorenzeqn, tspan, y0);
   
      plot(t,[y(:,1)+15 y(:,2) y(:,3)-40]);
      axis([0 30 -80 80])
      set(gca,'ytick',[-40 0 40],'yticklabel','y3|y2|y1')
      xlabel('t')
      title('Lorenz Attractor') 
   
      % ------------------------------------
   
      function ydot = lorenzeqn(t,y)
      
         A(1,3) = y(2);
         A(3,1) = -y(2);
         ydot = A*y;
    
      end
   
   end




Published with MATLAB® R2014a

|
  • print

コメント

コメントを残すには、ここ をクリックして MathWorks アカウントにサインインするか新しい MathWorks アカウントを作成します。