Stability of a Numerical Solution of Differential Equations

• Authors: W. E. Milne;R. R. Reynolds
• Venue: Journal of the ACM (JACM)
• Year: 1959

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 E0 at x = x0 is propagated through subsequent steps according to the second order difference equation (1 - sn+1/3) En+1 - (4sn/3) En - (1 + sn-1/3) En-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 yn = Aens, 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) yn+1 - (4s/3) yn - (1 + s/3) yn-1 = 0, (1) with the general solution yn = Ar1n + Br2n, (2) in which r1 and r2 are the roots of the quadratic equation (1 - s/3) r2 - (4s/3) r - (1 + s/3) = 0. From this equation we obtain the derivative dr/ds = (r2 + 4r + 1)2(12r2 + 12r + 12)-1, which is never negative. Hence the roots r1 = [2s/3 + (1 + s2/3)1/2] (1 - s/3)-1, r2 = [2s/3 - (1 + s2/3)1/2] (1 - s/3)-1, are monotone increasing functions for all real values of s, except for a discontinuity in r1 at s = 3. Moreover, the roots r1 and r2 are analytic within a circle of radius 31/2 with center at the origin in the complex s plane. Through terms of degree five in s the power series for r1 and r2 are respectively r1 = 1 + s + s2/2! + s3/3! + s4/4! + s5/72 + ··· = es + s5/180 + ··· , (3) and r2 = -1 + s/3 - s2/18 - s3/54 + 5s4/648 + 5s5/1944 + ··· (4) Obviously r1 is the desired root and r2 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 r1 and r2 for s ranging from -1 to +1 at steps of 0.1. It is evident that in this range r2 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 yk from the values already found, using Newton's “three eighths rule”, to obtain yk* = yk-3 + (3s/8)(yk + 3yk-1 + 3yk-2 + yk-3). Then we replace the originally computed yk by the arithmetic mean yk = (yk + yk*)/2. (5) From (2) and (5) we find that yk = Ark-31K(r1) + Brk-32K(r2), (6) in which K(r) is defined by the equation K(r) = [r3 + 1 + (3s/8) (r + 1)3]/2. (7) This function K(r) is the key to the problem. For by means of the series for r1 and r2 it can be shown that K(r1) = r13 + s5/96 + ··· (8) while K(r2) = s/2 - s2/4 + ··· . Hence equation (3) becomes yk = Ar1k + terms of degree 5 and higher + Brk-32(s/2) + terms of degree 2 and higher. (9) Comparing yk with yk we note that the desired solution is substantially unchanged, and agrees with the true solution eks 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(r1) and K(r2) in the range from s = -1 to s = +1. It is seen that in this interval the absolute value of K(r2) 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 En is a solution of equation (1), in the mth group of k steps the error En can be expressed by formula (2) in the form En = amr1j + bmr2j, (10) provided n = mk + j and j k. But if j = k the value of En can be expressed approximately as En = amr1k + bmrk-32K(r2). (11) To obtain this result one must replace K(r1) by its approximate value r13, as shown in equation (8). Similarly in the (m + 1)th group of steps we may let En = am+1r1j + bm+1r2j, where n = (m + 1)k + j and j k. The coefficients am+1 and bm+1 are connected to the coefficients am and bm by the following equations, written in matrix form: (r-11 r-12 1 1)(am+1 bm+1) = (rk-11 rk-12 r1k rk-32K(r2))(am bm). One may verify the above statement by noting that both members of the first equation are equal to Emk+k-1 and both members of the second equation are equal to Emk+k. Left-hand multiplication by the inverse of the matrix on the left leads to (am+1 bm+1) = M (am bm) (12) in which M = (u v 0 w), where u = r1k, v = r2kP, w = r2kQ, and P = - [r1 - r1K(r2)r-32] (r1 - r2)-1, Q = [r1 - r2K(r2)r32] (r1 - r2)-1. The quantities am and bm, and consequently (by equation (10)) the quantity En 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 = r2kQ. 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(-r2) (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 En 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 r2| 1, and Simpson's rule, if uncorrected, produces instability.Table III shows the difference E = ens - yn between the true solution ens and the computed solution yn 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 s = -0.07 if k s = -0.04 if k s = -0.01 if k 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.