Stability of a Numerical Solution of Differential Equations

In 1926 Milne [1] published a numerical method for the solution of ordinary differential equations. This method turns out to be unstable, as shown by Muhin [2], Hildebrand [3], Liniger [4], and others. Instability was not too serious in the day of desk calculators but is fatal in the modern era of high speed computers. The basic cause of the instability in this particular method is the use of Simpson's rule to perform the final integration. Simpson's rule integrates over two intervals, and under certain conditions can produce an error which alternates in sign from step to step and which increases in magnitude exponentially. It is the purpose of this paper to show that the occasional application of Newton's “three eighths” quadrature formula over three intervals can effectively damp out the unwanted oscillation without harm to the desired solution.

Let the given differential equation be dy/dx = ƒ(x, y), and let the step length for the independent variable x be denoted by h. The quantity s = h ∂ƒ/∂y plays a basic role in the analysis, for it may be shown that when Simpson's rule is used an error E₀ at x = x₀ is propagated through subsequent steps according to the second order difference equation (1 - s_n+1/3) E_n+1 - (4s_n/3) E_n - (1 + s_n-1/3) E_n-1 = 0. (See e.g., Hildebrand [3], p. 206, Milne [5], p. 68.)

While in general s is a variable, the special case where s is constant not only permits a simple analysis but also serves to explain the behavior in other cases. Cf. Hildebrand [3], p. 202. Accordingly we treat s as a constant and assume that our differential equation is dy/dx = Gy, the general solution of which, after n steps, is y_n = Ae^ns, where A is an arbitrary constant, s = hG, and x = nh.

When Simpson's rule is used, the corresponding difference equation is (1 - s/3) y_n+1 - (4s/3) y_n - (1 + s/3) y_n-1 = 0, (1) with the general solution y_n = Ar₁ⁿ + Br₂ⁿ, (2) in which r₁ and r₂ are the roots of the quadratic equation (1 - s/3) r² - (4s/3) r - (1 + s/3) = 0. From this equation we obtain the derivative dr/ds = (r² + 4r + 1)²(12r² + 12r + 12)^-1, which is never negative. Hence the roots r₁ = [2s/3 + (1 + s²/3)^1/2] (1 - s/3)^-1, r₂ = [2s/3 - (1 + s²/3)^1/2] (1 - s/3)^-1, are monotone increasing functions for all real values of s, except for a discontinuity in r₁ at s = 3.

Moreover, the roots r₁ and r₂ are analytic within a circle of radius 3^1/2 with center at the origin in the complex s plane. Through terms of degree five in s the power series for r₁ and r₂ are respectively r₁ = 1 + s + s²/2! + s³/3! + s⁴/4! + s⁵/72 + ··· = e^s + s⁵/180 + ··· , (3) and r₂ = -1 + s/3 - s²/18 - s³/54 + 5s⁴/648 + 5s⁵/1944 + ··· (4) Obviously r₁ is the desired root and r₂ is the unwanted root that produces the oscillation.

Quite apart from questions of stability the process of numerical integration with Simpson's rule requires that the quantity s/3 must be numerically less than unity and in practical computation should be considerably less than unity. Cf. Milne [5], p. 67. We shall therefore assume that |s| \lt 1. Table I shows to six decimal places the value of r₁ and r₂ for s ranging from -1 to +1 at steps of 0.1. It is evident that in this range r₂ is numerically less than one if s is positive, greater than one if s is negative. Thus the oscillating term increases exponentially with n if G is negative, decreases if G is positive.

Now suppose that after k steps of the process we recompute y_k from the values already found, using Newton's “three eighths rule”, to obtain y_k^* = y_k-3 + (3s/8)(y_k + 3y_k-1 + 3y_k-2 + y_k-3). Then we replace the originally computed y_k by the arithmetic mean y_k = (y_k + y_k^*)/2. (5) From (2) and (5) we find that y_k = Ar^k-3₁K(r₁) + Br^k-3₂K(r₂), (6) in which K(r) is defined by the equation K(r) = [r³ + 1 + (3s/8) (r + 1)³]/2. (7) This function K(r) is the key to the problem.

For by means of the series for r₁ and r₂ it can be shown that K(r₁) = r₁³ + s⁵/96 + ··· (8) while K(r₂) = s/2 - s²/4 + ··· . Hence equation (3) becomes y_k = Ar₁^k + terms of degree 5 and higher + Br^k-3₂(s/2) + terms of degree 2 and higher. (9) Comparing y_k with y_k we note that the desired solution is substantially unchanged, and agrees with the true solution e^ks through terms of degree 4 in s, while the unwanted solution has been decreased roughly by a factor of magnitude s/2.

Table I shows to six decimal places the values of K(r₁) and K(r₂) in the range from s = -1 to s = +1. It is seen that in this interval the absolute value of K(r₂) is always less than unity.

Consider now the propagation of a single error starting at n = 0 and modified after every group of k steps by means of formula (5). Since E_n is a solution of equation (1), in the mth group of k steps the error E_n can be expressed by formula (2) in the form E_n = a_mr₁^j + b_mr₂^j, (10) provided n = mk + j and j < k. But if j = k the value of E_n can be expressed approximately as E_n = a_mr₁^k + b_mr^k-3₂K(r₂). (11) To obtain this result one must replace K(r₁) by its approximate value r₁³, as shown in equation (8).

Similarly in the (m + 1)th group of steps we may let E_n = a_m+1r₁^j + b_m+1r₂^j, where n = (m + 1)k + j and j < k. The coefficients a_m+1 and b_m+1 are connected to the coefficients a_m and b_m by the following equations, written in matrix form: (r^-1₁ r^-1₂ 1 1)(a_m+1 b_m+1) = (r^k-1₁ r^k-1₂ r₁^k r^k-3₂K(r₂))(a_m b_m). One may verify the above statement by noting that both members of the first equation are equal to E_mk+k-1 and both members of the second equation are equal to E_mk+k.

Left-hand multiplication by the inverse of the matrix on the left leads to (a_m+1 b_m+1) = M (a_m b_m) (12) in which M = (u v 0 w), where u = r₁^k, v = r₂^kP, w = r₂^kQ, and P = - [r₁ - r₁K(r₂)r^-3₂] (r₁ - r₂)-1, Q = [r₁ - r₂K(r₂)r³₂] (r₁ - r₂)^-1.

The quantities a_m and b_m, and consequently (by equation (10)) the quantity E_n also, will approach zero as m becomes infinite provided both latent roots of M, namely u and w, are less than one in absolute value. For the case under consideration where s lies between -1 and 0 the value of u is always less than one, as we see from table I.

It remains to examine the other latent root w = r₂^kQ. The quantity Q is a function of s alone and its values for s in the interval form -1 to 0 are shown in table II. If we define q by the equation q = - log Q/log(-r₂) (13) it can be shown that for s between -1 and 0 and for k an integer less than q the latent root w will be less than one in absolute value. Hence to assure that the propagated error E_n will approach zero it is sufficient to choose a value of k which is less than q. For convenience of computers some values of q for s between -1 and 0 are supplied in table II.

It may be noted that the foregoing analysis does not strictly apply if k is less than 3 since formula (6) on which the reasoning depends was derived with the tacit assumption that k is not less than 3. Nevertheless machine tests indicate that the convergence for k less than 3 is just about what might be expected on the basis of the above analysis. However, it is unwise for other reasons than stability to use values of s numerically greater than 0.8, so that the accuracy of table II in this range is unimportant.

To illustrate the foregoing theory several computations were performed on the Alwac III-E at Oregon State College for the system dy/dx = -y, y(0) = 1. In this case s < 0, |r₂| > 1, and Simpson's rule, if uncorrected, produces instability.

Table III shows the difference E = e^ns - y_n between the true solution e^ns and the computed solution y_n after n steps of the computation. Six values of k are used in table III, k = ∞ (that is, no stabilization), k = 169, 39, 19, 5, and 3. Four values of s are used, namely s = -0.10, -0.07, -0.04, and -0.01. The number n of steps in the computations varies from 300 for larger values of -s to 2000 for the smallest. Not all computations were carried to the full number of steps shown at the left, hence some columns are partially blank.

The number of decimal places is indicated for each division of the table. For example at s = -0.10, k = 169, and n = 300, the entry -31 means -0.000031, while for s = -0.04, k = 169, n = 500, the entry -6 means -0.00000006.

From table II we obtain the integral parts of q corresponding to the given values of s and find that according to theory the solution should be stable for s = -0.10 if k < 21, for s = -0.07 if k < 30, for s = -0.04 if k < 52, and for s = -0.01 if k < 208. Hence in table III the three right-hand columns should be stable for all four values of s, four right-hand columns should be stable for s = -0.04, and five right-hand columns for s = -0.01. The computations appear to conform to the theory, since the error is negligible in these cases. Occasional errors of one unit in the last place are to be explained by the accidents of roundoff, for if they were due to instability they would increase with n. (The error -2 for k = 169 at n = 2000 is unexplained, but apparently is not due to instability, as it persisted without increasing through many steps between 1800 and 2000.)

The results shown in table III illustrate the normal situation occurring in practice, where if h is properly chosen and the machine operates correctly, the only source of error is roundoff.

We note that of two consecutive values of k the odd value is likely to give somewhat better results. For if k is even no stabilization occurs at odd values of n, and since Simpson's rule operates over two intervals the effect of stabilization only reaches the odd entries indirectly through the derivative y′.

It is the intent of the authors to treat differential systems of higher order in a future paper.