INTRODUCTION

In Section 2.6 we examined one of the simplest numerical methods for approximating solutions of first-order initial-value problems y′ = f(x, y), y(x₀) = y₀. Recall that the backbone of Euler’s method was the formula

(1)

where f is the function obtained from the differential equation y′ = f(x, y). The recursive use of (1) for n = 0, 1, 2, . . . yields the y-coordinates y₁, y₂, y₃, . . . of points on successive “tangent lines” to the solution curve at x₁, x₂, x₃, . . . or x_n = x₀ + nh, where h is a constant and is the size of the step between x_n and x_{n + 1}. The values y₁, y₂, y₃, . . . approximate the values of a solution y(x) of the IVP at x₁, x₂, x₃, . . . . But whatever advantage (1) has in its simplicity is lost in the crudeness of its approximations.

A Comparison

In Problem 4 in Exercises 2.6 you were asked to use Euler’s method to obtain the approximate value of y(1.5) for the solution of the initial-value problem y′ = 2xy, y(1) = 1. You should have obtained the analytic solution and results similar to those given in Tables 6.1.1 and 6.1.2.

In this case, with a step size h = 0.1, a 16% relative error in the calculation of the approximation to y(1.5) is totally unacceptable. At the expense of doubling the number of calculations, some improvement in accuracy is obtained by halving the step size to h = 0.05.

Errors in Numerical Methods

In choosing and using a numerical method for the solution of an initial-value problem, we must be aware of the various sources of errors. For some kinds of computation the accumulation of errors might reduce the accuracy of an approximation to the point of making the computation useless. On the other hand, depending on the use to which a numerical solution may be put, extreme accuracy may not be worth the added expense and complication.

One source of error always present in calculations is round-off error. This error results from the fact that any calculator or computer can represent numbers using only a finite number of digits. Suppose, for the sake of illustration, that we have a calculator that uses base 10 arithmetic and carries four digits, so that is represented in the calculator as 0.3333 and is represented as 0.1111. If we use this calculator to compute (x² − )/(x − ) for x = 0.3334, we obtain

With the help of a little algebra, however, we see that

so that when x = 0.3334, (x² − )/(x − ) ≈ 0.3334 + 0.3333 = 0.6667. This example shows that the effects of round-off error can be quite serious unless some care is taken. One way to reduce the effect of round-off error is to minimize the number of calculations. Another technique on a computer is to use double-precision arithmetic to check the results. In general, round-off error is unpredictable and difficult to analyze, and we will neglect it in the error analysis that follows. We will concentrate on investigating the error introduced by using a formula or algorithm to approximate the values of the solution.

Truncation Errors for Euler’s Method

In the sequence of values y₁, y₂, y₃, . . . generated from (1), usually the value of y₁ will not agree with the actual solution evaluated at x₁, namely, y(x₁), because the algorithm gives only a straight-line approximation to the solution. See Figure 2.6.2. The error is called the local truncation error, formula error, or discretization error. It occurs at each step; that is, if we assume that y_n is accurate, then y_{n + 1} will contain a local truncation error.

To derive a formula for the local truncation error for Euler’s method, we use Taylor’s formula with remainder. If a function y(x) possesses k + 1 derivatives that are continuous on an open interval containing a and x, then

where c is some point between a and x. Setting k = 1, a = x_n, and x = x_{n + 1} = x_n + h, we get

.

Euler’s method (1) is the last formula without the last term; hence the local truncation error in y_{n + 1} is

where x_n < c < x_{n + 1}.

The value of c is usually unknown (it exists theoretically), and so the exact error cannot be calculated, but an upper bound on the absolute value of the error is

In discussing errors arising from the use of numerical methods, it is helpful to use the notation O(hⁿ). To define this concept we let e(h) denote the error in a numerical calculation depending on h. Then e(h) is said to be of order hⁿ, denoted by O(hⁿ), if there is a constant C and a positive integer n such that |e(h)| ≤ Chⁿ for h is sufficiently small. Thus the local truncation error for Euler’s method is O(h²). We note that, in general, if e(h) in a numerical method is of order hⁿ, and h is halved, the new error is approximately C(h/2)ⁿ = Chⁿ/2ⁿ; that is, the error is reduced by a factor of ⁿ.

EXAMPLE 1 Bound for Local Truncation Error

Note in Example 1 that if h is halved to 0.05, then the error bound is 0.0480, about one-fourth as much as shown in (2). This is expected because the local truncation error for Euler’s method is O(h²).

In the above analysis we assumed that the value of y_n was exact in the calculation of y_{n + 1}, but it is not because it contains local truncation errors from previous steps. The total error in y_{n + 1} is an accumulation of the errors in each of the previous steps. This total error is called the global truncation error. A complete analysis of the global truncation error is beyond the scope of this text, but it can be shown that the global truncation error for Euler’s method is O(h).

We expect that, for Euler’s method, if the step size is halved the error will be approximately halved as well. This is borne out in Tables 6.1.1 and 6.1.2, where the absolute error at x = 1.50 with h = 0.1 is 0.5625 and with h = 0.05 is 0.3171, approximately half as large.

In general it can be shown that if a method for the numerical solution of a differential equation has local truncation error O(h^{α + 1}), then the global truncation error is O(h^α).

For the remainder of this section and in the subsequent sections we study methods that give significantly greater accuracy than Euler’s method.

Improved Euler’s Method

The numerical method defined by the formula

,(3)

where (4)

is commonly known as the improved Euler’s method. In order to compute y_n+1 for n = 0, 1, 2, . . . , from (3) we must, at each step, first use Euler’s method (4) to obtain an initial estimate . For example, with n = 0, (4) would give = y₀ + hf(x₀, y₀), and then knowing this value we use (3) to get y₁ = y₀ + h(f(x₀, y₀) + f(x₁, ))/2, where x₁ = x₀ + h. These equations can be readily visualized. In FIGURE 6.1.1 observe that m₀ = f(x₀, y₀) and m₁ = f(x₁, ) are slopes of the solid straight lines shown passing through the points (x₀, y₀) and (x₁, ), respectively. By taking an average of these slopes, that is, m_ave = (f(x₀, y₀) + f(x₁, ))/2, we obtain the slope of the parallel dashed skew lines. With the first step, rather than advancing along the line through (x₀, y₀) with slope f(x₀, y₀) to the point with y-coordinate obtained by Euler’s method, we advance instead along the red dashed line through (x₀, y₀) with slope m_ave until we reach x₁. It seems plausible from inspection of the figure that y₁ is an improvement over .

In general, the improved Euler’s method is an example of a predictor–corrector method. The value of given by (4) predicts a value of y(x_n), whereas the value of y_{n + 1} defined by formula (3) corrects this estimate.

EXAMPLE 2 Improved Euler’s Method

Use the improved Euler’s method to obtain the approximate value of y(1.5) for the solution of the initial-value problem y′ = 2xy, y(1) = 1. Compare the results for h = 0.1 and h = 0.05.

SOLUTION

With x₀ = 1, y₀ = 1, f(x_n , y_n) = 2x_ny_n, n = 0, and h = 0.1 we first compute (4):

= y₀ + (0.1)(2x₀ y₀) = 1 + (0.1)2(1)(1) = 1.2.

We use this last value in (3) along with x₁ = 1 + h = 1 + 0.1 = 1.1:

The comparative values of the calculations for h = 0.1 and h = 0.05 are given in Tables 6.1.3 and 6.1.4, respectively.

A brief word of caution is in order here. We cannot compute all the values of first and then substitute these values into formula (3). In other words, we cannot use the data in Table 6.1.1 to help construct the values in Table 6.1.3. Why not?

Truncation Errors for the Improved Euler’s Method

The local truncation error for the improved Euler’s method is O(h³). The derivation of this result is similar to the derivation of the local truncation error for Euler’s method. Since the local truncation error for the improved Euler’s method is O(h³), the global truncation error is O(h²). This can be seen in Example 2; when the step size is halved from h = 0.1 to h = 0.05, the absolute error at x = 1.50 is reduced from 0.0394 to 0.0108, a reduction of approximately ()² = .

6.1 Exercises Answers to selected odd-numbered problems begin on page ANS-14.

Given the initial-value problems in Problems 1–10, use the improved Euler’s method to obtain a four-decimal approximation to the indicated value. First use h = 0.1 and then use h = 0.05.

1. y′ = 2x − 3y + 1,     y(1) = 5; y(1.5)
2. y′ = 4x − 2y,     y(0) = 2; y(0.5)
3. y′ = 1 + y²,     y(0) = 0; y(0.5)
4. y′ = x² + y²,     y(0) = 1; y(0.5)
5. y′ = e^2y,     y(0) = 0; y(0.5)
6. y′ = x + y²,     y(0) = 0; y(0.5)
7. y′ = (x² y)²,     y(0) = 0.5; y(0.5)
8. y′ = xy + ,     y(0) = 1; y(0.5)
9. y′ = xy² − ,     y(1) = 1; y(1.5)
10. y′ = y − y²,     y(0) = 0.5; y(0.5)
11. Consider the initial-value problem y′ = (x + y − 1)², y(0) = 2. Use the improved Euler’s method with h = 0.1 and h = 0.05 to obtain approximate values of the solution at x = 0.5. At each step compare the approximate value with the exact value of the analytic solution.
12. Although it may not be obvious from the differential equation, its solution could “behave badly” near a point x at which we wish to approximate y(x). Numerical procedures may give widely differing results near this point. Let y(x) be the solution of the initial-value problem y′ = x² + y³, y(1) = 1.
    1. Use a numerical solver to obtain the graph of the solution on the interval [1, 1.4].
    2. Using the step size h = 0.1, compare the results obtained from Euler’s method with the results from the improved Euler’s method in the approximation of y(1.4).
13. Consider the initial-value problem y′ = 2y, y(0) = 1. The analytic solution is y = e^2x.
    1. Approximate y(0.1) using one step and Euler’s method.
    2. Find a bound for the local truncation error in y₁.
    3. Compare the actual error in y₁ with your error bound.
    4. Approximate y(0.1) using two steps and Euler’s method.
    5. Verify that the global truncation error for Euler’s method is O(h) by comparing the errors in parts (a) and (d).
14. Repeat Problem 13 using the improved Euler’s method. Its global truncation error is O(h²).
15. Repeat Problem 13 using the initial-value problem y′ = x − 2y, y(0) = 1. The analytic solution is

    y = x − + e^−2x.

16. Repeat Problem 15 using the improved Euler’s method. Its global truncation error is O(h²).
17. Consider the initial-value problem y′ = 2x − 3y + 1, y(1) = 5. The analytic solution is

    y(x) = + x + e^{−3(x − 1)}.

    1. Find a formula involving c and h for the local truncation error in the nth step if Euler’s method is used.
    2. Find a bound for the local truncation error in each step if h = 0.1 is used to approximate y(1.5).
    3. Approximate y(1.5) using h = 0.1 and h = 0.05 with Euler’s method. See Problem 1 in Exercises 2.6.
    4. Calculate the errors in part (c) and verify that the global truncation error of Euler’s method is O(h).
18. Repeat Problem 17 using the improved Euler’s method, which has a global truncation error O(h²). See Problem 1. You may need to keep more than four decimal places to see the effect of reducing the order of error.
19. Repeat Problem 17 for the initial-value problem y′ = e^−y, y(0) = 0. The analytic solution is y(x) = ln(x + 1). Approximate y(0.5). See Problem 5 in Exercises 2.6.
20. Repeat Problem 19 using the improved Euler’s method, which has a global truncation error O(h²). See Problem 5. You may need to keep more than four decimal places to see the effect of reducing the order of error.

Discussion Problem

21. Answer the question: “Why not?” that follows the three sentences after Example 2 on page 310.