Stepsize Selection in Explicit Runge-Kutta Methods for Moderately Stiff Problems

We present an algorithm for determining the stepsize in an explicit Runge-Kutta method that is suitable when solving moderately stiff differential equations. The algorithm has a geometric character, and is based on a pair of semicircles that enclose the boundary of the stability region in the left half of the complex plane. The algorithm includes an error control device. We describe a vectorized form of the algorithm, and present a corresponding MATLAB code. Numerical examples for Runge-Kutta methods of third and fourth order demonstrate the properties and capabilities of the algorithm.


Introduction
Stiff initial-value problems (IVPs) are often solved numerically using implicit A-stable Runge-Kutta (RK) methods.In such methods, there is no need to adjust the stepsize for the sake of stability, since the stability region of the method is the entire left half of the complex plane.These methods are particularly useful when the problem is very stiff.However, implementing an implicit RK method requires the solution of a nonlinear system of stage equations at each step, whereas an explicit RK method does not.It is, therefore, often feasible to use an explicit RK method to solve moderately stiff problems, wherein the stiffness constant λ is not too large.Often, the explicit RK pairs RK (3,4) [1] and RK (4,5) [2,3] are used to solve IVPs, so we will focus our attention on these methods.
The stability region of explicit RK methods is a bounded region S in the complex plane, and the stepsize h must be chosen such that, if the vector hλ is in the left half of the complex plane, then it lies within S.Moreover, λ can be complex, so that hλ does not necessarily lie along the real axis (even though h is always real).In this paper, we present a simple algorithm for determining h such that hλ lies within S, for any relevant λ and S pertaining to an explicit RK3 or RK4 method.

Relevant Concepts, Terminology and Notation
An m-dimensional IVP of the form where can be solved using an explicit RK method where i w denotes the approximate numerical solution at i x (i.e.

 
i i w y x  ), F is a function that defines the method, and From here on, the notation RKr indicates an explicit Runge-Kutta method of order r.
The stability function   R z of the RK method is obtained by applying the RK method to the Dahlquist equation to get, with , which can be written with z h  .As a simple example, consider the second-order method of Heun [4], which has Applying this method to the Dahlquist equation from which we identify Generally speaking, the stability function for explicit RK methods is a power series in z that represents an approximation to the exponential solution of the Dahlquist equation.Indeed, we see that in (10) is a low-order Taylor series for the exponential function For RK3 and RK4 we have [5]   and these are the stability functions that will interest us in the remainder of this paper.The region of stability of the RK method is defined as [6]    : and we denote the boundary of this region by ∂S.For an explicit RK method, ∂S is a closed contour in the complex plane (see Figure 1).We must consider complex values for z, because it is possible that the stiffness constant λ could be complex; this is particularly true for systems of ODEs, as in (1), where the stiffness constants are determined from the eigenvalues of the Jacobian   There are m eigenvalues, and those with negative real part are taken as the stiffness constants of the system.Now, assume 0   in (4).The exact solution is therefore expected to be an exponentially decreasing function of x.However, if h is such that   1 R z  , then the numerical solution (9) will increase with x, which is quite the opposite behaviour to what is expected.Furthermore, the numerical solution will increase without bound under iteration, whereas the exact solution tends to zero.This is referred to as an unstable solution, or instability with regard to stiff ODEs.It is therefore vital to choose h at each step of the RK method so that  lies within S.An algorithm for determining an appropriate stepsize h is the subject of the next and subsequent sections.

Theoretical Description of the Algorithm
We first consider the algorithm for a single stiffness constant (eigenvalue of J  where  is the magnitude of λ.Then  pendent on *    , and  .Hence, the smaller we choose ε, the closer the endpoint of hλ is to ∂S, particularly for large  . Note the following: (22) gives and if we demand that the relative difference * h h h  must be less than some value δ, then choosing ensures that this will be case.In other words, (24) provides a form of error control, since 1 r and δ are known a priori.
The algorithm is depicted in  .The segment ween arrowheads that ∂S cuts is the segment of length D, referred to earlier.
For a system of ODEs, in which there is more than one stiffness constant, we bet would apply the above algorithm for each such constant.This would yield a stepsize for each stiffness constant, and we would choose the minimum of these stepsizes as h in (3).

Implementation of the Algorithm
and let  denote the corresponding n vector of unit vectors, as i and hence, determine 1 NM ent-by is the N × M unit matrix, and ⊠ denotes the elem -element product of two matrices, sometimes product (the matrices must have f course).The structure of known as the Hadamard the same dimensions, o where We provide a MATLAB code nction for determining h for RK3.

5.
Some comments are in order: 1) It may occur that one of the , it should be clear that it can be applied to RK method for which semicircles with r ii 1 r and 2 r can be constructe (i.e. in the left half the complex plane, one of the semicircles contains S entirely, and the other is contained entirely within S).For methods of order greater than four, however, the stability function R(z) is method-specific so that the appropriate semicircles would also be method-specific.In our numerical examples in the next section, we will give semicircles that are universally applicable to all RK3 and RK4 methods.
3) We have not considered applying the algorithm to RK1 and RK2 because it does not seem possible to define a semicircle interior to S for these methods (see Fig- ure 1).Moreover, the low order of these methods probably mitigates against their use in solving moderately stiff IVPs, anyway.
4) The vectorized algorithm described above is not the only way to determine We could use the algorithm with 1 M  in a for-loop, providing a different λ for each pass through the loop, which would give a stepsize for each pass, and then taking the minimum of these stepsizes.The vectorized algorithm seems, to us, to be more elegant and may also be faster, generally speaking, al onstan though this might depend on computational platform.5) In principle, the algorithm can be applied for stiffness c ts of any magnitude, although RK3 and RK4 would most likely be used only for moderately stiff problems, wherein |λ|≲1000.

Numerical Examples
In our first example, we consider the stiffness constants which have magnitudes (rounded to nearest integer) of ly.The first of these lies ose to the real axis, the second is close to the diagonal, Applying the algorithm to the RK3 stability region, with (36) 1000, 647 and 910, respective cl and the third is close to the imaginary axis.

Conclusions
We have designed a ithm for determining stepsizes appropriate for stable solutions of stiff IVPs, when such solutions are computed using explicit RK methods.The gorithm, based on a pair of sem Clearly, the relative "error" in the stepsizes is very small.Of course, since it is proportional to tive, and possesses the facility to control the accuracy with which the stepsize is determined.The algorithm caters for complex stiffness constants, which can arise in IVP systems, and can be implemented in vectorized or for-loop form, with the former appearing to be slightly faster in terms of execution time.Features of the algorithm have been ably demonstrated with respect to the RK3 and RK4 methods.Indeed, we expect that the algorithm will be most useful when using RK3 and RK4 to solve moderately stiff problems, although it can easily be used with explicit RK methods of higher order.maller, which can be made even smaller by choosing ε s be done in a controlled way, using (24).rtheless, the calculations done here show that the ste determined by the algorithm are very close to al, even for a fairly large tolerance universally applicable to RK3 and RK4.

References
In our second example, we measure speed of computa d seco Neve psizes optim of .Note that the radii and 2 r used in these calculations are tion for different values of M, for both the vectorized algorithm and the for-loop algorithm mentione in #3 of our earlier comments.Our computational platform is MATLAB 6.5, Windows XP Pro, HP 30AA Mboard, Intel Celeron M 1.73 GHz, 1 MB L2 Cache, 2 GB RAM.Results are shown in Table 1.
In this table, all times are in nds, v indicates vec-

Figure 1 .
Figure 1.Stability regions for the indicated explicit RK methods.
sem uch that, in the left of the complex plane, 1 C is contained entirely within S, and S is contained entirely within 2 C , as shown in Figure 2. In this figure, S is the stabili egion of RK3 (for the sake of example), although the algorithm we describe here holds for the stability region of RK4, as well.Say λ is an eigenvalue of ty r J  , and λ lies in the left half of the complex plane (i.e.λ a stiffness constant).Define the unit vector is

Figure 3 .
Figure 3. Geometrical representation of the algorithm.The boundary of the stability region is indicated as ∂S.

Figure 3 .
In this figure, the curves labelled 1 r and 2 r are the semicircles, and ∂S indicates the boundary of the stability region.The arrow touching the 1 r semicircle indicates 1 r  , and the arrow touching the 2 r semicircle indicates 2 r  .The arrows between these two represent j z , as indicated.The arrow within S and closest to ∂S is It is clear that ∂S cuts 2 r c z , as shown.
we present a vectorized version of the atering for a system of ODEs which has more than one stiffness constant.The parameters 1 2 , r r and ε are input from which N and *  yields M distinct stiffness constants (note that M m  ).Let  denote a row vector containing these ess constants, as in M stiffn than the upper bounds in (39), as expected from (23).If we had desired a bound of 0