The functions ode23 and ode45 are the principal MATLAB and Simulink tools for solving nonstiff ordinary differential equations.

Contents

Single Step Methods

The two functions ode23 and ode45 are single step ODE solvers. They are also known as Runge-Kutta methods. Each step is almost independent of the previous steps. Two important pieces of information are passed from one step to the next. The step size $h$ expected to achieve a desired accuracy is passed from step to step. And, in a strategy known as FSAL, for First Same as Last, the final function value at the end of a successful step is used at the initial function value at the following step.

ODE23

The BS23 algorithm is due to Larry Shampine and his student Przemyslaw Bogacki. Przemyslaw is now a Professor at Old Dominion University. The " 23 " in the function name indicates that two simultaneous single-step formulas, one of second order and one of third order, are involved.

The method has three stages, but there are four slopes $s_i$ because, after the first step, the $s_1$ for one step is the $s_4$ from the previous step. The essentials are

$$ s_1 = f(t_n,y_n) $$

$$ s_2 = f\left(t_n+{h \over 2},y_n+{h \over 2} s_1\right) $$

$$ s_3 = f\left(t_n+{3 \over 4} h, y_n+{3 \over 4} h s_2\right) $$

$$ t_{n+1} = t_n + h $$

$$ y_{n+1} = y_n + {h \over 9} (2 s_1 + 3 s_2 + 4 s_3) $$

$$ s_4 = f(t_{n+1},y_{n+1}) $$

$$ e_{n+1} = {h \over 72} (-5 s_1 + 6 s_2 + 8 s_3 - 9 s_4) $$

Here is a graphical view of the steps.

   ode23steps

These graphics show the starting situation and the three stages. We start at a point $(t_n,y_n)$ with an initial slope $s_1 = f(t_n,y_n)$ and an estimate of a good step size, $h$. Our goal is to compute an approximate solution $y_{n+1}$ at $t_{n+1} = t_n+h$ that agrees with the true solution $y(t_{n+1})$ to within the specified tolerances.

The first stage uses the initial slope $s_1$ to take an Euler step halfway across the interval. The function is evaluated there to get the second slope, $s_2$. This slope is used to take an Euler step three-quarters of the way across the interval. The function is evaluated again to get the third slope, $s_3$. A weighted average of the three slopes is used to step all the way across the interval to get a tentative value for $y_{n+1}$.

$$ s = {1 \over 9} (2 s_1 + 3 s_2 + 4 s_3) $$

The function is evaluated once more to get $s_4$. The error estimate then uses all four slopes.

$$ e_{n+1} = {h \over 72} (-5 s_1 + 6 s_2 + 8 s_3 - 9 s_4) $$

If the error is within the specified tolerance, the step is successful, the tentative value of $y_{n+1}$ is accepted, and $s_4$ becomes the $s_1$ of the next step. If the error is too large, then the tentative $y_{n+1}$ is rejected and the step must be redone. In either case, the error estimate $e_{n+1}$ provides the basis for determining the step size $h$ for the next step.

Notice that $y_{n+1}$ is computed before $s_4$. The final function evaluation is used for the error estimate, but not for the step itself. But this $s_4$ is the $s_1$ in the next step.

Error Estimate

The coefficients in the error estimate $e_{n+1}$ result from the difference between the third order formula that is used to compute $y_{n+1}$ and an independent second order formula that involves the same function values. That second order result is not actually computed because its value is not needed.

Interpolant

All the solvers in the suite provide an interpolant that can generate values approximating the solution to the differential equation to the desired accuracy anywhere in the interval without requiring further evaluation of the function defining the ode. In the case of ode23 this interpolant happens to be "free". In turns out that the Hermite cubic polynomial defined by the four values

$$ y_n, \ F(t_n,y_n), \ y_{n+1}, \ F(t_{n+1}, \ y_{n+1}) $$

does the job. For other solvers in the suite, providing the accompanying interpolant is an important aspect of the algorithm derivation.

Runge-Kutta-Fehlberg

Before today's version of ode45, there was an earlier one. In a 1969 NASA report, Erwin Fehlberg introduced a so-called six stage Runge-Kutta method that requires six function evaluations per step. These function values can be combined with one set of coefficients to produce a fifth-order accurate approximation and with another set of coefficients to produce an independent fourth-order accurate approximation. Comparing these two approximations provides an error estimate and resulting step size control.

Notice that it takes six stages to get fifth order. It is not possible to get fifth order with only five function evaluations per step. The combinatorial complexity of the Taylor series in two variables for $F(t,y)$ overpowers the information available from the function evaluations. In fact, if you continue to investigate the development of Runge-Kutta methods, you will find that, for example, with ten stages it is only possible to achieve seventh order.

In the early 1970's Shampine and his colleague H. A. (Buddy) Watts at Sandia Laboratories published a Fortran code, RKF45, based on Fehlberg's algorithm. In 1977, we made RKF45 the ODE solver in our text book Computer Methods for Mathematical Computations, by Forsythe, Malcolm and Moler. This book is now out of print so I can't provide a decent link to it. But the Fortran source code for RKF45 is still available from netlib.

One thing that I will always remember about RKF45 is ${12 \over 13}$ The fourth function evaluation, $s_4$, is made at the point

$$ t_n + {12 \over 13} h $$

I doubt that you will ever come across ${12 \over 13}$ anyplace else in mathematical software.

RKF45 became the basis for the first version of ODE45 in MATLAB in the early 1980s and for early versions of Simulink. The Felhberg (4,5) pair did a terrific job for almost fifteen years until the late 1990s when Shampine and MathWorker Mark Reichelt modernized the suite and introduced a more efficient algorithm.

ODE45

The new ode45 introduced in the late 1990s is based on an algorithm of Dormand and Prince. It uses six stages, employs the FSAL strategy, provides fourth and fifth order formulas, has local extrapolation and a companion interpolant.

The new ode45 is so effective that it is the only function in the suite where the default value of the refine argument is set to 4. This means that the step size the algorithm naturally wants to choose is so large that the output is more widely spaced than most people prefer. So the interpolant is called up to produce more finely spaced output at no additional cost measured in function evaluations.

The codes for the two routines ode23 and ode45 are very similar. The most important place they differ is the portion that sets the key parameters. Here are the parameters in ode23.

   dbtype 200:209 ode23

200   % Initialize method parameters.
201   pow = 1/3;
202   A = [1/2, 3/4, 1];
203   B = [
204       1/2         0               2/9
205       0           3/4             1/3
206       0           0               4/9
207       0           0               0
208       ]; 
209   E = [-5/72; 1/12; 1/9; -1/8];

The parameter pow is used in the step size calculation. We see that ode23 is a third order method. The array A gives the fractions for each partial step. The array B provides the weights for the slopes. And the array E provides the coefficients in the error estimate. Here are the corresponding parameters for the Dormand-Prince algorithm used in ode45.

   dbtype 201:213 ode45

201   % Initialize method parameters.
202   pow = 1/5;
203   A = [1/5, 3/10, 4/5, 8/9, 1, 1];
204   B = [
205       1/5         3/40    44/45   19372/6561      9017/3168       35/384
206       0           9/40    -56/15  -25360/2187     -355/33         0
207       0           0       32/9    64448/6561      46732/5247      500/1113
208       0           0       0       -212/729        49/176          125/192
209       0           0       0       0               -5103/18656     -2187/6784
210       0           0       0       0               0               11/84
211       0           0       0       0               0               0
212       ];
213   E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];

Comparison

ode23 is a three-stage, third-order, Runge-Kutta method. ode45 is a six-stage, fifth-order, Runge-Kutta method. ode45 does more work per step than ode23, but can take much larger steps. For differential equations with smooth solutions, ode45 is often more accurate than ode23. In fact, it may be so accurate that the interpolant is required to provide the desired resolution. That's a good thing.

ode45 is the anchor of the differential equation suite. The MATLAB documentation recommends ode45 as the first choice. And Simulink blocks set ode45 as the default solver. But I have a fondness for ode23. I like its simplicity. I particularly like it for graphics. The natural step size that ode23 chooses is frequently just right for display purposes. For example, the plot of the Lorenz chaotic attractor at the end of my previous post is done with ode23 choosing the step size.

So give ode23 a try.

Get the MATLAB code

Published with MATLAB® R2014a