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 1 + 2 k 2 + 2 k 3 + k 4 ) + O ( h 5 ) k 1 = ƒ( x i , y i ) k 2 = ƒ( x i + h /2, y i + k 1 h /2) k 3 = ƒ( x i + h /2, y i + k 2 h /2) k 4 = ƒ( x i + h, y i + k 3 h ) (3) When the O ( h 5 ) 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 5 , (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 5 ) 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 ε 1 in y i at the i th step, which might be due either to previous truncation or round-off error. Then the calculated value y * i +1 (as opposed to the true value if ε i were zero) at the next step is given, after several applications of the mean value theorem, by y * i +1 = y i + ε i + h /6{ k 1 + ∂ƒ/∂ y ε i + 2 k 2 + 2[∂ƒ/∂ y + (∂ƒ/∂ y ) 2 h /2]ε i + 2 k 3 + 2[∂ƒ/∂ y (1 + ∂ƒ/∂ y h /2(1 + ∂ƒ/∂ y h /2))]ε i + k 4 + [∂ƒ/∂ y (1 + ∂ƒ/∂ y h (1 + ∂ƒ/∂ y h /2(1 + ∂ƒ/∂ y h /2)))]ε i } (5) where the various partial derivatives are evaluated within a rectangle given by x i ≦ x ≦ x i + h ; | y i - y | ≦ Qh + | ε i | where Q is given below. Consider the true solution y i t to the difference equation. During its evaluation, in order that all values of ƒ( x, y ) used in its calculation lie within a region D , the true solution should not approach the y -boundaries of D closer than | y i t ; - y | > Qh , where if Q is chosen as Q ≧ max x,y ε D | ƒ( x, y ) | ≧ max (| k 1 /2 |, | k 2 /2 |, | k 3 |) then certainly only values within D will be involved in the calculations. The same argument holds for the difference equation containing error ε i at any step. If such a solution does not approach within closer than Qh + | ε i | 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 y i +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 ε i present): η i +1 = y * i +1 - y i +1 = ε i [1 + h /6(∂ƒ/∂ y + 2∂ƒ/∂ y + 2∂ƒ/∂ y + ∂ƒ/∂ y ) + h 2 /6(∂ƒ/∂ y ∂ƒ/∂ y + ∂ƒ/∂ y ∂ƒ/∂ y + ∂ƒ/∂ y ∂ƒ/∂ y ) + h 3 /12(∂ƒ/∂ y ∂ƒ/∂ y ∂ƒ/∂ y + ∂ƒ/∂ y ∂ƒ/∂ y ∂ƒ/∂ y ) + h 4 /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, - M 2 < ∂ƒ/∂ y < - M 1 < 0, where M 2 > M 1 < 0, then for a maximum error ( truncation, or round-off, or both ) E in absolute value at each step the Kutta fourth-order numerical integration procedure has total error at the i-th step, i arbitrary, in the region D * , of | ε i | ≦ 2 E / hM 1 (7) where the step-size h is taken to be h < min ( M 1 / M 2 2 , 4 M 1 3 / M 2 4 . (Obviously if D extends to infinity, the boundaries of D and D * will coincide at infinity. Better (but more complex) bounds could certainly be obtained on Q and therefore D * .) PROOF: For h as given, the multiplier of | ε i | within the absolute value signs on the right-hand side of (8) below is positive and the inequality does indeed hold, by applying the bounds on ƒ v of the theorem to (6) above.In fact, for h as given η i +1 | = | y * i +1 - y i +1 | ≦ | ε i | | 1 - hM 1 + ( hM 2 ) 2 /2! - ( hM 1 ) 3 /3! + ( hM 2 ) 4 /4!| ≦ | ε i | |1 - hM 1 /2|. (8) The total error at the ( i + 1)th step is, in absolute value, | ε i +1 | = |ρ i +1 | + |η i +1 | ≦ E + | ε i | (1 - hM 1 /2) (9) where ρ i +1 is the sum of the round-off and truncation error in that step and η i +1 the propagated error, 1 - hM 1 /2 < 1, and E > | ρ i | for all i . For i = 0, ε 0 = 0, and therefore certainly | ε 0 | ≦ 2 E / hM 1 . (10) If for a given i , | ε i | ≦ 2 E / hM 1 , then from (9) | ε i +1 | ≦ E + 2 E / hM 1 (1 - hM 1 /2) ≦ 2 E / hM 1 , (11) and the proof holds by induction 1 on i . Note that, if ρ i consists only of truncation error, | ε i +1 | ≦ 2 Ch 5 / hM 1 , (12) where C > C i for all i (if sufficient bounded derivatives on ƒ( x, y ) are assumed), and lim h →0 | ε i +1 | = 0. (13) On the other hand, if ρ i consists of round-off error ρ * also, bounded by | ρ * | for all i , then lim h →0 | ε i +1 | ≦ 2(| ρ * | + Ch 5 )/ hM 1 = ∞, (14) provided | ρ * | 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 | η i +1 | = | y * i +1 - y i +1 | ≧ | ε i | (1 - hM 2 ) ≧ K | ε i |, (15) where now K = (1 - hM 2 ). If, for example, ρ k > R > 0, for all k , one obtains | ε i +1 | = ρ i +1 + | η i +1 | ≧ ρ i +1 + K (ρ i + K (ρ i -1 + ··· + K i +1 ρ 0 ) ···) (16) ≧ R (1 + K + K 2 + ··· + K i +1 ) ≧ R ( i + 2) K i +1 ∴lim h →0 | ε i +1 | ≧ ( i + 2) R (17) Therefore, in the case ρ k > 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 h k , 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 ≦ |∂ƒ/∂ y | < M . (18) Then the propagated error in D * at the ( i + 1)th step, due to a given (round-off or truncation error) ε i at the i th step, is | y * i +1 - y i +1 | ≦ | ε i | (1 + hM + ( hM ) 2 )/2 + ( hM ) 3 )/3! + ( hM ) 4 )/4!) (19) or | η i +1 | ≦ | ε i | e hM , (20) where η i +1 is the propagated error and E and h are given as before; and one obtains an error bound similar to most classical ones | ε i +1 | ≦ E (1 + K + K 2 + ··· + K i ) ≦ E K i +1 - 1/ K - 1 = E e ( i +1) h M - 1/ e h M - 1, (21) where K = e h M . In both the stable and unstable case, ignoring round-off, the process is obviously convergent in D as h → 0, since D * approaches D as h → 0 and | ε i +1 | → 0. For the case ƒ y < 0 equations (7) and (14) above give an algorithm for a bound on the stepsize h ; a similar analysis could be based on equation (21). Suppose the round-off error per single-precision step is bounded by | ρ *(1) ( h ) |, per n -precision step by | ρ *(n) ( h ) |. (Here n could just as well be number of digits as number of words.) Then, since | ε i +1 | ≦ 2(|ρ *( n ) ( h )| + Ch 5 )/ hM 1 , (22) for | ε i +1 | required to be less than B > 0, choose h < M 1 / M 2 2 such that 2 Ch 5 / hM 1 < B /2 (23) or h 4 < M 1 B /4 C . (23a) Then for that value of h , choose the value of n such that | ρ *( n ) ( h ) | < BhM 1 /4. (24) This will guarantee | ε i +1 | < B . Obviously, such an h is not necessarily the best possible. The computation of | ρ *( n ) ( h ) |, while intricate, could very well be done by the computer itself, if procedures for forming formal derivatives of machine functions were included in whatever compiling techniques were used with the routine. If the calculated solution approached the boundaries of D * , the step-size h would have to be decreased, and the values of h accordingly readjusted in order to continue to apply the theorem. If h and n are picked in this fashion, the numerical machine solution can be made as close to the true solution as desired, and the process will converge. As for “stability” of the solution process under round-off, as would appear intuitive, the propagated error can be guaranteed to go to zero with the step-size h only if the round-off error per step varies with h as O ( h k ) with k > 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.