# 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

## 댓글

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