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.