# Engineering

FIGURE 25.3 Comparison of the true solution with a numerical solution using Euler’s method for the integral of y9 5 22×3 1 12×2 2 20x 1 8.5 from x 5 0 to x 5 4 with a step size of 0.5. The initial condition at x 5 0 is y 5 1.

y

4

0 x4

True solution

h = 0.5

20

25.1 EULER’S METHOD 713

2. Round-off errors caused by the limited numbers of signi! cant digits that can be

retained by a computer.

The truncation errors are composed of two parts. The ! rst is a local truncation error

that results from an application of the method in question over a single step. The second

is a propagated truncation error that results from the approximations produced during

the previous steps. The sum of the two is the total, or global truncation, error.

Insight into the magnitude and properties of the truncation error can be gained by

deriving Euler’s method directly from the Taylor series expansion. To do this, realize that

the differential equation being integrated will be of the general form

y¿ 5 f(x, y) (25.3)

where y9 5 dyydx and x and y are the independent and the dependent variables, respectively. If the solution—that is, the function describing the behavior of y—has continuous deriva-

tives, it can be represented by a Taylor series expansion about a starting value (xi, yi), as in

[recall Eq. (4.7)]

yi11 5 yi 1 y¿i h 1 y–i

2! h2 1 p 1

y(n)i

n! hn 1 Rn (25.4)

where h 5 xi11 2 xi and Rn 5 the remainder term, de! ned as

Rn 5 y(n11)(j)

(n 1 1)! hn11 (25.5)

where j lies somewhere in the interval from xi to xi11. An alternative form can be de-

veloped by substituting Eq. (25.3) into Eqs. (25.4) and (25.5) to yield

yi11 5 yi 1 f(xi, yi)h 1 f ¿(xi, yi)

2! h2 1 p 1

f (n21)(xi, yi)

n! hn 1 O(hn11) (25.6)

where O(hn11) speci! es that the local truncation error is proportional to the step size

raised to the (n 1 1)th power.

By comparing Eqs. (25.2) and (25.6), it can be seen that Euler’s method corresponds

to the Taylor series up to and including the term f(xi, yi)h. Additionally, the comparison

indicates that a truncation error occurs because we approximate the true solution using

a ! nite number of terms from the Taylor series. We thus truncate, or leave out, a part of

the true solution. For example, the truncation error in Euler’s method is attributable to

the remaining terms in the Taylor series expansion that were not included in Eq. (25.2).

Subtracting Eq. (25.2) from Eq. (25.6) yields

Et 5 f ¿(xi, yi)

2! h2 1 p 1 O(hn11) (25.7)

where Et 5 the true local truncation error. For suf! ciently small h, the errors in the terms

in Eq. (25.7) usually decrease as the order increases (recall Example 4.2 and the ac-

companying discussion), and the result is often represented as

Ea 5 f ¿(xi, yi)

2! h2 (25.8)

714 RUNGE-KUTTA METHODS

or

Ea 5 O(h 2) (25.9)

where Ea 5 the approximate local truncation error.

EXAMPLE 25.2 Taylor Series Estimate for the Error of Euler’s Method

Problem Statement. Use Eq. (25.7) to estimate the error of the ! rst step of Example 25.1. Also use it to determine the error due to each higher-order term of the Taylor series

expansion.

Solution. Because we are dealing with a polynomial, we can use the Taylor series to obtain exact estimates of the errors in Euler’s method. Equation (25.7) can be written as

Et 5 f ¿(xi, yi)

2! h2 1

f –(xi, yi)

3! h3 1

f (3)(xi, yi)

4! h4 (E25.2.1)

where f 9(xi, yi) 5 the ! rst derivative of the differential equation (that is, the second de-

rivative of the solution). For the present case, this is

f ¿(xi, yi) 5 26x 2

1 24x 2 20 (E25.2.2)

and f 0(xi, yi) is the second derivative of the ODE

f –(xi, yi) 5 212x 1 24 (E25.2.3)

and f (3)(xi, yi) is the third derivative of the ODE

f (3)(xi, yi) 5 212 (E25.2.4)

We can omit additional terms (that is, fourth derivatives and higher) from Eq. (E25.2.1)

because for this particular case they equal zero. It should be noted that for other func-

tions (for example, transcendental functions such as sinusoids or exponentials) this would

not necessarily be true, and higher-order terms would have nonzero values. However, for

the present case, Eqs. (E25.2.1) through (E25.2.4) completely de! ne the truncation error

for a single application of Euler’s method.

For example, the error due to truncation of the second-order term can be calculated as

Et, 2 5 26(0.0)2 1 24(0.0) 2 20

2 (0.5)2 5 22.5 (E25.2.5)

For the third-order term:

Et, 3 5 212(0.0) 1 24

6 (0.5)3 5 0.5

and the fourth-order term:

Et, 4 5 212

24 (0.5)4 5 20.03125

These three results can be added to yield the total truncation error:

Et 5 Et, 2 1 Et, 3 1 Et, 4 5 22.5 1 0.5 2 0.03125 5 22.03125

25.1 EULER’S METHOD 715

which is exactly the error that was incurred in the initial step of Example 25.1. Note

how Et, 2 . Et, 3 . Et, 4, which supports the approximation represented by Eq. (25.8).

As illustrated in Example 25.2, the Taylor series provides a means of quantifying

the error in Euler’s method. However, there are limitations associated with its use for

this purpose:

1. The Taylor series provides only an estimate of the local truncation error—that is, the

error created during a single step of the method. It does not provide a measure of the

propagated and, hence, the global truncation error. In Table 25.1, we have included

the local and global truncation errors for Example 25.1. The local error was computed

for each time step with Eq. (25.2) but using the true value of yi (the second column

of the table) to compute each yi1l rather than the approximate value (the third column),

as is done in the Euler method. As expected, the average absolute local truncation

error (25 percent) is less than the average global error (90 percent). The only reason

that we can make these exact error calculations is that we know the true value a

priori. Such would not be the case in an actual problem. Consequently, as discussed

below, you must usually apply techniques such as Euler’s method using a number of

different step sizes to obtain an indirect estimate of the errors involved.

2. As mentioned above, in actual problems we usually deal with functions that are more

complicated than simple polynomials. Consequently, the derivatives that are needed

to evaluate the Taylor series expansion would not always be easy to obtain.

Although these limitations preclude exact error analysis for most practical problems,

the Taylor series still provides valuable insight into the behavior of Euler’s method. Ac-

cording to Eq. (25.9), we see that the local error is proportional to the square of the step

size and the ! rst derivative of the differential equation. It can also be demonstrated that

the global truncation error is O(h), that is, it is proportional to the step size (Carnahan

et al., 1969). These observations lead to some useful conclusions:

1. The error can be reduced by decreasing the step size.

2. The method will provide error-free predictions if the underlying function (that is, the

solution of the differential equation) is linear, because for a straight line the second

derivative would be zero.

This latter conclusion makes intuitive sense because Euler’s method uses straight-line

segments to approximate the solution. Hence, Euler’s method is referred to as a rst-

order method.

It should also be noted that this general pattern holds for the higher-order one-step

methods described in the following pages. That is, an nth-order method will yield perfect

results if the underlying solution is an nth-order polynomial. Further, the local truncation

error will be O(hn11) and the global error O(hn).

EXAMPLE 25.3 Effect of Reduced Step Size on Euler’s Method

Problem Statement. Repeat the computation of Example 25.1 but use a step size of 0.25.

716 RUNGE-KUTTA METHODS

Solution. The computation is repeated, and the results are compiled in Fig. 25.4a. Halving the step size reduces the absolute value of the average global error to 40 percent

and the absolute value of the local error to 6.4 percent. This is compared to global and

local errors for Example 25.1 of 90 percent and 24.8 percent, respectively. Thus, as

expected, the local error is quartered and the global error is halved.

Also, notice how the local error changes sign for intermediate values along the range.

This is due primarily to the fact that the ! rst derivative of the differential equation is a

parabola that changes sign [recall Eq. (E25.2.2) and see Fig. 25.4b]. Because the local

error is proportional to this function, the net effect of the oscillation in sign is to keep the

global error from continuously growing as the calculation proceeds. Thus, from x 5 0 to

x 5 1.25, the local errors are all negative, and consequently, the global error increases

FIGURE 25.4 (a) Comparison of two numerical solutions with Euler’s method using step sizes of 0.5 and 0.25. (b) Comparison of true and estimated local truncation error for the case where the step size is 0.5. Note that the “estimated” error is based on Eq. (E25.2.5).

y

4

0 x4

True solution

h = 0.5

h = 0.25

2

(a)

0

y

– 0.5

0 x42

True

Estimated

(b)

25.1 EULER’S METHOD 717

over this interval. In the intermediate section of the range, positive local errors begin to

reduce the global error. Near the end, the process is reversed and the global error again

in” ates. If the local error continuously changes sign over the computation interval, the net

effect is usually to reduce the global error. However, where the local errors are of the

same sign, the numerical solution may diverge farther and farther from the true solution

as the computation proceeds. Such results are said to be unstable.

The effect of further step-size reductions on the global truncation error of Euler’s

method is illustrated in Fig. 25.5. This plot shows the absolute percent relative error at

x 5 5 as a function of step size for the problem we have been examining in Examples

25.1 through 25.3. Notice that even when h is reduced to 0.001, the error still exceeds

0.1 percent. Because this step size corresponds to 5000 steps to proceed from x 5 0 to

x 5 5, the plot suggests that a ! rst-order technique such as Euler’s method demands

great computational effort to obtain acceptable error levels. Later in this chapter, we

present higher-order techniques that attain much better accuracy for the same computa-

tional effort. However, it should be noted that, despite its inef! ciency, the simplicity of

Euler’s method makes it an extremely attractive option for many engineering problems.

Because it is very easy to program, the technique is particularly useful for quick analy-

ses. In the next section, a computer algorithm for Euler’s method is developed.

FIGURE 25.5 Effect of step size on the global truncation error of Euler’s method for the integral of y9 5 22×3 1 12×2 2 20x 1 8.5. The plot shows the absolute percent relative global error at x 5 5 as a function of step size.

1

10

100

0.1 0.1

Step size

A b

s o

lu te

p e rc

e n

t re

la ti

v e e

rr o

r

0.01 0.0011

50

Steps

500 50005

718 RUNGE-KUTTA METHODS

25.1.2 Algorithm for Euler’s Method

Algorithms for one-step techniques such as Euler’s method are extremely simple to

program. As speci! ed previously at the beginning of this chapter, all one-step methods

have the general form

New value 5 old value 1 slope 3 step size (25.10)

The only way in which the methods differ is in the calculation of the slope.

Suppose that you want to perform the simple calculation outlined in Table 25.1. That

is, you would like to use Euler’s method to integrate y9 5 22×3 1 12×2 2 20x 1 8.5,

with the initial condition that y 5 1 at x 5 0. You would like to integrate out to x 5 4

using a step size of 0.5, and display all the results. A simple pseudocode to accomplish

this task could be written as in Fig. 25.6.

Although this program will “do the job” of duplicating the results of Table 25.1, it

is not very well designed. First, and foremost, it is not very modular. Although this is

not very important for such a small program, it would be critical if we desired to mod-

ify and improve the algorithm.

Further, there are a number of issues related to the way we have set up the iterations.

For example, suppose that the step size were to be made very small to obtain better ac-

curacy. In such cases, because every computed value is displayed, the number of output

values might be very large. Further, the algorithm is predicated on the assumption that

the calculation interval is evenly divisible by the step size. Finally, the accumulation of

x in the line x 5 x 1 dx can be subject to quantizing errors of the sort previously dis-

FIGURE 25.6 Pseudocode for a “dumb” version of Euler’s method.

‘set integration range

xi 5 0

xf 5 4

‘initialize variables

x 5 xi

y 5 1

‘set step size and determine

‘number of calculation steps

dx 5 0.5

nc 5 (xf 2 xi)/dx

‘output initial condition

PRINT x, y

‘loop to implement Euler’s method

‘and display results

DOFOR i 5 1, nc

dydx 5 22×3 1 12×2 2 20x 1 8.5

y 5 y 1 dydx ? dx

x 5 x 1 dx

PRINT x, y

END DO

25.1 EULER’S METHOD 719

cussed in Sec. 3.4.1. For example, if dx were changed to 0.01 and standard IEEE ” oat-

ing point representation were used (about seven signi! cant digits), the result at the end

of the calculation would be 3.999997 rather than 4. For dx 5 0.001, it would be 3.999892!

A much more modular algorithm that avoids these dif! culties is displayed in Fig. 25.7.

The algorithm does not output all calculated values. Rather, the user speci! es an output

interval, xout, that dictates the interval at which calculated results are stored in arrays, xpm

and ypm. These values are stored in arrays so that they can be output in a variety of ways

after the computation is completed (for example, printed, graphed, or written to a ! le).

The Driver Program takes big output steps and calls an Integrator routine that takes

! ner calculation steps. Note that the loops controlling both large and small steps exit on

logical conditions. Thus, the intervals do not have to be evenly divisible by the step sizes.

The Integrator routine then calls an Euler routine that takes a single step with Euler’s

method. The Euler routine calls a Derivative routine that calculates the derivative value.

Whereas such modularization might seem like overkill for the present case, it will

greatly facilitate modifying the program in later sections. For example, although the

program in Fig. 25.7 is speci! cally designed to implement Euler’s method, the Euler

module is the only part that is method-speci! c. Thus, all that is required to apply this

algorithm to the other one-step methods is to modify this routine.