Error Bounds for the Runge-Kutta Single-Step Integration Process

Abstract

Although numerous writers have stated that the class of single-step (“Runge-Kutta”—like) methods of numerical integration of ordinary differential equations are stable under calculational or round-off error, no one has given formal equations for the bounds on the propagated error to indicate this stability. Rutishauser [1] justifies the stability by noting that there is only one solution to the approximating difference equation, and Hildebrand [2] calculates a propagated error bound for the simplest (Euler) case. However, the latter bound does not indicate the stability for even that case. It is the purpose of this paper to investigate this “stability” of the Kutta fourth order procedure for integration of the ordinary differential equation dy / dx = ƒ( x, y ), (1) where ƒ( x, y ) possesses a continuous first-order partial derivative with respect to y throughout a region D in which the integration is to take place. (By alteration of the proof below, this condition can be replaced by a Lipschitz condition.) Since the Kutta process is the most complicated of such single-step procedures, it should be apparent that similar error bounds can be derived for the various other single-step methods of same or lower order (and in particular the Gill variant method, probably most often used in machine integration because of the storage savings). It is plausible that such error bounds can also be extended to the stable (extrapolation) multi-step methods, such as the Adams method, and to systems of ordinary differential equations. If the variational equation d η/ dx = ∂ƒ( x, y )/∂ y η (2) for the above ordinary differential equation has ƒ _v ( x, y ) < 0 throughout a region D the ordinary differential equation is termed stable in D , and for small enough variations in the initial conditions the absolute value of the propagated error decreases with growing x . Todd [3], Milne [4], Rutishauser [1], and others have shown that numerous multi-step numerical integration techniques are unstable in that even when the differential equation is stable, the difference equation will introduce spurious solutions for any step-size h . For the Kutta fourth order process, as seen below, this is not the case; for the stable differential equation the propagated error in the difference approximation remains bounded for small enough (but not too small!) step-size h ; and for a given value of x , bounds on the propagated error decrease to a minimum given as a function of the round-off error as the step-size is decreased. Similar statements can be proved (but are not proved here) for other single-step processes. For no round-off (an “infinite word-length machine”) the process converges as h goes to zero. In addition, an algorithm for determining the step-size as a function of the partial derivative is given below so as to keep the propagated error within a given bound. The classical Kutta procedure [5] gives the value of y at the ( i + 1)th step in terms of the value at the i th step and the step-size h as follows: y _{i +1} = y _i + h /6( k ₁ + 2 k ₂ + 2 k ₃ + k ₄ ) + O ( h ⁵ ) k ₁ = ƒ( x _i , y _i ) k ₂ = ƒ( x _i + h /2, y _i + k ₁ h /2) k ₃ = ƒ( x _i + h /2, y _i + k ₂ h /2) k ₄ = ƒ( x _i + h, y _i + k ₃ h ) (3) When the O ( h ⁵ ) term is neglected, the value of y _{i +1} , here called y ^t _{i +1} , is only an approximation. An error bound for the truncation error at each step is given by Bieberbach [6, 7] and Lotkin [8]. This guarantees that for h small enough the truncation error at a single step, starting at y _i , will be bounded by | y _{i +1} - y ^t _{i +1} | ≦ C _{i +1} h ⁵ , (4) where C _{i +1} ≧ 0 is a function only of i , the function ƒ( x, y ), and its partial derivatives of the first three orders; and y ^t _i is the true solution to the equations (3) truncated so that the O ( h ⁵ ) term does not appear. If the function and its derivaties of the first three orders exist and are bounded throughout a region, then all C _i would be bounded in the region. Suppose there exists an error ε ₁ in y _...