Stability of a Numerical Solution of Differential Equations—Part II

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

In Part I of this paper [1] the authors have shown that instability in Milne's method of solving differential equations numerically [2] can be avoided by the occasional use of Newton's “three eights” quadrature formula. Part I dealt with a single differential equation of first order. In Part II the analysis is extended to equations and systems of equations of higher order.A differential equation of order m or a system of equations of total order m can be expressed by m simultaneous equations of the first order dyi/dt = Fi(y1, …, ym, t), (i = 1, 2, …, m) (1) in m dependent variables yi. Let y denote a vector in m-space which has the components y1, &hellip, ym. Let z, with components z1, z2, …, zm, denote a small variation in the vector y such as might have been produced by a small error occurring at an earlier value of t. The differential system (1) can be written in vector form as dy/dt = F(y, t) (2) and the varied vector y + z satisfies dy/dt + dz/dt = F (y + z, t). (3) Let us assume that the functions Fi have continuous partial derivatives ∂Fi/∂yi = aij and let G denote the m × m matrix [aij]. From equations (2) and (3) it may be shown that to terms of the first order in the small quantity z we have dz/dt = Gz. (4)The matrix G is usually a variable depending on the y's and t, but it is possible to forecast the general behavior of the error by treating the simple case where G is constant. (See [3] of Part I.) Assuming that G is constant we first of all briefly review some well-known facts about the solution of the differential system (4).If we introduce a new vector variable x in place of the variable z by means of a nonsingular linear transformation x = Tz, it is seen that the system (4) is transformed into dx/dt = (TGT-1)x. (5) In particular the transformation T may be chosen so as to put the matrix TGT-1 into the classical canonical form. (Cf. e.g., Turnbull and Aitkin [3].) Then if the latent roots &lgr;1, &lgr;2, … &lgr;m, of the matrix G are all distinct, the new matrix TGT-1 has these latent roots in the principal diagonal and zeros elsewhere. In these case the differential system (5) separates into m independent equations dxi/dt = &lgr;ixi, (i = 1, 2, …, m), (6) where xi is the ith component of x. The solutions are xi = cie&lgr;it, in which the ci are arbitrary constants. If the latent roots are not all distinct there may occur ones instead of zeros in the diagonal just above the principal diagonal. Wherever a one occurs the latent root to its left equals the latent root below it. In such a case we have in addition to equations of type (6) certain nonhomogeneous linear equations. For example, if &lgr;1 = &lgr;2 = &lgr;3, while the remaining roots are distinct, we would have dx1/dt = &lgr;1x1 + x2, dx2/dt = &lgr;1x2 + x3, (7) dx3/dt = &lgr;1x3. The solutions of the system (7) are x3 = c3e&lgr;1t, x2 = (c2 + c3t)e&lgr;1t, x1 = (c1 + c2t + c3t2/2)e&lgr;1t.It should be noted that multiple roots do not always lead to equations of type (7). Such a case as illustrated by example 3 at the end of the paper. Turning now to numerical integration and the problem of stability we introduce a linear operator S defined by the equation Sƒ(t) = ƒ(tn+1) - ƒ(tn-1) - (h/3)[ƒ′ (tn+1) + 4ƒ′ (tn) + ƒ′ (tn-1)], (8) where ƒ′ (t) means df/dt and where tn = nh. If ƒ(tn+1) has been computed by means of Simpson's rule from ƒ(tn-1) and the values of ƒ′ (t), it is clear that Sƒ(t) = 0. Since in Milne's method the final values of the variables yi, are found by Simpson's rule we may assume1 that Syi = 0 for i = 1 to m. Then Sy = 0, and since z represents an error inherited from previous steps we have S(y + z) = 0, whence Sz = 0 also. Moreover, since x = Tz, it follows that Sx = STz = TSz = 0. Hence Sxi = 0, i = 1, …, m. (9) From equations (6), (8), and (9) it is clear that in case &lgr;i is not a multiple root the values xi satisfy the difference equation (1 - s/3)xi(tn+1) - (4s/3) xi (tn) - (1 + s/3) xi (tn-1) = 0, in which s = h&lgr;i. This is identical with equation (1) of Part I. Hence we see by equation (2) of Part I that xi is expressed in the form xi (tn) = Ar1n + Br2n, where r1 and r2 are roots of the quadratic equation (1 - s/3) r2 - (4s/3) r - (1 + s/3) = 0. (10)We now consider the following three cases: Case 1, where &lgr;i is real and simple. In this case the treatment of Part I applies without change to the stabilization of the particular component xi.Case 2, where the root &lgr;i is complex but not a multiple root. Consider equation (10) as the defining equation for a complex function r of a complex variable s. This two-valued function has branch points at s = ± i√3. We make branch cuts in the two-sheeted Riemann surface in the s-plane along the axis of imaginaries from s = + i√3 up to infinity, and from s = -i√3 down to infinity, as shown in figure 1 (here i is the unit of imaginaries).These branch cuts are mapped in the r-plane on the circle BCDC′ with center at (-2, 0) and radius √3. Corresponding points in the s- and r-planes are denoted by the same letters. The function which we have called r2 is mapped conformally on the interior of the circle BCDC′, the function r1 on the exterior. The interior of the unit circle in the s-plane (which is our only concern in this discussion) is mapped on the shaded regions in figure 2. The segment BD of the imaginary axis in the s-plane goes into the unit circle BADA′ with center at the origin in the r-plane. When &lgr;i is a complex latent root, but not a multiple root, we proceed exactly as in Part I, and find that the stability of the solution corresponding to s = h&lgr;i depends on the latent roots u and w of the matrix M. These latent roots are approximately expressed by the formulas u = r1k, w = r2kQ, just as in Part I.(We recall that the above simple forms for u and w were obtained by replacing K(r1) by r13. The resulting error in u and w can be shown to be of the order of s5. In the domain of s that is used in practice we may ignore these errors.) When &lgr;i is a complex latent root, but not a multiple root, we proceed exactly as in Part I, and find that the stability of the solution corresponding to s = h&lgr;i depends on the latent roots u and w of the matrix M. These latent roots are approximately expressed by the formulas u = r1k, w = r2kQ, just as in Part I.(We recall that the above simple forms for u and w were obtained by replacing K(r1) by r13. The resulting error in u and w can be shown to be of the order of s5. In the domain of s that is used in practice we may ignore these errors.)