### Stability of a Numerical Solution of Differential Equations—Part IIStability of a Numerical Solution of Differential Equations—Part II

Access Restriction
Subscribed

 Author Milne, W. E. ♦ Reynolds, R. R. Source ACM Digital Library Content type Text Publisher Association for Computing Machinery (ACM) File Format PDF Copyright Year ©1960 Language English
 Subject Domain (in DDC) Computer science, information & general works ♦ Data processing & computer science Abstract 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 $\textit{m}$ or a system of equations of total order $\textit{m}$ can be expressed by $\textit{m}$ simultaneous equations of the first order $\textit{dyi}/\textit{dt}$ = $\textit{Fi}(\textit{y}1,$ …, $\textit{ym},$ $\textit{t}),$ $(\textit{i}$ = 1, 2, …, $\textit{m})$ (1) in $\textit{m}$ dependent variables $\textit{yi}.$ Let y denote a vector in $\textit{m}-space$ which has the components $\textit{y}1,$ &hellip, $\textit{ym}.$ Let z, with components $\textit{z}1,$ $\textit{z}2,$ …, $\textit{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 $\textit{t}.$ The differential system (1) can be written in vector form as $\textit{d}y/\textit{dt}$ = F(y, $\textit{t})$ (2) and the varied vector y + z satisfies $\textit{d}y/\textit{dt}$ + $\textit{d}z/\textit{dt}$ = F (y + z, $\textit{t}).$ (3) Let us assume that the functions $\textit{Fi}$ have continuous partial derivatives $∂\textit{Fi}/∂\textit{yi}$ = $\textit{aij}$ and let G denote the $\textit{m}$ × $\textit{m}$ matrix $[\textit{aij}].$ From equations (2) and (3) it may be shown that to terms of the first order in the small quantity z we have $\textit{d}z/\textit{dt}$ = Gz. (4)The matrix G is usually a variable depending on the $\textit{y}'s$ and $\textit{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 $\textit{d}x/\textit{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 λ1, λ2, … $λ\textit{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 $\textit{m}$ independent equations $\textit{dxi}/\textit{dt}$ = $λ\textit{ixi},$ $(\textit{i}$ = 1, 2, …, $\textit{m}),$ (6) where $\textit{xi}$ is the $\textit{i}th$ component of x. The solutions are $\textit{xi}$ = $\textit{cie}λ\textit{it},$ in which the $\textit{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 λ1 = λ2 = λ3, while the remaining roots are distinct, we would have $\textit{dx}1/\textit{dt}$ = $λ1\textit{x}1$ + $\textit{x}2,$ $\textit{dx}2/\textit{dt}$ = $λ1\textit{x}2$ + $\textit{x}3,$ (7) $\textit{dx}3/\textit{dt}$ = $λ1\textit{x}3.$ The solutions of the system (7) are $\textit{x}3$ = $\textit{c}3\textit{e}λ1\textit{t},$ $\textit{x}2$ = $(\textit{c}2$ + $\textit{c}3\textit{t})\textit{e}λ1\textit{t},$ $\textit{x}1$ = $(\textit{c}1$ + $\textit{c}2\textit{t}$ + $\textit{c}3\textit{t}2/2)\textit{e}λ1\textit{t}.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 $\textit{S}$ defined by the equation $\textit{S}ƒ(\textit{t})$ = $ƒ(\textit{tn}+1)$ - $ƒ(\textit{tn}-1)$ - $(\textit{h}/3)[ƒ′$ $(\textit{tn}+1)$ + 4ƒ′ $(\textit{tn})$ + ƒ′ $(\textit{tn}-1)],$ (8) where ƒ′ $(\textit{t})$ means $\textit{df}/\textit{dt}$ and where $\textit{tn}$ = $\textit{nh}.$ If $ƒ(\textit{tn}+1)$ has been computed by means of Simpson's rule from $ƒ(\textit{tn}-1)$ and the values of ƒ′ $(\textit{t}),$ it is clear that $\textit{S}ƒ(\textit{t})$ = 0. Since in Milne's method the final values of the variables $\textit{yi},$ are found by Simpson's rule we may assume1 that $\textit{Syi}$ = 0 for $\textit{i}$ = 1 to $\textit{m}.$ Then $\textit{S}y$ = 0, and since z represents an error inherited from previous steps we have $\textit{S}(y$ + z) = 0, whence $\textit{S}z$ = 0 also. Moreover, since x = Tz, it follows that $\textit{S}x$ = $\textit{S}Tz$ = $T\textit{S}z$ = 0. Hence $\textit{Sxi}$ = 0, $\textit{i}$ = 1, …, $\textit{m}.$ (9) From equations (6), (8), and (9) it is clear that in case $λ\textit{i}$ is not a multiple root the values $\textit{xi}$ satisfy the difference equation (1 - $\textit{s}/3)\textit{xi}(\textit{tn}+1)$ - $(4\textit{s}/3)$ $\textit{xi}$ $(\textit{tn})$ - (1 + $\textit{s}/3)$ $\textit{xi}$ $(\textit{tn}-1)$ = 0, in which $\textit{s}$ = $\textit{h}λ\textit{i}.$ This is identical with equation (1) of Part I. Hence we see by equation (2) of Part I that $\textit{xi}$ is expressed in the form $\textit{xi}$ $(\textit{tn})$ = $\textit{Ar}1\textit{n}$ + $\textit{Br}2\textit{n},$ where $\textit{r}1$ and $\textit{r}2$ are roots of the quadratic equation (1 - $\textit{s}/3)$ $\textit{r}2$ - $(4\textit{s}/3)$ $\textit{r}$ - (1 + $\textit{s}/3)$ = 0. (10)We now consider the following three cases: $\textit{Case}$ 1, where $λ\textit{i}$ is real and simple. In this case the treatment of Part I applies without change to the stabilization of the particular component $\textit{xi}.\textit{Case}$ 2, where the root $λ\textit{i}$ is complex but not a multiple root. Consider equation (10) as the defining equation for a complex function $\textit{r}$ of a complex variable $\textit{s}.$ This two-valued function has branch points at $\textit{s}$ = ± $\textit{i}√3.$ We make branch cuts in the two-sheeted Riemann surface in the $\textit{s}-plane$ along the axis of imaginaries from $\textit{s}$ = + $\textit{i}√3$ up to infinity, and from $\textit{s}$ = $-\textit{i}√3$ down to infinity, as shown in figure 1 (here $\textit{i}$ is the unit of imaginaries).These branch cuts are mapped in the $\textit{r}-plane$ on the circle $\textit{BCDC}′$ with center at (-2, 0) and radius √3. Corresponding points in the $\textit{s}-$ and $\textit{r}-planes$ are denoted by the same letters. The function which we have called $\textit{r}2$ is mapped conformally on the $\textit{interior}$ of the circle $\textit{BCDC}′,$ the function $\textit{r}1$ on the $\textit{exterior}.$ The interior of the unit circle in the $\textit{s}-plane$ (which is our only concern in this discussion) is mapped on the shaded regions in figure 2. The segment $\textit{BD}$ of the imaginary axis in the $\textit{s}-plane$ goes into the unit circle $\textit{BADA}′$ with center at the origin in the $\textit{r}-plane.$ When $λ\textit{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 $\textit{s}$ = $\textit{h}λ\textit{i}$ depends on the latent roots $\textit{u}$ and $\textit{w}$ of the matrix $\textit{M}.$ These latent roots are approximately expressed by the formulas $\textit{u}$ = $\textit{r}1\textit{k},$ $\textit{w}$ = $\textit{r}2\textit{k}\textit{Q},$ just as in Part I.(We recall that the above simple forms for $\textit{u}$ and $\textit{w}$ were obtained by replacing $\textit{K}(\textit{r}1)$ by $\textit{r}13.$ The resulting error in $\textit{u}$ and $\textit{w}$ can be shown to be of the order of $\textit{s}5.$ In the domain of $\textit{s}$ that is used in practice we may ignore these errors.) When $λ\textit{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 $\textit{s}$ = $\textit{h}λ\textit{i}$ depends on the latent roots $\textit{u}$ and $\textit{w}$ of the matrix $\textit{M}.$ These latent roots are approximately expressed by the formulas $\textit{u}$ = $\textit{r}1\textit{k},$ $\textit{w}$ = $\textit{r}2\textit{k}\textit{Q},$ just as in Part I.(We recall that the above simple forms for $\textit{u}$ and $\textit{w}$ were obtained by replacing $\textit{K}(\textit{r}1)$ by $\textit{r}13.$ The resulting error in $\textit{u}$ and $\textit{w}$ can be shown to be of the order of $\textit{s}5.$ In the domain of $\textit{s}$ that is used in practice we may ignore these errors.) ISSN 00045411 Age Range 18 to 22 years ♦ above 22 year Educational Use Research Education Level UG and PG Learning Resource Type Article Publisher Date 1960-01-01 Publisher Place New York e-ISSN 1557735X Journal Journal of the ACM (JACM) Volume Number 7 Issue Number 1 Page Count 11 Starting Page 46 Ending Page 56

#### Open content in new tab

Source: ACM Digital Library