Stability of a Numerical Solution of Differential Equations

Publisher Website
Google Scholar
Abstract
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 0 at x = x 0 is propagated through subsequent steps according to the second order difference equation (1 - s n +1 /3) E n +1 - (4 s 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 - (4 s /3) y n - (1 + s /3) y n -1 = 0, (1) with the general solution y n = Ar 1 n + Br 2 n , (2) in which r 1 and r 2 are the roots of the quadratic equation (1 - s /3) r 2 - (4 s /3) r - (1 + s /3) = 0. From this equation we obtain the derivative dr / ds = ( r 2 + 4 r + 1) 2 (12 r 2 + 12 r + 12) -1 , which is never negative. Hence the roots r 1 = [2 s /3 + (1 + s 2 /3) 1/2 ] (1 - s /3) -1 , r 2 = [2 s /3 - (1 + s 2 /3) 1/2 ] (1 - s /3) -1 , are monotone increasing functions for all real values of s , except for a discontinuity in r 1 at s = 3. Moreover, the roots r 1 and r 2 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 1 and r 2 are respectively r 1 = 1 + s + s 2 /2! + s 3 /3! + s 4 /4! + s 5 /72 + ··· = e s + s 5 /180 + ··· , (3) and r 2 = -1 + s /3 - s 2 /18 - s 3 /54 + 5 s 4 /648 + 5 s 5 /1944 + ··· (4) Obviously r 1 is the desired root and r 2 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 1 and r 2 for s ranging from -1 to +1 at steps of 0.1. It is evident that in this range r 2 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 + (3 s /8)( y k + 3...