Relative Global Error Control in the RKQ Algorithm for Systems of Ordinary Di ff erential Equations

We generalize the RKrvQz algorithm to solve nonstiff initial-value problems in ordinary differential equations. The algorithm can now be applied to systems of nonstiff initial-value problems (IVPs) in ordinary differential equations, and both relative error and absolute error can be controlled, locally and globally. We demonstrate the algorithm by solving the simple harmonic oscillator for moderate and strict tolerances.


Introduction
Recently, we described the RKrvQz algorithm (Prentice, 2011) for solving nonstiff initial-value problems (IVPs) in ordinary differential equations (ODEs).This algorithm uses three explicit Runge-Kutta (RK) methods, of orders r, v and z, to control both local and global errors in a stepwise manner.In that paper, we considered control of absolute error only, neglecting relative error control, and we considered the application of the algorithm to scalar problems, rather than systems of ODEs.The motivation for developing RKrvQz is that, in local extrapolation, the RKv solution is not only used to estimate the local error in the RKr solution, but the RKv solution is propagated in the RKr method.Any global error in the RKv solution is thus also propagated in the resulting RKr solution.This global error can accumulate to the point where the RKr solution is globally inaccurate (relative to some desired level of accuracy), even though its local error has been controlled.RKrvQz represents an attempt to control the global error in the RKr solution in a stepwise manner, i.e. as the RK iteration proceeds, rather than using an 'after-the-fact' reintegration procedure.
x ∈ x 0 , x f y (x 0 ) = y 0 , that is to say, systems of ODEs, and we discuss the inclusion of relative error control in the algorithm.

Relevant Concepts, Terminology and Notation
The current paper is based almost entirely on our previous work (Prentice, 2011), but for ease of reference we will present here concepts, notation and terminology relevant to our discussion as it pertains to systems of ODEs.Also, we will designate our previous paper as PR1, since we will refer to it several more times in the text.To a large extent, the remainder of this section is excerpted from PR1, with appropriate modifications.Throughout the remainder of the paper, quantities in normal font are scalars, and quantities in boldface font are n × 1 vectors, except α and F r y , which are n × n matrices.Additionally, we refer the reader to Hairer et al (2000), Butcher (2003), Iserles (2009), Kincaid & Cheney (2002), LeVeque (2007), and many references therein, for discussions of Runge-Kutta methods.

Runge-Kutta methods
The most general definition of a Runge-Kutta (RK) method for systems is (1) Such a method is said to have m stages (the k q ), and each stage is an n × 1 vector.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.As indicated earlier, we will focus our attention on 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 (of course, F (x i , w i ) is an n × 1 vector).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 w r i+1 = w r i + h i F r x i , w r i . (2) We may regard RKr as being defined by F r , although it is understood that, for any r, there are, generally speaking, numerous possible choices for F r .We denote the jth component of F r by F r j .The superscripts in (2) are labels, not exponents.The stepsize h i is given by h and carries the subscript because it may vary from step to step.
Note that , wherein the first index in the subscript indicates iteration number, and the second index indicates component.

Error propagation
We define the global error in a numerical solution generated by RKr at x i+1 by and the local error at x i+1 by Note the use of the exact value y i in the bracketed term in (4).Again, the superscripts are labels.The errors Δ r i+1 and ε r i+1 are n × 1 vectors.
We have previously shown (Prentice, 2009) that where is the Jacobian of F r (x, y) evaluated at (x i , ξ i ) , and ξ i is a vector of constants arising in the residual term of a first-order Taylor expansion of F r x i , w r i about the point (x i , y i ) ; see (Prentice, 2009) for detail.Equation (5) provides a "master" 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 the n × 1 vector ε 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 δ (we will say more about the nature of this tolerance later).Moreover, assume that, using stepsize h i , we find for the In other words, the magnitude of the local error ε r i+1, j exceeds the desired tolerance.We remedy the situation by determining a new stepsize h * i, j from and we repeat the RK computation with this new stepsize.
We would carry out this process of determining a new stepsize for each component of ε r i+1 that exceeds δ in magnitude; each such process would yield a stepsize h * i, j ; we would then choose Often, we introduce a so-called 'safety factor' σ, as in where σ < 1, so that the new stepsize is slightly smaller than that given by ( 9).This is an attempt to cater for the possibility that β r i+1 may have been underestimated, due to the assumptions made in deriving ( 7).The choice of the value of σ is subjective, although a representative value is 0.8.Hence, we have On the other hand, if we find that the estimated error does not exceed the tolerance in any component, then no stepsize adjustment is necessary, and we proceed directly to the next step, using the already existing value of the stepsize.
Furthermore, since the higher-order solution w v i is available, we use w v i (in place of w r i in ( 10)) 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 ( 7).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, and we denote this algorithm by RKrv.

Absolute and relative error control
If the tolerance δ in ( 8) is a constant, then the form of error control is absolute, i.e. we are demanding that ε r i+1, j must be less than some absolute limit δ.Alternatively, we could demand This means that we require ε r i+1, j to be less than some limit, relative to the magnitude of w i+1, j ≈ y i+1, j .Of course, if w i+1, j is close to zero, then the corresponding stepsize h * i, j will be very small, and if w i+1, j = 0, then h * i, j cannot be computed at all.To counteract this possibility, we actually demand where δ A and δ R are known as the absolute and relative tolerances, respectively, and δ i+1 denotes a node-dependent tolerance.We then have two cases: , then the two cases are equivalent.Often, we use There is a good practical reason for implementing relative error control.If w i+1, j > 1, then (9) gives, with . So h * i, j determined in the relative sense is larger than h * i, j determined in the absolute sense, particularly if w i+1, j 1.This implies that fewer nodes x i would be required on x 0 , x f , resulting in a more efficient algorithm.
In (11), we have defined the node-dependent tolerance δ i+1, j .The local error control algorithm described in section 2.3 is easily adapted to cater for relative error control simply by replacing δ with δ i+1, j .Since δ i+1, j can differ for each component of w i+1 , we can write this tolerance as a vector Consequently, in the next section, when we write we mean that each component of ε r i+1 is less than the corresponding component of δ i+1 , and the notation means that there is at least one component of ε r i+1 that is greater than its corresponding component of δ i+1 .

The RKrvQz Algorithm, with Absolute and Relative Error Control, for Systems of ODEs
We now describe the RKrvQz algorithm for systems, incorporating relative error control.This is a generalization of the algorithm presented in PR1, and we will be economical in our discussion; the reader is referred to PR1 for detail.
We have three RK methods (RKr, RKv and RKz) at our disposal, with r < v z.Let w v i+1 denote the approximate solution at x i+1 obtained with RKv, and similarly for w z i+1 .Let w rv i+1 denote the approximate solution at x i+1 obtained with RKr, using w v i as its input.Vol. 3, No. 4; November 2011 We have, using (3) and ( 5), since r z.The notation α rv i in the above simply indicates that w v i has been used as input for RKr, and the form of α rv i is no different from α r i in ( 6).Local error control via local extrapolation, as described in sections 2.3 and 2.4, is carried out using w rv i+1 and w v i+1 (i.e. with RKrv), yielding a stepsize h * i .Note that here h * i is determined with (9), so there is at least one component j such that If it so happens that a new stepsize is not necessary then, for the purposes of what follows, we set h * i = h i−1 .Assuming now that a safety factor σ has been used, so that the new stepsize is σh * i , we have the solutions The LHS of ( 15) is the 'new' local error, arising from the use of σh * i as the stepsize.Now, if where the LHS is known from ( 14), then global error control is not necessary, and we proceed with the next RK iteration.However, if then it means that at least one component of Δ v i has become unacceptably large.We respond by setting because RKz is of much higher order than RKv, and we recalculate w rv i+1 and w v i+1 using h * i .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.We have referred to this process as quenching in PR1.
Note that the safety factor σ ensures that ε r i+1 is strictly less than δ i+1 in (15), so that the global error component α rv i Δ v i can be accommodated somewhat.The extent of this accommodation is determined by σ r+1 .Say σ = 0.8 and r = 3, so that σ r+1 = 0.41.Then, assuming that all components of β r i+1 σh * i r+1 and α rv i Δ v i have the same sign, α rv i Δ v i can be accommodated up to a magnitude of 0.59δ i+1 , before quenching is needed.
Finally, we emphasize that, at each node x i+1 , it is w rv i+1 that is presented as the solution to the IVP, since this is the numerical solution for which both local and global error control has been implemented.

Comments
(1) In PR1 we also considered the computation of Published by Canadian Center of Science and Education to be estimated.This is merely a refinement, and is not necessary for the implementation of RKrvQz; furthermore, it does require additional computational effort.However, such an estimate of α rv i Δ v i is likely to be very good, since v z, and together with (14), will provide a very good estimate of ε r i+1 , which may be more reliable than the estimate of ε r i+1 obtained from w rv i+1 − w v i+1 (the estimate of ε r i+1 is, of course, required for local error control via RKrv).
(2) We have used the same tolerance δ i+1 for both local and global error control.It is not necessary to do so although, since global error is, roughly speaking, an accumulation of local error, it is probably wise to ensure that the global tolerance is not smaller than the local tolerance.
(3) The algorithm presented in PR1 is obtained from the current algorithm when δ R = 0 and n = 1.

Numerical Example
A particularly suitable example is the simple harmonic oscillator which has solution Since the solution is oscillatory, there will be occasions when the solution is close to zero, which requires absolute error control, as per ( 14).Also, since the maximum magnitude of the solution is 1000, there will be occasions when relative error control is required.
We use the explicit methods RK3, RK4 and RK8, as referenced in PR1.We solve (18) on x ∈ [0, 20] with RK34 (local extrapolation only), and RK34Q8 (local extrapolation with global quenching), for the cases δ A = δ R = 10 −5 and δ A = δ R = 10 −10 .We refer to these cases as Case I and Case II, respectively.For both cases we use δ = 0.8.We show some performance parameters in Table 1.In this table we show the maximum global error in each component of ( 18) and the number of quenches required in each case.Note that here the magnitude of the global error |Δ i | is calculated as It is clear that RK34, despite the implementation of local extrapolation, does not satisfy the imposed tolerances, and in one instance has a maximum error almost 700 times larger than the desired tolerance.On the other hand, RK34Q8 always achieves the desired level of accuracy.Error curves are shown in Figures 1 and 2, wherein the effects of the quenching procedure are clear.In each figure, there are two plots; it is understood that these plots share a common legend and a common x-axis.In the RK34Q8 plots, the quenches occur at those values of x where the error exhibits a sharp decrease.
We have estimated the absolute local error ε 8 i+1 of the RK8 solution using Richardson extrapolation (see PR1 for details), even though this requires extra computational effort.We then estimate the global error in the RK8 solution using with Δ 8 0 = 0 and For Case I, we estimate that the maximum magnitude of either component in Δ 8 is 13 × 10 −12 ; the actual value is 9 × 10 −12 .
For Case I, we estimate the maximum magnitude of either component in Δ 8 to be 5 × 10 −12 , while the actual value is 3 × 10 −12 .In both cases, our estimate is good and slightly larger than the actual value.Of course, using (19), these estimates can be computed in situ, i.e. as the iteration proceeds.This global error is measure of the quality of the solution w 8 which is used for the quenching process.If one or more of the components of the estimated value of Δ 8 is considered to be too large relative to δ A and/or δ R , we propose that δ A and/or δ R be increased relative to the estimated value of Δ 8 .This amounts to reducing the level of accuracy imposed on the problem, but would still yield reliable, albeit less accurate, results.If one or more of the components of Δ 8 become comparable to δ A and/or δ R , then the entire quenching process will be compromised, because then w 8 is no more accurate than w 3 , which defeats the purpose of quenching.This strategy can be summarized as max j=1,2 where γ < 1, η A > 1 and η R > 1 are user-defined, and the index j indicates component.For example, in Case II (δ A = δ R = 10 −10 ), if we had set γ = 0.03, we would have estimated Δ 8 1 > 0.03δ A at x ∼ 10.4.We could then have made the adjustments δ A → 2δ A , δ R → 2δ R , say.On the remainder of the interval, Δ 8 j ≯ γδ A = 6 × 10 −12 , and so no further adjustments to the tolerances would have been needed.It is reasonable to believe that a tolerance of δ A = δ R = 2 × 10 −10 would have been considered acceptable, in these circumstances.

Conclusion
We have extended the functionality of the RKrvQz algorithm, which now can be applied to IVPs in the form of systems of ODEs, and for which relative and/or absolute tolerances on local and global error can be imposed.The simple harmonic oscillator has served as a useful example for demonstrating this updated version of RKrvQz.

Figure 1 .
Figure 1.Global error curves for both components of the oscillator problem, using RK34 and RK34Q8, for Case I

Table 1 .
Performance parameters for RK34 and RK34Q8 applied to the oscillator problem