11.5. Numerical Differential Equations¶

Reading Assignment

Please read chapter 8 of Physical Modeling in MATLAB, by Allen B. Downey [DOWNEY11].

An ordinary differential equation (ODE) is an equation of the form

$y'(t) = \frac{dy}{dt} = f(t, y), \text{ with } y(t_0) = y_0.$

We are given the rate of change (the derivative function $f(t, y)$ ) of $y(t)$ from which we solve for $y(t)$ . Because we are also given the initial value $y(t_0) = y_0$ , this class of ODE problem is called an initial value problem (IVP). We will consider cases of the function $y(t)$ that represents a single equation and also a system of equations.

Differential equations are important to engineers and scientists because changes to the state of dynamic systems are described by differential equations. Examples of these dynamic systems include objects that move; wires carrying fluctuating electrical current; reservoirs of fluid with measurable pressure, temperature, and volume; or anything where the state of the system can be quantified numerically.

Numerical methods are especially needed to solve higher order and nonlinear ODEs. As is the case with integration, our ability to properly evaluate some ODEs using analytic methods is less than satisfying [STRANG14]. So numerical methods for solving ODEs is an important tool in the engineer’s computational toolbox.

This is our third encounter with ODEs. In Throwing a Ball Up, we used MATLAB’s Symbolic Math Toolbox to find the algebraic equations from the differential equations associated with the effect of gravity on a ball thrown upwards. Then in Systems of Linear ODEs we considered systems of first order linear ODEs that may be solved analytically with eigenvalues and eigenvectors. Now, we will consider numerical methods to find and plot the solutions to ODEs. We will cover the basic algorithmic concepts of how numerical ODE solvers work and review the efficient suite of ODE solver provided by MATLAB.

Recall the definition of a derivative:

$\frac{dy}{dt} = \lim_{h \rightarrow 0} \frac{y(t + h) - y(t)}{h}.$

Considered in terms of discrete samples of $y$ , the derivative becomes

$y_k' = \frac{y_{k+1} - y_k}{h},$

where $h$ is the distance between sample points, $h = t_{k+1} - t_k$ .

Rearranged to find $y_{k+1}$ from the previous value, $y_k$ , and the slope of the function at $t_k$ , we have

(11.1)¶ $y_{k+1} = y_k + h\,y'_k$

This equation is the basis for the simplest ODE solver, Euler’s method. However, due to concerns about excessive computation and round-off errors, we are limited on how small $h$ can be. So the accuracy of (11.1) is not acceptable for most applications.

Since differential equations are given in terms of an equation for the derivative of $y$ , $y' = f(t, y)$ , and the initial value, $y_0$ , the general discrete solution comes from the fundamental theorem of calculus, which applies the antiderivative of a function to a definite integral.

(11.2)¶ $y_{k+1} = y_k + \int_{t_k}^{t_{k+1}} f(t, y) dt$

Then the numerical approximation comes from the algorithm used to calculate the definite integral.

(11.3)¶ $I_k = \int_{t_k}^{t_{k+1}} f(t, y) dt$

At each step, the ODE solver uses a straight line with a slope of $I_k/h$ to advance $y_k$ to $y_{k+1}$ .

11.5.1. Euler’s Method¶

Euler’s method uses a Reimann integral to approximate (11.3) with constant height rectangular subintervals. It is a reasonable solution if $f(t, y)$ is constant or nearly constant. However, the results from Euler’s methods applied to most ODEs are not accurate enough for most applications unless the step size is small enough relative to the rate of change of the $f(t, y)$ function.

$\int_{t_k}^{t_{k+1}} f(t, y) dt \approx h\,f(t_k, y_k)$

$y_{k+1} = y_k + h\,f(t_k, y_k)$

Euler’s method is a perfect match if we let $f(t, y) = a$ , then

$y_{k+1} = y_k + h\,a.$

Then by induction,

$\spalignsys{y_1 = {y_0 + h\,a} {}; y_2 = {y_1 + h\,a} {= y_0 + 2\,h\,a}; \. \vdots; y_k = {y_0 + k\,h\,a} {}}.$

The eulerODE function listed below is an implementation of Euler’s method. The first three arguments passed to eulerODE are the same as is used with the MATLAB ODE solvers. An additional argument for the number of values to find is also passed to eulerODE. First, define a function or function handle. The function defines the derivative of the solution. It should take two arguments, $t$ and $y$ . Although in many cases, the function only operates on one variable—usually the second, y, variable. The function should accept and return vector variables.

$\frac{dy}{dt} = f(t, y)$

Next, specify a tspan argument, which is a row vector of the initial and final $t$ values ([t_0 t_final]). Then, the value of $y$ at the initial time is specified $y(t_0)$ . The final argument is $n$ , the number of values to find. The returned data from all of the ODE solvers are the $t$ evaluation points and the $y$ solution.

>> f = @(t, y) 5
>> [t, y] = eulerODE(f, [0 10], 5, 4)
t =
         0    3.3333    6.6667   10.0000
y =
    5.0000   21.6667   38.3333   55.0000

Fig. 11.18 shows a plot comparing the relative accuracy of Euler’s method to other fixed step size algorithms.

function [t, y] = eulerODE(f, tspan, y0, n)
% EULERODE - ODE solver using Euler's method.
%
%   f - vectorized derivative function, y' = f(t, y)
%   tspan - row vector [t0 tf] for initial and final time
%   y0 - initial value of y
%   n - number of y values to find
    h = (tspan(2) - tspan(1))/(n-1);
    t = linspace(tspan(1), tspan(2), n);
    y = zeros(1, n);
    y(1) = y0;
    for k = 1:n-1
        y(k+1) = y(k) + h*f(t(k), y(k));
    end

11.5.2. Heun’s Method¶

Heun’s method uses trapezoid rule integration to approximate (11.3).

(11.4)¶ $y_{k+1} = y_k + \frac{h}{2}\left[f(t_k, y_k) + f(t_{k+1}, y_{k+1})\right]$

However, (11.4) makes use of $y_{k+1}$ , which is not yet know. So Euler’s method is first used to make an initial estimate of $y_{k+1}$ . Because the algorithm makes an initial prediction of $y_{k+1}$ and then updates the estimate with a later equation, Heun’s method belongs to a class of algorithms called predictor–corrector methods.

$y_{k+1} = y_k + h\,f(t_k, y_k)$

The complete Heun method approximation equation is given by (11.5).

(11.5)¶ $y_{k+1} = y_k + \frac{h}{2}\left[f(t_k, y_k) + f(t_{k+1}, (y_k + h\,f(t_k, y_k))\right]$

The heunODE function listed below is called with the same arguments as the eulerODE function.

function [t, y] = heunODE(f, tspan, y0, n)
% HEUNODE - ODE solver using Heun's method.
%
%   f - vectorized derivative function, y' = f(t, y)
%   tspan - row vector [t0 tf] for initial and final time
%   y0 - initial value of y
%   n - number of y values to find
    h = (tspan(2) - tspan(1))/(n-1);
    t = linspace(tspan(1), tspan(2), n);
    y = zeros(1, n);
    y(1) = y0;
    for k = 1:n-1
        fk = f(t(k), y(k));
        yk1 = y(k) + h*fk;
        y(k+1) = y(k) + (h/2)*(fk + f(t(k+1), yk1));
    end

11.5.3. Runge-Kutta Method¶

The Runge-Kutta method uses Simpson’s rule integration to approximate (11.3). There are actually several algorithms that are considered to be part of the Runge-Kutta family of algorithms. The so called classic algorithm used here is an order–4 algorithm with a comparatively simple expression. Runge-Kutta algorithms are weighted sums of coefficients calculated from evaluations of $f(t, y)$ (the derivative function). Later coefficients make use of previous coefficients, so the coefficients at each evaluation point are first calculated sequentially before the solution value can be determined.

Runge-Kutta algorithms are the most accurate of the fixed step size algorithms. They are only surpassed by adaptive step size algorithms, which typically make use of Runge-Kutta equations.

$\begin{array}{rl} s_1 &= f(t_k, y_k) \\ \hfill s_2 &= f(t_k + \frac{h}{2}, y_k + \frac{h}{2}\,s_1) \\ \hfill s_3 &= f(t_k + \frac{h}{2}, y_k + \frac{h}{2}\,s_2) \\ \hfill s_4 &= f(t_k + h, y_k + h\,s_3) \\ \hfill y_{k + 1} &= y_k + \frac{h}{6}\, \left( s_1 + 2s_2 + 2s_3 + s_4 \right) \end{array}$

The rk4ODE function listed below is called with the same arguments as the eulerODE and heunODE functions.

function [t, y] = rk4ODE(f, tspan, y0, n)
% RK4ODE - ODE solver using the classic Runge-Kutta method.
%
%   f - vectorized derivative function, y' = f(t, y)
%   tspan - row vector [t0 tf] for initial and final time
%   y0 - initial value of y
%   n - number of y values to find
    h = (tspan(2) - tspan(1))/(n - 1);
    t = linspace(tspan(1), tspan(2), n);
    y = zeros(1, n);
    y(1) = y0;
    for k = 1:n-1
        s1 = f(t(k), y(k));
        s2 = f(t(k)+h/2, y(k)+s1*h/2);
        s3 = f(t(k)+h/2, y(k)+s2*h/2);
        s4 = f(t(k+1), y(k)+s3*h/2);
        y(k+1) = y(k) + (h/6)*(s1 + 2*s2 + 2*s3 + s4);
    end

11.5.4. Algorithm Comparison¶

We will use the well known decaying exponential ODE to compare the three methods that use a fixed step size. We will also reference this simple ODE when considering stability constraints on the step size in Stability Constraints.

The ODE

$y' = \frac{dy(t)}{dt} = -a\,y,\; \text{where}\; a > 0,\: \text{and}\: y(0) = c$

has the solution $y(t) = c\,e^{-at}$ .

For our comparison, $a = 1$ , and $y(0) = 1$ .

We can see in Fig. 11.18 that, as expected, the classic Runge-Kutta method preforms the best followed by Heun’s method and then Euler’s method.

>> f = @(t, y) -y;
>> [t, y_euler] = eulerODE(f, [0 4], 1, 8);
>> y_exact = exp(-t);
>> [t, y_heun] = heunODE(f, [0 4], 1, 8);
>> [t, y_rk4] = rk4ODE(f, [0 4], 1, 8);

Fig. 11.18 Comparison of fixed interval ODE solvers.

Note

All of the algorithms shown in Fig. 11.18 have very good performance if $h$ is significantly reduced. Finding noticeable differences in the plot is only achieved by using a larger step size.

11.5.5. Stability Constraints¶

Stability and stability constraints are important to us because we do not want our numerical solution to give unstable and wrong answers. A user without knowledge of the stability constraints could specify a step size for the fixed step size algorithm that is too large to get a stable numerical solution. The stability constraint is especially concerning with stiff systems where one term in the solution may place a much higher constraint on the step size than other terms in the solution.

We return to the Euler method and an exponentially decaying ODE to find a stability constraint related to $h$ , the distance between consecutive evaluation points of ODE solvers [SHEN15]. With $y' = -a\,y$ , where $a > 0$ , $y(0) = 1$ , the exact solution is $y(t) = e^{-a\,t}$ . A property of $y(t)$ that needs to be preserved is that $y \to 0$ as $t \to \infty$ .

Referring to Euler’s Method, The forward Euler solution is $y_{k+1} = y_k - h\,a\,y_k$ , or $y_{k+1} = (1 - ha)\,y_k$ .

By induction,

$y_k = (1 - ha)^k\,y_0$

A stability constraint is needed to ensure that the property of $y \to 0$ as $k \to \infty$ is preserved.

$\begin{array}{rl} &\abs{1 - h\,a} < 1 \\ -1 &< (1 - h\,a) < 1 \\ -2 &< -h\,a < 0 \\ &h < \frac{2}{a} \end{array}$

In the following example $2/a = 2/0.5 = 4$ , but $h = 5$ and the results are unstable.

>> fexp = @(t, y) -0.5*y;
>> [t, y_unstable] = eulerODE(fexp, [0 20], 1, 5)
 t =
      0     5    10    15    20
 y_unstable =
      1.0000   -1.5000    2.2500   -3.3750    5.0625

The other ODE solvers evaluate the function at points between the set evaluation points, so the stability constraint for Euler’s method is a sufficient constraint for them.

11.5.6. Implicit Backward Euler Method¶

Implicit ODE solvers look at (11.2) from the perspective of finding $y_k$ from the yet unknown $y_{k+1}$ . So the terms of (11.2) are reordered and the limits of the integral are reversed. The resulting algorithm removes the stability constraint on the step size.

(11.6)¶ $y_k = y_{k+1} + \int_{t_{k+1}}^{t_k} f(t, y) dt$

Again, we will use the ODE $y' = -a\,y$ , where $a > 0$ , and $y(0) = 1$ . The exact solution is $y(t) = e^{-a\,t}$ . As before, we want to find any stability constraints that are required to ensure that $y \to 0$ as $k \to \infty$ . We can solve the decaying exponential ODE problem numerically by applying Euler’s method to (11.6).

(11.7)¶ $\begin{array}{rl} y_k &= y_{k+1} + h\,a\,y_{k+1} \\ y_k &= (1 + h\,a)\,y_{k+1} \\ y_{k+1} &= \frac{1}{1 + h\,a}\,y_k \\ y_{k+1} &= \left(\frac{1}{1 + h\,a}\right)^k y_0 \\ \end{array}$

To satisfy the stability constraint, we require that

$\left\vert\frac{1}{1 + h\,a}\right\vert < 1,$

which always holds because both $h$ and $a$ are positive numbers. So the implicit backward Euler method is unconditionally stable.

11.5.7. Adaptive Algorithms¶

Numerical ODE solvers need to use appropriate step sizes between evaluation points to achieve high levels of accuracy while also quickly finding solutions across the evaluation interval. Adaptive ODE solvers reduce the step size, $h$ , where needed to meet the accuracy requirements. To improve performance, they also increase the step size when the accuracy requirements can be satisfied with a larger step size.

One of two strategies can be used to determine the step size adjustments. As was done with adaptive numerical methods for computing integrals in Recursive Adaptive Integral, the solution from two half steps could be compared to the solution from a full step. Based on the comparison, the step size can be reduced, left the same, or increased. However, points of an ODE solution are found sequentially, so the strategy of splitting a wide region into successively smaller regions would be computationally slow. Instead, the current step size is used to find the next solution ( $y_{k+1}$ ) using two algorithms where one algorithm is regarded as more accurate than the other. The guiding principle is that if the solutions computed by the two algorithms are close to the same, then the step size is small enough or could even be increased. But if the difference between the solutions is larger than a threshold, then the step size should be reduced. If the solutions differ significantly, then the current evaluations are recalculated with a smaller step size. The algorithm advances to the next evaluation point only when the differences between the solutions is less than a specified tolerance.

11.5.7.1. The RKF45 Method¶

The Runge-Kutta-Fehlberg method (denoted as RKF45) uses fourth and fifth order Runge-Kutta ODE solvers to adjust the step size at each evaluation point. The fourth order equation is more complex than is used by the fourth order classic Runge-Kutta algorithm. But the computed coefficients for the fourth order equation are also used for the fifth order equation, which only requires one additional coefficient. The $s_2$ coefficient is used in finding later coefficients, but is not part of either the fourth or fifth order solution equations. In the rkf45ODE function, the fourth order solution values (RK4) are stored in the z array, and the fifth order solution values (RKF45) are stored in the y array, which is returned from the function. The coefficients and solution values are calculated using the following equations.

$\begin{array}{rl} &s_1 = h\,f\left(t_k, y_k\right) \\ \hfill &s_2 = h\,f\left(t_k + \frac{1}{4}h, y_k + \frac{1}{4}s_1\right) \\ \hfill &s_3 = h\,f\left(t_k + \frac{3}{8}h, y_k + \frac{3}{32}s_1 + \frac{9}{32}s_2\right) \\ \hfill &s_4 = h\,f\left(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}s_1 - \frac{7200}{2197}s_2 + \frac{7296}{2197}s_3\right) \\ \hfill &s_5 = h\,f\left(t_k + h, y_k + \frac{439}{216}s_1 - 8s_2 + \frac{3680}{513}s_3 - \frac{845}{4104}s_4\right) \\ \hfill &s_6 = h\,f\left(t_k + \frac{1}{2}h, y_k - \frac{8}{27}s_1 + 2s_2 - \frac{3544}{2565}s_3 + \frac{1859}{4104}s_4 - \frac{11}{40}s_5\right) \\ \hfill &z_{k+1} = y_k + \frac{25}{216}s_1 + + \frac{1408}{2565}s_3 + \frac{2197}{4101}s_4 - \frac{1}{5}s_5 \\ \hfill &y_{k+1} = y_k + \frac{16}{135}s_1 + \frac{6656}{12825}s_3 + \frac{28561}{56430}s_4 - \frac{9}{50}s_5 + \frac{2}{55}s_6 \end{array}$

There are variances in the literature on the implementation details of the RKF45 method. The rkf45ODE function listed below is adapted from the description of the algorithm in the numerical methods text by Mathews and Fink [MATHEWS05].

Some adaptive ODE algorithms make course changes to the step size such as either cutting it in half or doubling it. The rkf45ODE function first finds a scalar value, $s$ , from a tolerance value and the difference between the solutions from the two algorithms. The new step size comes from multiplying the current step size by the scalar. The scalar will be one if the difference between the solutions is half of the tolerance. The algorithm advances to the next step if the difference between the two solutions is less than the tolerance value, which occurs when the scalar is greater than 0.84. Otherwise, the current evaluation is recalculated with a smaller step size. The multiplier is restricted to not be more than two. So the step size will never increase to more than twice its current size. The tolerance and the threshold for advancing to the next point are tunable parameters.

$s = \left(\frac{tolerance}{2\,\norm{z_{k+1}- y_{k+1}}} \right)^{\frac{1}{4}}$

The initial value ( $y_0$ ) passed to the rkf45ODE function can be either a scalar for a single equation ODE or a column vector for a system of ODEs. Since the number of steps to be taken is not known in advance, the initial step size is a function argument instead of the number of solution values to find.

function [t, y] = rkf45ODE(f, tspan, y0, h0)
% rk5ODE - Adaptive interval ODE solver using the
%          Runge-Kutta-Fehlberg method.
%
%   f - vectorized derivative function, y' = f(t, y)
%   tspan - row vector [t0 tf] for initial and final time
%   y0 - initial value of y, scalar or column vector
%   h0 - intial step size
    tol = 2e-5;
    h = h0;
    nEst = ceil((tspan(2) - tspan(1))/h);
    m = size(y0, 1);
    t = zeros(1, 2*nEst);  % estimated preallocation
    y = zeros(m, 2*nEst);
    z = zeros(m, 2*nEst);
    t(1) = tspan(1);
    y(:,1) = y0;      % RKF45 results
    z(:,1) = y0;      % RK4 results
    k = 1;
    while t(k) < tspan(2)
        s1 = h*f(t(k), y(:,k));
        s2 = h*f(t(k) + h/4, y(:,k) + s1/4);
        s3 = h*f(t(k) + h*3/8, y(:,k) + s1*3/32 + s2*9/32);
        s4 = h*f(t(k) + h*12/13, y(:,k) + s1*1932/2197 ...
            - s2*7200/2197 + s3*7296/2197);
        s5 = h*f(t(k) + h, y(:,k) + s1*439/216 - s2*8 ...
            + s3*3680/513 - s4*845/4104);
        s6 = h*f(t(k) + h/2, y(:,k) - s1*8/27 + s2*2 ...
            - s3*3544/2565 + s4*1859/4104 - s5*11/40);
        z(:,k+1) = y(:,k) + s1*25/216 + s3*1408/2565 ...
            + s4*2197/4101 - s5/5;  % RK4
        y(:,k+1) = y(:,k) + s1*16/135 + s3*6656/12825 ...
            + s4*28561/56430 - s5*9/50 + s6*2/55; % RKF5
        % no divide by 0, h at most doubles
        algDiff = max(tol/32, norm(z(:,k+1) - y(:,k+1)));
        s = (tol/(2*algDiff))^0.25;
        h = h*s;
        if s > 0.84
            t(k+1) = t(k) + h;
            k = k + 1;
        end
    end
    % delete unused preallocated memory
    t(k+1:end) = []; y(:,k+1:end) = [];

Fig. 11.19 shows that the adaptive Runge-Kutta-Fehlberg (RKF45) method is more accurate than the classic Runge-Kutta method (RK4) with the initial step size of RKF45 being the same as the fixed step size used by RK4. For most ODEs, the ability of RKF45 to reduce the step size where needed likely contributes more to its improved accuracy than can be attributed to using a higher order equation.

>> f = @(t,y)-y;
>> [t, y_rkf45] = rkf45ODE(f, [0 4], 1, 0.5);
>> y_exact = exp(-t);
>> hold on
>> plot(t, y_exact)
>> plot(t, y_rkf45, ':')
>> [t, y_rk4] = rk4ODE(f, [0 4], 1, 8);
>> plot(t, y_rk4, '--')
>> hold off, axis tight, xlabel('t')
>> ylabel('y(t)'), legend('exact', 'RKF45', 'RK4')

Fig. 11.19 The adjustable step size of the RKF45 algorithm gives accurate results with good performance.

11.5.7.2. MATLAB’s `ode45`¶

Although the adaptive rkf45ODE function has good results, it is not as accurate as the MATLAB ode45 function. Using the same comparison test shown in Fig. 11.19, the ode45 function had a maximum error of $6{\times}10^{-6}$ , and a plot of its solution appears to be directly on top of the exact solution plot. MATLAB’s ode45 is the ODE solver that should be tried first for most problems. It has good accuracy and performance for non-stiff ODEs. Like the RKF45 method, it uses an adaptive step size algorithm based on the fourth and fifth order Runge-Kutta formula.

The ode45 function is called with the same arguments as the ODE solvers previously described, except neither the initial step size nor the number of evaluation points is given. The MATLAB ODE solvers also offer additional flexibility regarding the evaluation points. If the tspan argument contains more than two values, then tspan is taken to be a vector listing the evaluation points to use.

As before, we often save two return values from ode45, [t, y] = ode45(f, tspan, y0). The return values in are column vectors rather than row vectors as was the case with the simple functions previously used. If no return values are specified, ode45 will make a plot of the $y$ values against the $t$ –axis.

The MATLAB ODE solvers have several additional options and tunable parameters. Rather than passing those individually to the solver, a data structure holding the options is passed as an optional argument. The data structure is created, and modified using the odeset command with name–value pairs. Refer to the MATLAB documentation for odeset for a full list of the options and their descriptions.

>> options = odeset('RelTol', 1e-5);          % create
>> options = odeset(options, 'AbsTol',1e-7);  % modify
>> [t, y] = ode45(f, tspan, y0, options);     % use

11.5.8. The Logistic Equation¶

This example problem is a model of growth slowed down by competition. The competition is within the population for limited resources. The growth gives a $y$ term in the ODE—the growth is proportional to the population size. The competition ( $y$ against $y$ ) adds a $-y^2$ term.

$\frac{dy}{dt} = a\,y - b\,y^2.$

The solution is modeled with positive and negative values for $t$ . The analytic solution is: ([STRANG14] page 55)

$y(t) = \frac{a}{b + d\,e^{-t}}.$

Where, $d = \frac{a}{y(0)} - b$ .

In the distant past, ( $t = -\infty$ ), the population size is modeled as zero. The steady state ( $t = \infty$ ) population size is $y(t) = K = \frac{a}{b}$ , which is the sustainable population size.

The half way point, can be modeled as $y(0) = a/2b$ .

We will compute the positive half of the solution and then reflect the values to get the negative time values.

% File: logistic_diffEq.m
% Logistic Equation: a = 1, b = 1
f = @(t, y) y - y.^2;
tspan = [0 5];
y0 = 1/2;
[t, y] = ode45(f, tspan, y0);
figure, plot(t,y)
t1 = flip(-t);
y1 = flip(1-y);
t2 = [t1;t];
y2 = [y1;y];
analytic = 1./(1 + exp(-t2));
figure, hold on
plot(t2, y2, 'o')
plot(t2, analytic)
hold off, title('Logistic Equation')
xlabel('t'), ylabel('y')
legend('ode45 Numeric Solution','Analytic Solution', ...
    'Location', 'northwest');

Fig. 11.20 Numerical and analytical solution to the logistic ODE.

11.5.9. Stiff ODEs¶

A problem is stiff if the solution contains multiple terms with different step size requirements. A slower moving component may dominate the solution; however, a less significant component of the solution may change much faster. Adaptive step size algorithms must take small steps to maintain accurate results. Because stiff ODEs problems have solutions with multiple terms, they come from either second or higher order ODEs or from systems of ODEs.

We will now look at an example of a moderately stiff system of ODE equations. Then we will see how implicit solvers solve stiff systems of ODE and are unconditionally stable. Then we will loot at MATLAB’s stiff ODE solvers, such as ode15s.

11.5.9.1. Stiff ODE Systems¶

Here we show a system of ODEs that have a moderately stiff solution. The reader is referred to Systems of Linear ODEs for a review of systems of first order, linear ODEs.

$\spalignsys{x' = -25\,x + 24\,y{,} {\qquad} x(0) = 2; y' = 24\,x - 25\,y{,} {\qquad} y(0) = 1}$

It is more convenient to express systems of ODEs using vectors and matrices.

$\bm{u}(t) = \vector{x(t); y(t)} \qquad \mathbf{A} = \mat{-25 24; 24 -25} \qquad \bm{u'} = \mathbf{A}\,\bm{u} \qquad \bm{u}_0 = \vector{1; 2}$

The eigenvectors, eigenvalues, and $c$ coefficients are found with some help from MATLAB.

$\bm{x}_1 = \frac{1}{\sqrt{2}}\vector{1; -1} \quad \bm{x}_2 = \frac{1}{\sqrt{2}}\vector{1; 1} \quad \spaligntabular{l}{{$\lambda_1 = -49$}; {$\lambda_2 = -1$}} \quad \spaligntabular{l}{{$c_1 = -\sqrt{2}\cdot 0.5$}; {$c_2 = \sqrt{2}\cdot 1.5$}}$

The solution to the system of ODEs then comes from the general solution for systems of first order, linear ODEs as given in Systems of Linear ODEs.

$\bm{u}(t) = \spalignsys{x(t) = -0.5\,e^{-49\,t} + 1.5\,e^{-t}; y(t) = 0.5\,e^{-49\,t} + 1.5\,e^{-t}}$

This system is stiff because the exact solution has terms with largely different exponent values coming from the eigenvalues of $\bf{A}$ . The stability constraints on the $e^{-49\,t}$ terms require $h < 2/49 \approx 0.0408$ . From the beginning of the evaluation, those terms will move towards zero much faster than the $e^{-t}$ terms. We can see from Fig. 11.21 and Fig. 11.22 that an adaptive algorithm such as RKF45 or MATLAB’s ode45 will significantly reduce the step size from the beginning of the evaluations until the $e^{-49\,t}$ term has sufficiently gone to zero. Using a longer evaluation time, the errors between the RKF45 solution and the exact solution go to near zero and the step sizes get much larger.

>> A = [-25 24; 24 -25];
>> fA = @(t, u) A*u;
>> u0 = [1;2];
>> [t, u] = rkf45ODE(fA, [0 1], u0, 0.04);
>> x = @(t) -0.5*exp(-49*t) + 1.5*exp(-t);
>> y = @(t) 0.5*exp(-49*t) + 1.5*exp(-t);

>> figure, hold on
>> plot(t, u(1,:), 'o')
>> plot(t, u(2,:), '+')
>> plot(t, x(t), ':')
>> plot(t, y(t), '--')
>> hold off
>> axis tight
>> legend('u(1,:)', 'u(2,:)', 'x(t)', 'y(t)')
>> xlabel('t')

Fig. 11.21 Stiff ODE system—the $\pm 0.5\,e^{-49\,t}$ terms quickly goes to zero and the $x(t)$ and $y(t)$ solutions approximately merge to $1.5\,e^{-t}$ .

>> U = [x(t); y(t)];
>> figure;
>> plot(t, vecnorm(U - u))
>> hold on
>> plot(t(1:end-1), diff(t), 'o:')
>> axis tight
>> xlabel('t')
>> legend('error vector', 'h', 'Location', 'southeast')
>> hold off

Fig. 11.22 Error and step size for stiff ODE system. First curve—the norm of the error vector between the exact solution and that found from the RKF45 solver. Second curve—the step size, $h$ , for the stiff ODE system.

11.5.9.2. Implicit Solvers for Systems of ODEs¶

The implicit backward Euler method is used here to solve stiff systems of ODEs. The initial discussion parallels that of Implicit Backward Euler Method but the equations are extended to systems of ODEs. As with the implicit solver for a single ODE equation, the implicit backward Euler (IBE) method for systems of ODEs is also unconditionally stable.

We will again use a decaying exponential ODE such as described in Stiff ODE Systems. As was the case in Stiff ODE Systems, the $\bf{A}$ matrix is required be a symmetric matrix with real, negative eigenvalues.

$\bm{u'}(t) = \mathbf{A}\,\bm{u}(t) \qquad \bm{u}(0) = \bm{u}_0$

The solution is $\bm{u}(t) = \bm{u}_0 e^{\mathbf{A}\,t}$ .

Notice that we do not have a negative sign in front of the $\bf{A}$ matrix as we previously saw with the $a$ scalar ODE in Stability Constraints. But it is a decaying exponential that goes to zero when $t$ goes to infinity because the eigenvalues of $\bf{A}$ are negative.

We get the starting equation for the numerical solution by applying (11.6).

$\bm{u}_k = \bm{u}_{k+1} - h\,\mathbf{A}\,\bm{u}_{k+1} = (\mathbf{I} - h\,\mathbf{A})\,\bm{u}_{k+1}$

(11.8)¶ $\bm{u}_{k+1} = (\mathbf{I} - h\,\mathbf{A})^{-1}\,\bm{u}_k$

(11.8) is an implicit equation for solving systems of ODEs numerically, and is implemented in the IBEsysODE function listed below. Notice that the implicit approach requires solving an equation at every step. In this case, the equation is a linear system with a fixed step size, so the IBEsysODE function uses LU decomposition once and calls on MATLAB’s triangular solver with the left-divide operator to quickly find a solution in each loop iteration. If this were a non-linear problem, then numerical methods such as Newton’s method or the secant method described in Numerical Roots of a Function would be needed, which could add significant computation. For this reason, implicit methods are not recommended for non-stiff or non-linear problems [SHEN15].

function [t, u] = IBEsysODE(A, tspan, u0, n)
% IBESYSODE - ODE solver using the Implicit Euler Backward
%             method for systems of ODEs.
%
%   A - matrix defining the system coefficients
%   tspan - row vector [t0 tf] for initial and final time
%   y0 - initial value of y
%   n - number of y values to find
    h = (tspan(2) - tspan(1))/n;
    t = linspace(tspan(1), tspan(2), n);
    m = size(A, 1);
    u = zeros(m, n);
    u(:,1) = u0;
    [L, U, P] = lu(eye(m) - h*A);
    for k = 1:n-1
        u(:,k+1) = U\(L\(P*u(:,k)));
    end

Fig. 11.23 shows the numerical solution that the IBEsysODE function found to our stiff system of ODEs. Since this is a fixed step size solution, the step size was decreased to 0.005 to accurately track the quickly changing $e^{-49\,t}$ terms in the solution. Of course, the small step is not needed once the $e^{-49\,t}$ term was sufficiently close to zero. Remember that there is a difference between accuracy and stability. A lack of stability causes the solution to be wrong as $k \to \infty$ . A lack of accuracy just means that the exact solution is not well modeled by the numerical solution, which is often solved by reducing the step size.

>> A = [-25 24; 24 -25];
>> u0 = [1;2];
>> [t, u] = IBEsysODE(A, [0, 1], u0, 50);
>> x = @(t) -0.5*exp(-49*t) + 1.5*exp(-t);
>> y = @(t) 0.5*exp(-49*t) + 1.5*exp(-t);
>> figure, hold on
>> plot(t, u(1,:), 'o')
>> plot(t, u(2,:), '+')
>> plot(t, x(t), ':')
>> plot(t, y(t), '--')
>> hold off
>> axis tight
>> legend('u(1,:)', 'u(2,:)', 'x(t)', 'y(t)')
>> xlabel('t')

Fig. 11.23 With enough data points, the IBEsysODE function is able to to reasonably well track the rapidly changing solution of the $e^{-49\,t}$ terms.

To evaluate the stability of the implicit backward Euler method, we start by using induction to express the steady state solution in terms of the initial value.

$\bm{u}_k = \left[(\mathbf{I} - h\,\mathbf{A})^{-1}\right]^k\,u_0$

The stability constraint is that $\bm{u}_k \to 0$ as $k \to \infty$ . So we need that

$\norm{(I - h\,A)^{-1}}_2 < 1, \; \text{ for all }h.$

Since our $\bf{A}$ matrix is symmetric, we can use the property of $l_2$ matrix norms that $\norm{\mathbf{A}}_2 = \max_i \abs{\lambda_i(\mathbf{A})}$ , where $\lambda_i(\mathbf{A})$ means the $i^{th}$ eigenvalue of $\bf{A}$ . (See Matrix Norms for more information of matrix norms.)

$\max_i \abs{\lambda_i\left((I - h\,A)^{-1}\right)} < 1$

Now we can make use of the eigenvalue property for symmetric matrices that $\lambda_i(\mathbf{A}^{-1}) = \frac{1}{\lambda_i(\mathbf{A})}$ . (See Properties of Eigenvalues and Eigenvectors)

(11.9)¶ $\begin{array}{rl} & \max_i \abs{\frac{1}{\lambda_i(I - h\,\mathbf{A})}} < 1 \\ \hfill & \min_i \abs{\lambda_i(I - h\,\mathbf{A})} = \min_i \abs{1 - h\,\lambda_i(\mathbf{A})} > 1, \end{array}$

(11.9) holds for any symmetric $\bf{A}$ matrix with negative eigenvalues. So the implicit backward Euler method is unconditionally stable for systems of ODEs.

Unconditionally stable ODE solving algorithms remove the concern about having an unstable numerical solution because the step size is too large, but as we saw in Fig. 11.23, accuracy requirements may still require a reduced step size. The more complex and robust algorithms of MATLAB’s ODE solvers for stiff systems are able to adjust the step size enough to find accurate solutions without being excessively slowed by a stiff systems.

11.5.9.3. MATLAB’s `ode15s`¶

MATLAB has several stiff ODE solvers. The ode15s should be tried first for stiff ODEs. The primary authors of MATLAB’s suite of ODE solvers, Shampine and Reichelt, describe ode15s as an implicit quasi-constant step size algorithm. They further give the following suggestion regarding picking which of their ODE solver to use [SHAMPINE97].

$\ldots$ Except in special circumstances, ode45 should be the code tried first. If there is reason to believe the problem to be stiff, or if the problem turns out to be unexpectedly difficult for ode45, the ode15s code should be tried.

Fig. 11.24 shows how ode15s handles the stiff system of ODEs used in the previous example.

>> A = [-25 24; 24 -25];
>> fA = @(t, u) A*u;
>> u0 = [1;2];
>> [t, y_ode15s] = ode15s(fA, [0 1], u0);
>> t = t'; u = y_ode15s';
>> x = @(t) -0.5*exp(-49*t) + 1.5*exp(-t);
>> y = @(t) 0.5*exp(-49*t) + 1.5*exp(-t);
>> figure, hold on
>> plot(t, u(1,:)', 'o')
>> plot(t, u(2,:)', '+')
>> plot(t, x(t), ':')
>> plot(t, y(t), '--')
>> hold off
>> legend('u(1,:)', 'u(2,:)', 'x(t)', 'y(t)')
>> xlabel('t'), ylabel('y')

Fig. 11.24 The ode15s ODE solver is designed for stiff systems, so the less significant, but rapidly changing terms do not significantly reduce the step size.

11.5.10. MATLAB’s Suite of ODE Solvers¶

MATLAB’s ode45 and ode15s are the nonstiff and stiff ODE solvers to use for most problems. If those do not give satisfactory results then ode23 and ode113 might be tried for nonstiff problems. For stiff problems, ode23s should be the next function tried after ode15s. Additional documentation and research may be needed to select the best solver to use for some problems. The MathWorks’ web pages have quite a bit of documentation on the suit of ODE solvers. The MATLAB Guide by Higham and Higham [HIGHAM17] has extended coverage MATLAB’s suite of ODE solvers. Table 11.1 lists MATLAB’s suite of ODE solvers.

Table 11.1 The MATLAB ODE solvers¶
Solver	Problem type	Description
`ode45`	Nonstiff	Explicit Runge-Kutta, orders 4 and 5
`ode23`	Nonstiff	Explicit Runge-Kutta, orders 2 and 3
`ode113`	Nonstiff	Explicit linear multistep, orders 1–13
`ode15s`	Stiff	Implicit linear multistep, orders 1–5
`ode15i`	Fully implicit	Implicit linear multistep, orders 1–5
`ode23s`	Stiff	Modified Rosenbrock pair, orders 2 and 3
`ode23t`	Mildly stiff	Implicit trapezoid rule, orders 2 and 3
`ode23tb`	Stiff	Implicit Runge-Kutta, orders 2 and 3

Online Resources

This MathWorks web page has some good tutorial videos about numerical solutions to differential equations. Please watch at least the first three videos (Euler – ODE1, Midpoint Method – ODE2, and Classical Runge-Kutta – ODE4). Although, the ODE solvers in MATLAB are slightly more complex than the examples in the first three videos, they explain the basic concepts.

MATLAB has several ODE solvers, they are reviewed on this web page.

11.5. Numerical Differential Equations¶

11.5.1. Euler’s Method¶

11.5.2. Heun’s Method¶

11.5.3. Runge-Kutta Method¶

11.5.4. Algorithm Comparison¶

11.5.5. Stability Constraints¶

11.5.6. Implicit Backward Euler Method¶

11.5.7. Adaptive Algorithms¶

11.5.7.1. The RKF45 Method¶

11.5.7.2. MATLAB’s `ode45`¶

11.5.8. The Logistic Equation¶

11.5.9. Stiff ODEs¶

11.5.9.1. Stiff ODE Systems¶

11.5.9.2. Implicit Solvers for Systems of ODEs¶

11.5.9.3. MATLAB’s `ode15s`¶

11.5.10. MATLAB’s Suite of ODE Solvers¶

Table of Contents

Previous topic

Next topic

This Page

11.5. Numerical Differential Equations¶

11.5.1. Euler’s Method¶

11.5.2. Heun’s Method¶

11.5.3. Runge-Kutta Method¶

11.5.4. Algorithm Comparison¶

11.5.5. Stability Constraints¶

11.5.6. Implicit Backward Euler Method¶

11.5.7. Adaptive Algorithms¶

11.5.7.1. The RKF45 Method¶

11.5.7.2. MATLAB’s ode45¶

11.5.8. The Logistic Equation¶

11.5.9. Stiff ODEs¶

11.5.9.1. Stiff ODE Systems¶

11.5.9.2. Implicit Solvers for Systems of ODEs¶

11.5.9.3. MATLAB’s ode15s¶

11.5.10. MATLAB’s Suite of ODE Solvers¶

11.5.7.2. MATLAB’s `ode45`¶

11.5.9.3. MATLAB’s `ode15s`¶