# Engineering

FIGURE P25.20 A spherical tank.

H

r

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:

d 2y

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:

dy

dt 5 y

dv

dt 5 2g 2

cd

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

2 rACd

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

755

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-

size control.

26.1 STIFFNESS

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

ODE is

dy

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),

dy

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

dt h

Substituting Eq. (26.3) gives

yi11 5 yi 2 ayih

or

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.

3

y

2

1

0 42 t0

1

0 0.020.010

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

dt h

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

dy

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.

Solution.

(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

2ti11

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.

1.5

y

1

0.5

0 0.0060.004

h = 0.0015

h = 0.0005

Exact

(a)

t0 0.002

2

y

1

0 0.40.3

Exact

h = 0.05

(b)

t0 0.20.1

Systems of ODEs can also be stiff. An example is

dy1

dt 5 25y1 1 3y2 (26.6a)

dy2

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-

tion complexity.

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

differential equations.

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

0 i11)

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 )

2 h

(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.

y

xi

(a)

xxi + 1

y

xi

(b)

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

j21 i11

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.

y

xxi+1

xi–1

xi

(a)

(b)

Slope = f (xi+1, yi+1) 0

y

xxi+1xi

Slope = f (xi, yi) + f (xi+1, yi+1)

2

0

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

dy

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:

# yi11

yi

dy 5 #

xi11

xi

f(x, y) dx

The left side can be integrated and evaluated using [recall Eq. (25.21)]:

yi11 5 yi 1 # xi11

xi

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

# xi11

xi

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)

2 h

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

1