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

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 de/dx = ∂ƒ(x, y)/∂y e (2) for the above ordinary differential equation has ƒv(x, y) Qh, where if Q is chosen as Q g maxx,yeD | ƒ(x, y) | g max (| k1/2 |, | k2/2 |, | k3 |) then certainly only values within D will be involved in the calculations.The same argument holds for the difference equation containing error ei at any step. If such a solution does not approach within closer than Qh + | ei | to the y-boundary of D (such a region will be called D* below), then the true solution will not approach closer than Qh to the boundary. In the latter case the rectangle within which the partial derivatives are to be evaluated will always be within the given region D.This gives for the propagated error at the (i + 1)th step (if yi+1 is the value at step i + 1 given no error at step i, and y*i+1 is the value at step i + 1 with error ei present): ei+1 = y*i+1 - yi+1 = ei[1 + h/6(∂ƒ/∂y + 2∂ƒ/∂y + 2∂ƒ/∂y + ∂ƒ/∂y) + h2/6(∂ƒ/∂y ∂ƒ/∂y + ∂ƒ/∂y ∂ƒ/∂y + ∂ƒ/∂y ∂ƒ/∂y) + h3/12(∂ƒ/∂y ∂ƒ/∂y ∂ƒ/∂y + ∂ƒ/∂y ∂ƒ/∂y ∂ƒ/∂y) + h4/24(∂ƒ/∂y ∂fnof;/∂y ∂ƒ/∂y ∂ƒ/∂y)] (6) where the partial derivaties are written out to indicate that each is evaluated at what may be a different point in the rectangle. One can now prove the following:THEOREM: If ∂ƒ/∂y is continuous, negative, and bounded from above and below throughout a region D in the(x, y) plane, -M2 M1 | ri | for all i.For i = 0, e0 = 0, and therefore certainly | e0 | l 2E/hM1. (10) If for a given i, | ei | l 2E/hM1, then from (9) | ei+1 | l E + 2E/hM1 (1 - hM1/2) l 2E/hM1, (11) and the proof holds by induction1 on i.Note that, if ri consists only of truncation error, | ei+1 | l 2Ch5/hM1, (12) where C > Ci for all i (if sufficient bounded derivatives on ƒ(x, y) are assumed), and limh→0 | ei+1 | = 0. (13)On the other hand, if ri consists of round-off error r* also, bounded by | r* | for all i, then limh→0 | ei+1 | l 2(| r* | + Ch5)/hM1 = ∞, (14) provided | r* | is bounded from below by a positive non-zero value, which is generally true for fixed word-length computers; and the error bound is of no value in the limit. In that case, for h small enough | ei+1 | = | y*i+1 - yi+1 | g | ei | (1 - hM2) g K | ei |, (15) where now K = (1 - hM2).If, for example, rk > R > 0, for all k, one obtains | ei+1 | = ri+1 + | ei+1 | g ri+1 + K(ri + K(ri-1 + ··· + Ki+1r0) ···) (16) g R(1 + K + K2 + ··· + Ki+1) g R(i + 2)Ki+1 ∴limh→0 | ei+1 | g (i + 2)R (17)Therefore, in the case rk > R > 0, the error can increase without limit when (i + 1)h is held constant. Thus, as perhaps seems obvious from the start, decreasing the step-size will in this case, speaking roughly, improve the propagated error until the point at which the round-off and truncation error have the same order of magnitude. If on the other hand the round-off error can be made a function of the step-size of order hk, where k > 1, then the total error will go to zero as h → 0, or will remain bounded in the case k = 1. Users of variable wordlength digital computing machines should therefore increase their word-length as h is decreased in order to avoid this lower limit due to round-off error. Users of fixed word-length computers must use multi-precision. This result is not surprising; a specific algorithm for doing it is given below.If ƒy is allowed also to be zero or positive, but bounded throughout D, 0 l |∂ƒ/∂y| 0, choose h 1. The Kutta process is stable in the sense that for a stable differential equation, for a given h, the error will be bounded by (7), where E is given without round-off error by (12) and Bieberbach's bound (4). In the latter case, as h → 0, the truncation error goes to zero, which certainly is not true for many multi-step methods.