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 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 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