FIGURE P25.20 A spherical tank.
754 RUNGE-KUTTA METHODS
from its equilibrium position (m), and g 5 gravitational acceleration
(9.81 m/s2). Solve these equations for the positions and velocities of
the three jumpers given the initial conditions that all positions and
velocities are zero at t 5 0. Use the following parameters for your
calculations: m1 5 60 kg, m2 5 70 kg, m3 5 80 kg, k1 5 k3 5 50,
and k2 5 100 (N/m).
25.24 Given the initial conditions, y(0) 5 1 and y9(0) 5 0, solve
the following initial-value problem from t 5 0 to 4:
dt 2 1 4y 5 0
Obtain your solutions with (a) Euler’s method and (b) the fourth-
order RK method. In both cases, use a step size of 0.125. Plot both
solutions on the same graph along with the exact solution y 5 cos 2t.
25.25 Use the following differential equations to compute the
velocity and position of a soccer ball that is kicked straight up in the
air with an initial velocity of 40 m/s:
dt 5 y
dt 5 2g 2
m y Zy Z
where y 5 upward distance (m), t 5 time (s), y 5 upward velocity
(m/s), g 5 gravitational constant (5 9.81 m/s2), cd 5 drag coef” –
cient (kg/m), and m 5 mass (kg). Note that the drag coef” cient is
related to more fundamental parameters by
cd 5 1
where r 5 air density (kg/m3), A 5 area (m2), and Cd 5 the di-
mensionless drag coef” cient. Use the following parameter values
for your calculation: d 5 22 cm, m 5 0.4 kg, r 5 1.3 kg/m3, and
Cd 5 0.52.
25.26 Three linked bungee jumpers are depicted in Fig. P25.26. If
the bungee cords are idealized as linear springs (i.e., governed by
Hooke’s law), the following differential equations based on force
balances can be developed
m1 d 2×1
dt 2 5 m1g 1 k2(x2 2 x1) 2 k1x1
m2 d 2×2
dt 2 5 m2g 1 k3(x3 2 x2) 1 k2(x1 2 x2)
m3 d 2×3
dt 2 5 m3g 1 k3(x2 2 x3)
where mi 5 the mass of jumper i (kg), kj 5 the spring constant for
cord j (N/m), xi 5 the displacement of jumper i measured downward
FIGURE P25.26 Three individuals connected by bungee cords.
x1 = 0
(a) Unstretched (b) Stretched
x2 = 0
x3 = 0
26 C H A P T E R 26
Stiffness and Multistep Methods
This chapter covers two areas. First, we describe stiff ODEs. These are both indi-
vidual and systems of ODEs that have both fast and slow components to their solution.
We introduce the idea of an implicit solution technique as one commonly used remedy
for this problem. Then we discuss multistep methods. These algorithms retain informa-
tion of previous steps to more effectively capture the trajectory of the solution. They
also yield the truncation error estimates that can be used to implement adaptive step-
Stiffness is a special problem that can arise in the solution of ordinary differential equa-
tions. A stiff system is one involving rapidly changing components together with slowly
changing ones. In many cases, the rapidly varying components are ephemeral transients
that die away quickly, after which the solution becomes dominated by the slowly varying
components. Although the transient phenomena exist for only a short part of the integra-
tion interval, they can dictate the time step for the entire solution.
Both individual and systems of ODEs can be stiff. An example of a single stiff
dt 5 21000y 1 3000 2 2000e2t (26.1)
If y(0) 5 0, the analytical solution can be developed as
y 5 3 2 0.998e21000t 2 2.002e2t (26.2)
As in Fig. 26.1, the solution is initially dominated by the fast exponential term
(e21000t). After a short period (t , 0.005), this transient dies out and the solution becomes
dictated by the slow exponential (e2t).
Insight into the step size required for stability of such a solution can be gained by
examining the homogeneous part of Eq. (26.1),
dt 5 2ay (26.3)
756 STIFFNESS AND MULTISTEP METHODS
If y(0) 5 y0, calculus can be used to determine the solution as
y 5 y0e 2at
Thus, the solution starts at y0 and asymptotically approaches zero.
Euler’s method can be used to solve the same problem numerically:
yi11 5 yi 1 dyi
Substituting Eq. (26.3) gives
yi11 5 yi 2 ayih
yi11 5 yi(1 2 ah) (26.4)
The stability of this formula clearly depends on the step size h. That is, 01 2 ah 0 must be less than 1. Thus, if h . 2ya, 0 yi 0 n q as i n q. For the fast transient part of Eq. (26.2), this criterion can be used to show that the step
size to maintain stability must be , 2y1000 5 0.002. In addition, it should be noted that, whereas this criterion maintains stability (that is, a bounded solution), an even smaller step
size would be required to obtain an accurate solution. Thus, although the transient occurs for
only a small fraction of the integration interval, it controls the maximum allowable step size.
Super! cially, you might suppose that the adaptive step-size routines described at the
end of the last chapter might offer a solution for this dilemma. You might think that they
would use small steps during the rapid transients and large steps otherwise. However,
this is not the case, because the stability requirement will still necessitate using very
small steps throughout the entire solution.
FIGURE 26.1 Plot of a stiff solution of a single ODE. Although the solution appears to start at 1, there is actually a fast transient from y 5 0 to 1 that occurs in less than 0.005 time unit. This transient is perceptible only when the response is viewed on the fi ner timescale in the inset.
0 42 t0
26.1 STIFFNESS 757
Rather than using explicit approaches, implicit methods offer an alternative remedy.
Such representations are called implicit because the unknown appears on both sides of
the equation. An implicit form of Euler’s method can be developed by evaluating the
derivative at the future time,
yi11 5 yi 1 dyi11
This is called the backward, or implicit, Euler’s method. Substituting Eq. (26.3) yields
yi11 5 yi 2 ayi11 h
which can be solved for
yi11 5 yi
1 1 ah (26.5)
For this case, regardless of the size of the step, 0 yi 0 n 0 as i n q. Hence, the approach is called unconditionally stable.
EXAMPLE 26.1 Explicit and Implicit Euler
Problem Statement. Use both the explicit and implicit Euler methods to solve
dt 5 21000y 1 3000 2 2000e2t
where y(0) 5 0. (a) Use the explicit Euler with step sizes of 0.0005 and 0.0015 to solve
for y between t 5 0 and 0.006. (b) Use the implicit Euler with a step size of 0.05 to
solve for y between 0 and 0.4.
(a) For this problem, the explicit Euler’s method is
yi11 5 yi 1 (21000yi 1 3000 2 2000e 2ti)h
The result for h 5 0.0005 is displayed in Fig. 26.2a along with the analytical solu-
tion. Although it exhibits some truncation error, the result captures the general shape
of the analytical solution. In contrast, when the step size is increased to a value just
below the stability limit (h 5 0.0015), the solution manifests oscillations. Using
h . 0.002 would result in a totally unstable solution, that is, it would go in! nite
as the solution progressed.
(b) The implicit Euler’s method is
yi11 5 yi 1 (21000yi11 1 3000 2 2000e 2ti11)h
Now because the ODE is linear, we can rearrange this equation so that yi11 is isolated
on the left-hand side,
yi11 5 yi 1 3000h 2 2000he
1 1 1000h
The result for h 5 0.05 is displayed in Fig. 26.2b along with the analytical solution.
Notice that even though we have used a much bigger step size than the one that
758 STIFFNESS AND MULTISTEP METHODS
induced instability for the explicit Euler, the numerical solution tracks nicely on
the analytical result.
FIGURE 26.2 Solution of a “stiff” ODE with (a) the explicit and (b) implicit Euler methods.
h = 0.0015
h = 0.0005
h = 0.05
Systems of ODEs can also be stiff. An example is
dt 5 25y1 1 3y2 (26.6a)
dt 5 100y1 2 301y2 (26.6b)
For the initial conditions y1(0) 5 52.29 and y2(0) 5 83.82, the exact solution is
y1 5 52.96e 23.9899t
2 0.67e2302.0101t (26.7a)
y2 5 17.83e 23.9899t
1 65.99e2302.0101t (26.7b)
Note that the exponents are negative and differ by about 2 orders of magnitude. As with
the single equation, it is the large exponents that respond rapidly and are at the heart of
the system’s stiffness.
26.2 MULTISTEP METHODS 759
An implicit Euler’s method for systems can be formulated for the present example as
y1, i11 5 y1, i 1 (25y1, i11 1 3y2, i11)h (26.8a)
y2, i11 5 y2, i 1 (100y1, i11 2 301y2, i11)h (26.8b)
Collecting terms gives
(1 1 5h)y1, i11 2 3hy2, i11 5 y1, i (26.9a)
2100hy1, i11 1 (1 1 301h)y2, i11 5 y2, i (26.9b)
Thus, we can see that the problem consists of solving a set of simultaneous equations
for each time step.
For nonlinear ODEs, the solution becomes even more dif! cult since it involves
solving a system of nonlinear simultaneous equations (recall Sec. 6.6). Thus, although
stability is gained through implicit approaches, a price is paid in the form of added solu-
The implicit Euler method is unconditionally stable and only ! rst-order accurate. It
is also possible to develop in a similar manner a second-order accurate implicit trapezoi-
dal rule integration scheme for stiff systems. It is usually desirable to have higher-order
methods. The Adams-Moulton formulas described later in this chapter can also be used
to devise higher-order implicit methods. However, the stability limits of such approaches
are very stringent when applied to stiff systems. Gear (1971) developed a special series
of implicit schemes that have much larger stability limits based on backward difference
formulas. Extensive efforts have been made to develop software to ef! ciently implement
Gear’s methods. As a result, this is probably the most widely used method to solve stiff
systems. In addition, Rosenbrock and others (see Press et al., 2007) have proposed
implicit Runge-Kutta algorithms where the k terms appear implicitly. These methods have
good stability characteristics and are quite suitable for solving systems of stiff ordinary
26.2 MULTISTEP METHODS
The one-step methods described in the previous sections utilize information at a single
point xi to predict a value of the dependent variable yi11 at a future point xi11 (Fig. 26.3a).
Alternative approaches, called multistep methods (Fig. 26.3b), are based on the insight
that, once the computation has begun, valuable information from previous points is at
our command. The curvature of the lines connecting these previous values provides
information regarding the trajectory of the solution. The multistep methods explored in
this chapter exploit this information to solve ODEs. Before describing the higher-order
versions, we will present a simple second-order method that serves to demonstrate the
general characteristics of multistep approaches.
26.2.1 The Non-Self-Starting Heun Method
Recall that the Heun approach uses Euler’s method as a predictor [Eq. (25.15)]:
y0i11 5 yi 1 f(xi, yi)h (26.10)
760 STIFFNESS AND MULTISTEP METHODS
and the trapezoidal rule as a corrector [Eq. (25.16)]:
yi11 5 yi 1 f(xi, yi) 1 f(xi11, y
2 h (26.11)
Thus, the predictor and the corrector have local truncation errors of O(h2) and O(h3),
respectively. This suggests that the predictor is the weak link in the method because it
has the greatest error. This weakness is signi! cant because the ef! ciency of the iterative
corrector step depends on the accuracy of the initial prediction. Consequently, one way
to improve Heun’s method is to develop a predictor that has a local error of O(h3). This
can be accomplished by using Euler’s method and the slope at yi, and extra information
from a previous point yi21, as in
y0i11 5 yi21 1 f(xi, yi)2h (26.12)
Notice that Eq. (26.12) attains O(h3) at the expense of employing a larger step size, 2h. In
addition, note that Eq. (26.12) is not self-starting because it involves a previous value of the
dependent variable yi 2 1. Such a value would not be available in a typical initial-value problem.
Because of this fact, Eqs. (26.11) and (26.12) are called the non-self-starting Heun method.
As depicted in Fig. 26.4, the derivative estimate in Eq. (26.12) is now located at the
midpoint rather than at the beginning of the interval over which the prediction is made.
As demonstrated subsequently, this centering improves the error of the predictor to O(h3).
However, before proceeding to a formal derivation of the non-self-starting Heun, we will
summarize the method and express it using a slightly modi! ed nomenclature:
Predictor: y0i11 5 y m i21 1 f(xi, y
m i )2h (26.13)
Corrector: y j i11 5 y
m i 1
f(xi, y m i ) 1 f(xi11, y
j21 i11 )
(for j 5 1, 2, p , m) (26.14)
FIGURE 26.3 Graphical depiction of the fundamental difference between (a) one-step and (b) multistep methods for solving ODEs.
xxi + 1
xxi + 1xi – 1xi – 2
26.2 MULTISTEP METHODS 761
where the superscripts have been added to denote that the corrector is applied iteratively
from j 5 1 to m to obtain re! ned solutions. Note that ymi and y m i21 are the ! nal results
of the corrector iterations at the previous time steps. The iterations are terminated at any
time step on the basis of the stopping criterion
Zea Z 5 ` y j i11 2 y
y j i11
` 100% (26.15) When ea is less than a prespeci! ed error tolerance es, the iterations are terminated. At this
point, j 5 m. The use of Eqs. (26.13) through (26.15) to solve an ODE is demonstrated
in the following example.
EXAMPLE 26.2 Non-Self-Starting Heun Method
Problem Statement. Use the non-self-starting Heun method to perform the same com- putations as were performed previously in Example 25.5 using Heun’s method. That is,
FIGURE 26.4 A graphical depiction of the non-self-starting Heun method. (a) The midpoint method that is used as a predictor. (b) The trapezoidal rule that is employed as a corrector.
Slope = f (xi+1, yi+1) 0
Slope = f (xi, yi) + f (xi+1, yi+1)
762 STIFFNESS AND MULTISTEP METHODS
integrate y9 5 4e0.8x 2 0.5y from x 5 0 to x 5 4 using a step size of 1.0. As with Example
25.5, the initial condition at x 5 0 is y 5 2. However, because we are now dealing with a
multistep method, we require the additional information that y is equal to 20.3929953 at
x 5 21.
Solution. The predictor [Eq. (26.13)] is used to extrapolate linearly from x 5 21 to x 5 1.
y01 5 20.3929953 1 [4e 0.8(0)
2 0.5(2)] 2 5 5.607005
The corrector [Eq. (26.14)] is then used to compute the value:
y11 5 2 1 4e0.8(0) 2 0.5(2) 1 4e0.8(1) 2 0.5(5.607005)
2 1 5 6.549331
which represents a percent relative error of 25.73 percent (true value 5 6.194631). This
error is somewhat smaller than the value of 28.18 percent incurred in the self-starting Heun.
Now, Eq. (26.14) can be applied iteratively to improve the solution:
y21 5 2 1 3 1 4e0.8(1) 2 0.5(6.549331)
2 1 5 6.313749
which represents an et of 21.92%. An approximate estimate of the error can also be
determined using Eq. (26.15):
0 ea 0 5 ` 6.313749 2 6.549331 6.313749
` 100% 5 3.7% Equation (26.14) can be applied iteratively until ea falls below a prespeci! ed value of
es. As was the case with the Heun method (recall Example 25.5), the iterations converge
on a value of 6.360865 (et 5 22.68%). However, because the initial predictor value is
more accurate, the multistep method converges at a somewhat faster rate.
For the second step, the predictor is
y02 5 2 1 [4e 0.8(1)
2 0.5(6.360865) ] 2 5 13.44346 et 5 9.43%
which is superior to the prediction of 12.08260 (et 5 18%) that was computed with the
original Heun method. The ! rst corrector yields 15.76693 (et 5 6.8%), and subsequent
iterations converge on the same result as was obtained with the self-starting Heun method:
15.30224 (et 5 23.1%). As with the previous step, the rate of convergence of the corrector
is somewhat improved because of the better initial prediction.
Derivation and Error Analysis of Predictor-Corrector Formulas. We have just em- ployed graphical concepts to derive the non-self-starting Heun. We will now show how
the same equations can be derived mathematically. This derivation is particularly interest-
ing because it ties together ideas from curve ! tting, numerical integration, and ODEs.
The exercise is also useful because it provides a simple procedure for developing higher-
order multistep methods and estimating their errors.
The derivation is based on solving the general ODE
dx 5 f(x, y)
26.2 MULTISTEP METHODS 763
This equation can be solved by multiplying both sides by dx and integrating between
limits at i and i 1 1:
dy 5 #
f(x, y) dx
The left side can be integrated and evaluated using [recall Eq. (25.21)]:
yi11 5 yi 1 # xi11
f(x, y) dx (26.16)
Equation (26.16) represents a solution to the ODE if the integral can be evaluated.
That is, it provides a means to compute a new value of the dependent variable yi11 on
the basis of a prior value yi and the differential equation.
Numerical integration formulas such as those developed in Chap. 21 provide one
way to make this evaluation. For example, the trapezoidal rule [Eq. (21.3)] can be used
to evaluate the integral, as in
f(x, y) dx 5
f(xi, yi) 1 f(xi11, yi11)
2 h (26.17)
where h 5 xi1 1 2 xi is the step size. Substituting Eq. (26.17) into Eq. (26.16) yields
yi11 5 yi 1 f(xi, yi) 1 f(xi11, yi11)
which is the corrector equation for the Heun method. Because this equation is based on
the trapezoidal rule, the truncation error can be taken directly from Table 21.2,
Ec 5 2 1
12 h3y(3)(jc) 5 2