Stepwise Global Error Control in an Explicit Runge-Kutta Method Using Local Extrapolation with High-Order Selective Quenching

,


Introduction
Initial-value problems (IVPs) of the form y = f (x, y) x ∈ x 0 , x f y (x 0 ) = y 0 are often solved numerically using a Runge-Kutta method.Usually, some form of local error control is implemented in a step-by-step manner.However, local error control does not guarantee global error control.Global error could still accumulate and become unacceptably large.Global error control, if implemented, requires a reintegration processsolving the problem again using a smaller stepsize -and is not a stepwise process.
In this paper, we introduce an idea for achieving both local and global error control in a stepwise fashion, so that reintegration is not needed.For simplicity, we restrict our considerations to a nonstiff scalar problem (as opposed to a system of differential equations), for which the solution does not vary significantly in magnitude on x 0 , x f (this allows us to consider absolute error control only, rather than relative error control).
We describe relevant concepts in Section 2; we show how global error can grow, despite local error control, in Section 3; in Section 4 we discuss our approach for stepwise global error control; we describe this approach in algorithmic form in Section 5; and in Section 6 we give a numerical example, as a demonstration of our algorithm.

Relevant Concepts, Terminology and Notation
In this section, we discuss concepts relevant to the rest of paper, including appropriate terminology and notation.The reader is referred to Hairer et al (2000), Butcher (2003), Iserles (2009), Kincaid & Cheney (2002), LeVeque (2007), and many references therein, for discussions of Runge-Kutta methods and local error control in such methods.

Runge-Kutta Methods
The most general definition of a Runge-Kutta (RK) method is a pq k q p = 1, 2, ..., m w i+1 = w i + h i m p=1 b p k p ≡ w i + h i F (x i , w i ; h i ) .
(1) Such a method is said to have m stages (the k q ).We note that if a pq = 0 for all p q, then the method is said to be explicit; otherwise, it is known as an implicit RK method.We will restrict our considerations here to explicit methods.The number of stages is related to the order r of the method, and for explicit methods we always have r m.In the second line of (1), we have implicitly defined the function F. The symbol w is used here and throughout to indicate the approximate numerical solution, whereas the symbol y will be used to denote the exact solution.As a refinement to our notation, we will denote a Runge-Kutta method of order r as RKr and, for such a method, we write In a sense, RKr is defined by F r , although it is understood that, for any r, there are, generally speaking, numerous possible choices for F r .The superscripts in (2) are labels, not exponents.The stepsize h i is given by and carries the subscript because it may vary from step to step.

Error Propagation
We define the global error in a numerical solution generated by RKr at x i+1 by Δ r i+1 ≡ w r i+1 − y i+1 , and the local error at x i+1 by Note the use of the exact value y i in the bracketed term in (3).Again, the superscripts are labels.
We have previously shown (Prentice, 2009) that where ξ i ∈ y i , y i + Δ r i .Equation ( 4) provides the relationship between local and global errors in RKr.We will assume that Δ 0 = 0 (i.e. the initial value is known exactly).We see that the global error at any node x i+1 is the sum of a local error term and a term incorporating the global error at the previous node.For RKr, it is known that On the RHS of these expressions, the superscripts are exponents, and h is a parameter representative of the stepsizes h i .

Local Error Control via Local Extrapolation
Consider two RK methods of order r and order v, i.e.RKr and RKv, with r < v. Let w r i+1 denote the approximate solution at x i+1 obtained with RKr, and similarly for w v i+1 .Let the local error at x i+1 in the RKr method be given by ε r i+1 = β r i+1 h r+1 i , and similarly for Once we have estimated the local error, we can perform error control.Assume that we require that the local error at each step must be less than a user-defined tolerance δ.Moreover, assume that, using stepsize h i , we find In other words, the magnitude of the local error ε r i+1 exceeds the desired tolerance.We remedy the situation by determining a new stepsize h * i from and we repeat the RK computation with this new stepsize.This, of course, gives This procedure is then carried out on the next step, and so on.If the estimated error does not exceed the tolerance, then no stepsize adjustment is necessary, and we proceed to the next step.
Often, we introduce a so-called 'safety factor' σ, as in , where σ < 1, so that the new stepsize is slightly smaller than that given by ( 6).This is an attempt to cater for the possibility that β r i+1 may have been underestimated, due to the assumptions made in deriving (5).The choice of the value of σ is subjective, although a representative value is 0.85.
There is a very important point to be noted.Our procedure for determining β r i+1 hinged on the requirement w r i , w v i = y i .However, we only have the exact solution at the initial point x 0 ; at all subsequent nodes, the solution is approximate.How do we meet the requirement w r i , w v i = y i ?Since the higher-order solution w v i is available, we simply use w v i as input to generate both w r i+1 (using RKr), and w v i+1 (using RKv).In other words, we are assuming that w v i is accurate enough, relative to w r i , to be regarded as the exact value -an assumption entirely consistent with the assumption made in deriving ( 5).This means that we determine the higher-order solution at each node, and this solution is used as input for both RK methods in computing solutions at the next node.This form of local error control is known as local extrapolation.

The Problem
We assume that w v i is used to generate w r i+1 and w v i+1 .Such value of w r i+1 (and associated quantities) will be denoted w rv i+1 .Hence, we have Hence, We see that the presence of global error in the higher-order solution does not affect the expression for β r i+1 obtained under the assumption w r i , However, the expression for Δ rv i+1 informs of a potential problem: we have where Δ v i is the global error in w v i .In (7), we see that a subtractive cancellation ensures that the Δ v i term does not enter directly into the estimate for β r i+1 .Nevertheless, even if β r i+1 h r+1 i δ, we could still have Δ rυ i+1 > δ, perhaps substantially so, if Δ v i is large.Moreover, we should certainly expect that Δ v i could become large under iteration (i.e. as i increases), since global error is essentially an accumulation of local errors.The point here is that, even if local error control is effective, the global error Δ rv i+1 could become large, and could grow in an uncontrolled fashion.Usually, the approach to controlling global error involves estimating the global error after the RK solution has been obtained on the entire interval of integration x 0 , x f , and then repeating the RK computation on x 0 , x f , using a suitably reduced stepsize.Such an approach is termed reintegration.
Our ambition is to develop an algorithm through which the global error in w rv i+1 can be estimated and controlled in a stepby-step manner, without the need for reintegration, while at the same time controlling local error in the usual manner of local extrapolation.To achieve this end, we will use a third RK method of very high order, and we will exploit the safety factor mentioned previously.
Our motivation for developing this algorithm is both academic and practical: academically speaking, it is natural to ask whether or not global error can be controlled in a stepwise manner, given that local error is controlled in such fashion; practically speaking, real-time computations which require the use of an accurate solution, such as a feedback loop in a control system, cannot make use of a reintegration process, and so are reliant on the quality of output generated by stepwise algorithms.

The Solution
Say we have three explicit methods available (RKr, RKv and RKz), with r < v z, so that RKz is of much higher order than RKr and RKv.We would suggest v = r + 1 and z = v + 2, at least.Let h * i denote the stepsize for which where δ is a user-imposed tolerance.Of course, since the safety factor σ is less than unity, we have The quantity σh * i is the de facto stepsize used, as and when required, in local extrapolation.We implement local extrapolation in the usual way, using RKr and RKv, propagating w v i at each step.Simultaneously, we implement RKz at the same nodes.At each node we have Additionally, we can propagate w z i in RKr, which gives where w rz i+1 is the solution obtained using RKr with w z i as input.These expressions give since r z.We thus have a reliable estimate of the components of Δ rv i+1 = ε r i+1 + α rv i Δ v i .Now, say a suitable stepsize adjustment has been made, and the solutions w rv i+1 , w rz i+1 , w v i+1 , w z i+1 at x i+1 = x i + σh * i have been computed.We have If this is the case, we proceed to the next step.
If, however, we find we must conclude that Δ v i has become too large.We then set , because RKz is of much higher order than RKv, and we recalculate w rv i+1 and w v i+1 .In other words, w v i is replaced with a much more accurate value.This will yield so that w rv i+1 and w v i+1 will now have relatively small global error ∝ Δ z i accumulated from previous iterations.Effectively, we have greatly reduced, or quenched, the global error present in w rv i+1 and w v i+1 , due to RKv.We only select to perform a quench when we encounter the case in (10); the case (9) does not require a modification of w v i (and, hence, w rv i+1 and w v i+1 ).It may occur, and often does, that a stepsize adjustment via local extrapolation is not required, simply because the stepsize is already small enough.In such case we must, nevertheless, still test the inequality (10) and perform a quench, if necessary.
Usually, in local extrapolation, we would use w rv i+1 and w v i+1 to estimate the local error ε r i+1 but, if w rz i+1 is available, we should rather use (8) to estimate ε r i+1 , because such estimate is more reliable.The importance of the safety factor is apparent: σ ensures that ε r i+1 is strictly less than δ, so that the global error component α rv i Δ v i can be accommodated to some degree.The extent of this accommodation is determined by σ r+1 .Say σ = 0.85 and r = 3, so that σ r+1 = 0.522.Then, assuming both components of Δ rv i+1 have the same sign, α rv i Δ v i can be accommodated up to a magnitude of 0.478δ, before quenching is needed.
Lastly, we emphasize that, at each node x i+1 , it is w rv i+1 that is presented as the solution to the IVP -this is the numerical solution for which both local and global error control has been performed.

The RKrvQz Algorithm
We describe the sequential execution of the algorithm, which we designate RKrvQz, on a generic subinterval [x i , x i+1 ] : 1. Use w v i and w z i in RKr, RKv and RKz to generate w rv i+1 , w rz i+1 , w v i+1 and w z i+1 .Use h i = h i−1 as the stepsize (we discuss the case of h 0 in the Appendix).
2. Estimate ε r i+1 using w rv i+1 − w v i+1 or w rz i+1 − w z i+1 .The former is the usual local extrapolation approach, whereas the latter is more reliable.

If ε r
i+1 > δ, determine a new stepsize h i = σh * i , and repeat steps 1 and 2 using this new stepsize.Then go to step 5.

Estimate
6.If Δ rv i+1 > δ, set w v i = w z i (quenching) and repeat steps 1 and 2 using the stepsize determined in step 3. Then go to step 8.

If Δ rv i+1
δ, go to step 8. 8.We now have the numerical solutions w rv i+1 , w v i+1 , w rz i+1 and w z i+1 , at x i+1 = x i + h i , with the magnitude of both the local and global error in w rv i+1 less than the tolerance δ.

Numerical Example
By way of example, we apply RK34Q8 (r = 3, v = 4, z = 8) to the IVP y = ln 1000 100 y x ∈ [0, 100] The coefficient in the differential equation has been chosen so that y does not exceed 1000 on the interval of integration, i.e. y does not vary substantially, so that absolute error control is suitable.The exact solution is, of course, y (x) = e ln 1000 100 x .
As we shall see, the RK3 global error in this problem is a rapidly increasing function of x, and so it is an ideal problem for demonstrating the capabilities of RKrvQz.The RK methods used in RK34Q8 are RK3 (Kincaid & Cheney, 2002), the 'classical' RK4 (LeVeque, 2007) and Fehlberg's RK8 (Butcher, 2003).The tableaux for these methods are given in Tables 1 and 2.
In Figure 1 we show ε 3 i and α 34 i Δ 3 i for δ = 10 −4 , 10 −8 .The local errors are always less than δ, with the 'zigzag' shape of the curves reflecting stepsize adjustments.Nevertheless, we see that α 34 i Δ 3 i increases monotonically and, at x = 100, it is about 100 times larger than δ.
Global errors obtained with local error control only (the sum of the errors shown in Figure 1), and with RK34Q8, are shown in Figure 2.For the former, the global error increases significantly, even though the local error has been controlled -a potent example of the problem discussed in Section 3. The RK34Q8 global error, however, remains bounded by δ on the entire interval of integration -a vivid demonstration of the capability of RK34Q8.Here, the zigzag nature of the error curves is due to quenching (the propagation of w 8 i ).In all calculations shown in these figures, we used σ = 0.85, and the initial stepsize h 0 was determined using the procedure described in the Appendix.
It may be prudent to keep track of the global error in RK8.For this, we estimate ε 8 i+1 using Richardson extrapolation (Butcher, 2003) We then assume With Δ 8 0 = 0 and Δ 8 1 = ε 8 1 , we can estimate Δ 8 i+1 as RK34Q8 proceeds.If we detect that Δ 8 i+1 is approaching δ in magnitude, then we could increase δ.This would mean that the tolerance on the global error increases occasionally, as RK34Q8 proceeds, but this is better than incorrectly estimating ε r i+1 and α rv i Δ v i .A possible condition for increasing δ would be Δ 8 i+1 ∼ 0.01δ, with δ being increased to 2δ, perhaps.We must state here that this is purely a speculation on our part; for the example considered here, Δ 8 i+1 < 0.01δ on [0, 100] always.Indeed, in Figure 3 we show actual and estimated values of Δ 8 i for δ = 10 −4 , 10 −8 .It is clear that in both cases, Δ 8 i+1 δ.Also, the estimates are good, particularly for δ = 10 −4 .
In Figure 4, we show global errors in RK34Q8 with σ = 0.9.The increases in the safety factor results in more quenches, in comparison with the error curves in Figure 2.

Comments
A few comments, pertaining mostly to possible future research, are appropriate: 1.The use of RKz means that RKrvQz requires greater computational effort than standard local error control via local extrapolation.This, of course, is the price we must pay for controlling global error, in addition to local error.However, in comparison with reintegration, the extra effort may not be all that significant.In reintegration, we would need to use RKr and RKv to obtain solutions on the interval of integration using a smaller stepsize (the RKv solution would be needed to confirm the quality of the RKr solution), after having performed local error control on the entire interval.We suspect that the additional computational effort in using RKr and RKv a second time in a reintegration process would probably not be all that different from the computational effort involved in using RKz in RKrvQz.In Section 3, we gave our motives for developing RKrvQz, and we are sure that any extra effort in using RKz is a small price to pay for achieving simultaneous stepwise local and global error control.We are quite sure that, whenever and wherever possible, accuracy must take precedence over efficiency.This notwithstanding, we discuss some possible improvements in the efficiency of RKrvQz in #2 and #3 below.
2. The safety factor σ can be used to control the magnitude of α rv i Δ v i .Instead of demanding (9), we rather demand which is more stringent.This implies and if σ is close to unity, α rv i Δ v i must necessarily be small, otherwise quenching must be performed.Since α rv i Δ v i is small it is not unreasonable to assume that α v i Δ v i will also be small, so that the estimate of ε r i+1 using w rv i+1 − w v i+1 will be reliable.This simply means that it would not be necessary to determine w rz i+1 at each node, which might improve the efficiency of the algorithm.
3. Efficiency could also be improved by using an embedded RK pair for RKr and RKv.This would reduce the total number of stage evaluations at each step.A well-known example of such a method is Fehlberg's embedded RK(4,5) pair (Fehlberg, 1970).However, we must note that for any given r and v there is no guarantee that an embedded RK(r, v) pair exists.
4. We could use two tolerances δ 1 and δ 2 , imposed on the local and global error, respectively, with δ 1 < δ 2 .In other words, we do not impose the same level of accuracy on the global error as we do on the local error.This might contradict the objectives in #2 above, though.5.In our numerical example, we estimated ε z i+1 by means of Richardson's extrapolation, and then used (4) with F z y (x i , ξ i ; h i ) ≈ f y x i , w z i .A more accurate, and more efficient, estimate might be obtained by using a higherorder method RKz , with z > z.A high-order embedded RK(z, z ) pair might be useful here, such as Fehlberg's RK(7,8), although error control via local extrapolation with this particular pair is not accurate when f is a function of x only (Hairer et al, 2000).Verner has offered an embedded RK(5,6) pair (Verner, 1978), but this would restrict r and v to 2 and 3, respectively, whereas we would probably prefer to have r = 3 or 4.
6.We have considered IVPs for which the solution does not vary considerably in magnitude on the interval of integration.This has allowed us to consider absolute error control only (wherein a uniform tolerance is used).For solutions that vary significantly in magnitude, we would need to implement relative error control.This requires using a node-dependent tolerance where the tolerances δ A and δ R are user-defined, and δ A is included to cater for those occasions when |y i | ∼ 0.
7. Applying RKrvQz to a system of differential equations would require error control to be applied to each component of the system.This would lead to a value of h * i for each component -we would choose the smallest.We would also test the inequality (9) for each component; if at least one of the components failed the test, we would perform a quench in all components.
8. We anticipate that it should not be difficult to modify RKrvQz for implicit RK methods, which would be suitable for stiff problems.In this regard, RKr, RKv and RKz would all be implicit, A-stable RK methods.For example, we could use the well-known second-order Implicit Midpoint Rule, the fourth-order Kuntzmann-Butcher method, and the sixth-order Kuntzmann-Butcher method in place of RKr, RKv and RKz, respectively.We could denote such an algorithm by IRK24Q6, where the 'I' indicates 'implicit'.

Conclusion
We have developed a numerical algorithm for solving initial-value problems in ordinary differential equations, designated RKrvQz, that is capable of controlling both local and global errors in the numerical solution, in a stepwise manner.The algorithm uses local extrapolation to control the local error, and so-called quenching to retard the build-up of global error.Three Runge-Kutta methods are used in RKrvQz : RKr and RKv are used for local error control, and RKz is used for the estimation of global error and the quenching procedure.A numerical example with a rapidly increasing global error has demonstrated the effectiveness of the algorithm.In this exploratory paper, we have restricted our work to scalar problems, with absolute error control.The extension of RKrvQz to systems, and the incorporation of relative error control must be the subject of future research.