INTRODUCTION

In Section 2.1 we saw that we could glean qualitative information from a first-order DE about its solutions even before we attempted to solve the equation. In Sections 2.2–2.5 we examined first-order DEs analytically; that is, we developed procedures for actually obtaining explicit and implicit solutions. But many differential equations possess solutions and yet these solutions cannot be obtained analytically. In this case we “solve” the differential equation numerically; this means that the DE is used as the cornerstone of an algorithm for approximating the unknown solution. It is common practice to refer to the algorithm as a numerical method, the approximate solution as a numerical solution, and the graph of a numerical solution as a numerical solution curve.

In this section we are going to consider only the simplest of numerical methods. A more extensive treatment of this subject is found in Chapter 6.

Using the Tangent Line

Let us assume that the first-order initial-value problem

y′ = f(x, y),     y(x₀) = y₀(1)

possesses a solution. One of the simplest techniques for approximating this solution is to use tangent lines. For example, let y(x) denote the unknown solution of the first-order initial-value problem , y(2) = 4. The nonlinear differential equation cannot be solved directly by the methods considered in Sections 2.2, 2.4, and 2.5; nevertheless, we can still find approximate numerical values of the unknown y(x). Specifically, suppose we wish to know the value of y(2.5). The IVP has a solution, and, as the flow of the direction field in FIGURE 2.6.1(a) suggests, a solution curve must have a shape similar to the curve shown in blue.

The direction field in Figure 2.6.1(a) was generated so that the lineal elements pass through points in a grid with integer coordinates. As the solution curve passes through the initial point (2, 4), the lineal element at this point is a tangent line with slope given by f(2, 4) = 0.1 + 0.4(2)² = 1.8. As is apparent in Figure 2.6.1(a) and the “zoom in” in Figure 2.6.1(b), when x is close to 2 the points on the solution curve are close to the points on the tangent line (the lineal element). Using the point (2, 4), the slope f(2, 4) = 1.8, and the point-slope form of a line, we find that an equation of the tangent line is y = L(x), where . This last equation, called a linearization of y(x) at x = 2, can be used to approximate values y(x) within a small neighborhood of x = 2. If y₁ = L(x₁) denotes the value of the y-coordinate on the tangent line and y(x₁) is the y-coordinate on the solution curve corresponding to an x-coordinate x₁ that is close to x = 2, then y(x₁) ≈ y₁. If we choose, say, x₁ = 2.1, then y₁ = L(2.1) = 1.8(2.1) + 0.4 = 4.18, and so y(2.1) ≈ 4.18.

Euler’s Method

To generalize the procedure just illustrated, we use the linearization of the unknown solution y(x) of (1) at x = x₀:

L(x) = f(x₀, y₀)(x − x₀) + y₀.(2)

The graph of this linearization is a straight line tangent to the graph of y = y(x) at the point (x₀, y₀). We now let h be a positive increment of the x-axis, as shown in FIGURE 2.6.2. Then by replacing x by x₁ = x₀ + h in (2) we get

L(x₁) = f(x₀, y₀)(x₀ + h − x₀) + y₀     or     y₁ = y₀ + hf(x₀, y₀),

where y₁ = L(x₁). The point (x₁, y₁) on the tangent line is an approximation to the point (x₁, y(x₁)) on the solution curve. Of course the accuracy of the approximation y₁ ≈ y(x₁) depends heavily on the size of the increment h. Usually we must choose this step size to be “reasonably small.” We now repeat the process using a second “tangent line” at (x₁, y₁).^* By replacing (x₀, y₀) in the above discussion with the new starting point (x₁, y₁), we obtain an approximation y₂ ≈ y(x₂) corresponding to two steps of length h from x₀, that is, x₂ = x₁ + h = x₀ + 2h and

y(x₂) = y(x₀ + 2h) = y(x₁ + h) ≈ y₂ = y₁ + hf(x₁, y₁).

Continuing in this manner, we see that y₁, y₂, y₃, . . ., can be defined recursively by the general formula

(3)

where x_n = x₀ + nh, n = 0, 1, 2, . . . . This procedure of using successive “tangent lines” is called Euler’s method.

EXAMPLE 1 Euler’s Method

In Example 2 we apply Euler’s method to a differential equation for which we have already found a solution. We do this to compare the values of the approximations y_n at each step with the true values of the solution y(x_n) of the initial-value problem.

EXAMPLE 2 Comparison of Approximate and Exact Values

Consider the initial-value problem y′ = 0.2xy, y(1) = 1. Use Euler’s method to obtain an approximation to y(1.5) using first h = 0.1 and then h = 0.05.

SOLUTION

With the identification f(x, y) = 0.2xy, (3) becomes

y_{n + 1} = y_n + h(0.2x_n y_n),

where x₀ = 1 and y₀ = 1. Again with the aid of computer software we obtain the values in Tables 2.6.3 and 2.6.4.

TABLE 2.6.3 h = 0.1

TABLE 2.6.4 h = 0.05

≡

In Example 1, the true values were calculated from the known solution (verify). Also, the absolute error is defined to be

The relative error and percentage relative error are, in turn,

    and    

By comparing the last two columns in Tables 2.6.3 and 2.6.4, it is clear that the accuracy of the approximations improve as the step size h decreases. Also, we see that even though the percentage relative error is growing with each step, it does not appear to be that bad. But you should not be deceived by one example. If we simply change the coefficient of the right side of the DE in Example 2 from 0.2 to 2, then at x_n = 1.5 the percentage relative errors increase dramatically. See Problem 4 in Exercises 2.6.

A Caveat

Euler’s method is just one of many different ways a solution of a differential equation can be approximated. Although attractive for its simplicity, Euler’s method is seldom used in serious calculations. We have introduced this topic simply to give you a first taste of numerical methods. We will go into greater detail and discuss methods that give significantly greater accuracy, notably the fourth-order Runge–Kutta method, in Chapter 6. We shall refer to this important numerical method as the RK4 method.

Numerical Solvers

Regardless of whether we can actually find an explicit or implicit solution, if a solution of a differential equation exists, it represents a smooth curve in the Cartesian plane. The basic idea behind any numerical method for ordinary differential equations is to somehow approximate the y-values of a solution for preselected values of x. We start at a specified initial point (x₀, y₀) on a solution curve and proceed to calculate in a step-by-step fashion a sequence of points (x₁, y₁), (x₂, y₂), . . ., (x_n, y_n) whose y-coordinates y_i approximate the y-coordinates y(x_i) of points (x₁, y(x₁)), (x₂, y(x₂)), . . ., (x_n, y(x_n)) that lie on the graph of the usually unknown solution y(x). By taking the x-coordinates close together (that is, for small values of h) and by joining the points (x₁, y₁), (x₂, y₂), . . ., (x_n, y_n) with short line segments, we obtain a polygonal curve that appears smooth and whose qualitative characteristics we hope are close to those of an actual solution curve. Drawing curves is something well suited to a computer. A computer program written to either implement a numerical method or to render a visual representation of an approximate solution curve fitting the numerical data produced by this method is referred to as a numerical solver. There are many different numerical solvers commercially available, either embedded in a larger software package such as a computer algebra system or as a stand-alone package. Some software packages simply plot the generated numerical approximations, whereas others generate both hard numerical data as well as the corresponding approximate or numerical solution curves. As an illustration of the connect-the-dots nature of the graphs produced by a numerical solver, the two red polygonal graphs in FIGURE 2.6.3 are numerical solution curves for the initial-value problem y′ = 0.2xy, y(0) = 1, on the interval [0, 4] obtained from Euler’s method and the RK4 method using the step size h = 1. The blue smooth curve is the graph of the exact solution of the IVP. Notice in Figure 2.6.3 that even with the ridiculously large step size of h = 1, the RK4 method produces the more believable “solution curve.” The numerical solution curve obtained from the RK4 method is indistinguishable from the actual solution curve on the interval [0, 4] when a more typical step size of h = 0.1 is used.

Using a Numerical Solver

Knowledge of the various numerical methods is not necessary in order to use a numerical solver. A solver usually requires that the differential equation be expressed in normal form dy/dx = f(x, y). Numerical solvers that generate only curves usually require that you supply f(x, y) and the initial data x₀ and y₀ and specify the desired numerical method. If the idea is to approximate the numerical value of y(a), then a solver may additionally require that you state a value for h, or, equivalently, require the number of steps that you want to take to get from x = x₀ to x = a. For example, if we want to approximate y(4) for the IVP illustrated in Figure 2.6.3, then, starting at x = 0, it takes four steps to reach x = 4 with a step size of h = 1; 40 steps is equivalent to a step size of h = 0.1. Although it is not our intention here to delve into the many problems that one can encounter when attempting to approximate mathematical quantities, you should be at least aware of the fact that a numerical solver may break down near certain points or give an incomplete or misleading picture when applied to some first-order differential equations in the normal form. FIGURE 2.6.4 illustrates the numerical solution curve obtained by applying Euler’s method to a certain first-order initial value problem dy/dx = f(x, y), y(0) = 1. Equivalent results were obtained using three different commercial numerical solvers, yet the graph is hardly a plausible solution curve. (Why?) There are several avenues of recourse when a numerical solver has difficulties; three of the more obvious are decrease the step size, use another numerical method, or try a different numerical solver.

2.6 Exercises Answers to selected odd-numbered problems begin on page ANS-3.

In Problems 1 and 2, use Euler’s method to obtain a four-decimal approximation of the indicated value. Carry out the recursion of (3) by hand, first using h = 0.1 and then using h = 0.05.

1. y′ = 2x − 3y + 1,     y(1) = 5; y(1.2)
2. y′ = x + y²,     y(0) = 0; y(0.2)

In Problems 3 and 4, use Euler’s method to obtain a four-decimal approximation of the indicated value. First use h = 0.1 and then use h = 0.05. Find an explicit solution for each initial-value problem and then construct tables similar to Tables 2.6.3 and 2.6.4.

3. y′ = y,     y(0) = 1; y(1.0)
4. y′ = 2xy,     y(1) = 1; y(1.5)

In Problems 5–10, use a numerical solver and Euler’s method to obtain a four-decimal approximation of the indicated value. First use h = 0.1 and then use h = 0.05.

5. y′ = e^-y,     y(0) = 0; y(0.5)
6. y′ = x² + y²,     y(0) = 1; 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)

In Problems 11 and 12, use a numerical solver to obtain a numerical solution curve for the given initial-value problem. First use Euler’s method and then the RK4 method. Use h = 0.25 in each case. Superimpose both solution curves on the same coordinate axes. If possible, use a different color for each curve. Repeat, using h = 0.1 and h = 0.05.

11. y′ = 2(cos x)y,     y(0) = 1
12. y′ = y(10 − 2y),     y(0) = 1

Discussion Problem

13. Use a numerical solver and Euler’s method to approximate y(1.0), where y(x) is the solution to y′ = 2xy², y(0) = 1. First use h = 0.1 and then h = 0.05. Repeat using the RK4 method. Discuss what might cause the approximations of y(1.0) to differ so greatly.

*This is not an actual tangent line since (x₁, y₁) lies on the first tangent and not on the solution curve.