A class of high-order Runge-Kutta-Chebyshev stability polynomials

# A class of high-order Runge-Kutta-Chebyshev stability polynomials

Stephen O’Sullivan School of Mathematical Sciences, Dublin Institute of Technology, Kevin Street, Dublin 8, Ireland
###### Abstract

The analytic form of a new class of factorized Runge-Kutta-Chebyshev (FRKC) stability polynomials of arbitrary order is presented. Roots of FRKC stability polynomials of degree are used to construct explicit schemes comprising forward Euler stages with internal stability ensured through a sequencing algorithm which limits the internal amplification factors to . The associated stability domain scales as along the real axis. Marginally stable real-valued points on the interior of the stability domain are removed via a prescribed damping procedure.

By construction, FRKC schemes meet all linear order conditions; for nonlinear problems at orders above 2, complex splitting or Butcher series composition methods are required. Linear order conditions of the FRKC stability polynomials are verified at orders 2, 4, and 6 in numerical experiments. Comparative studies with existing methods show the second-order unsplit FRKC2 scheme and higher order (4 and 6) split FRKCs schemes are efficient for large moderately stiff problems.

###### keywords:
Stiff equations, Stability and convergence of numerical methods, Method of lines
###### Msc:
[2010] 65L04, 65L20, 65M20
journal: Journal of Computational Physics

## 1 Introduction

Runge-Kutta-Chebyshev methods are explicit numerical integration schemes with extended stability domains derived from the optimality properties of Chebyshev polynomials (Markov, 1890, 1892). These methods are commonly applied to moderately stiff systems of semi-discrete equations of the form

 w′=f(t,w), (1)

yielding an approximate solution at time , defined on a spatial mesh of spacing at points , with . Such systems arise naturally through application of the method of lines to parabolic systems. Runge-Kutta-Chebyshev methods may be broadly categorized as factorized or recursive in nature.

Factorized Runge-Kutta-Chebyshev methods are formed from a sequence of forward Euler stages. These methods were first suggested by Saulâev (1960); Guillou and Lago (1960) and were subsequently considered by Gentzsch and Schluter (1978) and van der Houwen (1996). They have been applied at first-order and extended to second-order via Richardson extrapolation by various authors (Alexiades et al., 1996; O’Sullivan and Downes, 2006, 2007; O’Sullivan and O’Sullivan, 2011). Based on a strategy proposed by Lebedev (2000), the DUMKA stability polynomials exist at orders 2, 3, and 4 (Medovikov, 1998).

Recursive Runge-Kutta-Chebyshev methods were first described by van Der Houwen and Sommeijer (1980) and rely on (three-term) recursion to generate a solution. They were introduced at second-order by van Der Houwen and Sommeijer (1980) and, subsequently, other second, third, and fourth-order methods have been developed (Verwer, 1996; Sommeijer et al., 1998; Abdulle and Medovikov, 2001; Abdulle, 2002; Martin-Vaquero and Janssen, 2009). We note that alternative approaches with second-order accuracy involving Legendre polynomials have recently been proposed by Meyer et al. (2012, 2014). At orders above 2, for both factorized and recursive methods, composition techniques relying on Butcher series theory (Hairer et al., 1993; Butcher, 2008) are typically used to satisfy the full set of order conditions (Medovikov, 1998; Abdulle, 2002).

This paper is organized as follows. In Section 2, the analytic form of the class of FRKC stability polynomials is presented. The construction of stable time-marching schemes based on the roots of these polynomials is outlined. Section 3 is given over to the derivation of the polynomial through consideration of associated recurrence relations. In Section 4, numerical tests are presented confirming the order and efficiency properties of FRKC methods. Conclusions are presented in Section 5.

## 2 High-order factorized Runge-Kutta-Chebyshev

### 2.1 General prescription

Eq. 1 may be written in autonomous form by appending to the vector of dependent variables for the system

 w′=f(w). (2)

Parentheses may be used in the remainder of this work to differentiate exponents from indices. We proceed by considering order extended stability explicit Runge-Kutta schemes over stages

 WL=W0+TL∑l=1alf(Wl−1), (3)

where corresponds to the approximate solution at time level , and yields at a time later. The timestep related to each stage is then given by .

The FRKC polynomial of rank , and degree , is given by

 BNM(z)=dN0+2N∑k=1dNkCkM(z), (4)

where denotes the the Chebyshev polynomial of the first kind of degree . The corresponding optimal real stability range is , where , , and (with rapidly converging to 1 as increases). In this limit, the polynomials generate , and of the optimal intervals for order respectively (see Van Der Houwen (1977) and Abdulle (2001) for estimates of the optimal values for ). The limiting step size is , where are the negative-definite eigenvalues for the Jacobian of Eq. 2. We note that the form of Eq. 4 is consistent with the known result that Chebyshev expansions of stability polynomials exist to arbitrary order (Bakker, 1971). Furthermore, following from a proposition by Lomax (1968), Riha (1972) confirmed the existence and uniqueness of optimal stability polynomials with local maxima with value unity. A full derivation of the FRKC polynomial expression given by Eq. 4 is provided in Section 3.

The order coefficients , which we refer to collectively as the order pattern, are determined by requiring that the (undamped) stability polynomial , consisting of shifted Chebyshev polynomials, satisfies the linear order conditions

 RNM(n)(0)=1,n=1,⋯,N. (5)

This requirement is met by solving the -dimensional linear system111The identity is useful here.

 ⎡⎢ ⎢ ⎢ ⎢⎣C(1)M(1)…C(1)NM(1)⋮⋱⋮C(N)M(1)…C(N)NM(1)⎤⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢⎣dN1⋮dNN⎤⎥ ⎥ ⎥⎦=⎡⎢ ⎢⎣(M2αM)1⋮(M2αM)N⎤⎥ ⎥⎦, (6)

coupled with the constraint

 dN0=1−2N∑k=1dNk. (7)

Following identification of the roots of the FRKC polynomial , the damped order scheme corresponding to Eq. 3 is determined via

 al=1M2αM 11−ζl. (8)

In order to ensure a stable scheme for small perturbations from the real axis in the spectrum of Eq. 2, it is necessary to introduce a suitable damping procedure. We find an effective prescription for the damped order scheme is given by

 al=1(1−ν)M2αM 1−μl1−(1−2μl)ζl, (9)

where the damping is parameterized by the small positive quantity , resulting in the real extent of the stability interval being reduced to . The value of is regulated by means of the reference damping parameter , such that maxima in along the real axis are scaled by approximately .

For the case , with , and for various values of , Fig. 1 illustrates the effect of the damping procedure. It is clear that the undamped polynomials are marginally stable at points on the interior of the stability domain along the real axis. (In fact, for sub-optimal , internal marginally stable points occur at locations for even values of , or locations for odd values of .) Examples of the order patterns for with are given in A.

The -tuple has cardinality and regulates the implementation of damping in the scheme while preserving the nominal order of accuracy. The values of are obtained by tuning the damped stability polynomial to meet the linear order conditions given in Eq. 5. We describe the procedure for the determination of the damping coefficients in Section 2.2.

### 2.2 Identification of damping parameter L-tuple

The elementary symmetric polynomial, , is defined as the sum of all possible products formed from unrepeated elements drawn from the first elements of an -tuple . The definition is extended by setting and . We associate the roots , in order of increasing real component , with the damping coefficients by cycling through the damping coefficients a total of times. Newton-Raphson iterations then converge rapidly to the linear order conditions given by Eq. 5. The effects of the damping procedure on the stability domain are shown in Fig. 1. The stage intervals given by Eq. 9 are complex in general, however, with , , the standard first-order super-timestepping scheme (Alexiades et al., 1996; O’Sullivan and Downes, 2006) is recovered with . For , either one or two values of have negative real parts.

The presented prescription implements conjugate pairs separately thereby necessitating full complex arithmetic. Other than some penalty in the additional computational demand required, we find no practical disadvantage to preserving this model of treating each factor as distinct. We note that Lebedev (2000, 1994) proposed a scheme which treats stages in pairs and, when applied to conjugate pairs, removes the need for complex arithmetic.

### 2.3 Internal stability

Schemes comprising a high number of stages are internally unstable if the sequencing of the stages is allowed to admit uncontrolled growth of numerical errors (Van Der Houwen, 1977; Lebedev and Finogenov, 1976; Verwer et al., 1990; Hairer and Wanner, 1996; Ketcheson et al., 2013). Lebedev and Finogenov (1973) first suggested sequencing stages to manage uncontrolled growth of internal instabilities (see also (Marchuk and Lebedev, 1986)). Here, we present a straightforward algorithm for sequencing stages which limits the maximum amplification factor of internal instabilities to .

We define , where are discrete values spanning the spectrum of Eq. 12. The -tuple is then ordered by holding the normed quantity

 ∥∥ ∥∥max⎛⎝l∏j=1vj,k,L∏j=l+1vj,k⎞⎠∥∥ ∥∥1 (10)

to a minimum value while is increased from 1 to . This procedure suppresses the growth of the internal stability functions , for , over , and provides excellent internal stability properties with high numbers of stages at low computational cost. In Fig. 2, we plot the maximum internal stability function for the test cases , with , and . The optimization may be enhanced by concentrating the points towards the bounds of the interval. (In this work a logistic function over a range of 15 is employed to generate the sample points.) We observe the maximum internal amplification factor scales approximately as . Hence, the internal stability properties are well within the acceptable limits of modern computing precision for any practical problem.

Consistent with these findings, we note that internal amplification factors of are quoted in the literature for RKC methods with stages (Hundsdorfer and Verwer, 2003), and furthermore, a quadratic dependence on stage number is suggested by Sommeijer et al. (1998). ROCK2 methods are reported to demonstrate amplification factors of at 200 stages by Hundsdorfer and Verwer (2003), suggesting internal instability growth rates 150 time larger than for RKC and FRKC2 schemes.

We note that the SERK scheme is also limited in stage number, requiring 600 digits of precision for 320 stages, albeit principally due to severely ill-conditioned matrix systems used in calculating the stability polynomials by means of the Remez algorithm (Martin-Vaquero and Janssen, 2009). A subsequent revision of the SERK methodology has demonstrated a stability range which is four times larger Martín-Vaquero et al. (2014).

## 3 Factorized Runge-Kutta-Chebyshev polynomial derivation

We proceed by considering the one-dimensional diffusion equation

 ∂w∂t=∂2w∂x2. (11)

The semi-discrete form of Eq. 11 may be written , where is a tridiagonal matrix with diagonal entries -2, subdiagonal entries 1, corresponding to a second-order central discretization of the spatial derivative. The eigenvalues of are negative with a maximum magnitude of 4. Application of the numerical scheme given by Eq. 3 yields

 wn+1=L∏l=1(I+τlh2D)wn. (12)

The FRKC polynomial may be derived by consideration of the canonical scheme given by Eq. 12 over an extended timestep , spanning time levels to , and consisting of segments, with each segment comprising stages. We write the solution state corresponding to as , and assume and ; more complex states may be constructed by superposition. The solution state corresponding to is then obtained from . To aid the following discussion, Fig. 3 is provided to graphically represent solution states at different segment levels for the particular case . A reference point value of the solution state , at spatial index , is shown as a black node.

To proceed, we assume schemes consisting of segments () are known which generate the solution states, . For , the solution state spans nodes from a given point profile . Successive states regenerate this pattern, but spanning nodes, with non-zero values interspersed by zero valued nodes. We refer to the sequence of patterns over increasing values of as a pattern flow. Illustrations of sample pattern flows are given in Fig. 3.

Using Eq. 8, the components of the states may be recast as . Over a single timestep, Eq. 12 then takes the simplified form

 ¯¯¯¯¯¯Wm=mN∏l=1¯¯¯¯¯Dml¯¯¯¯¯¯W0, (13)

where is a tridiagonal matrix with diagonal entries and subdiagonal entries 1/2. In terms of the elementary symmetric polynomials we have

 ¯¯¯¯¯¯Wmj=m∑l=0cmj,lσml, (14)

where are coefficients dependent on the scheme Eq. 13. By induction, these coefficients have the properties

 (15)

The ()-tuples fully determine through the roots of the associated polynomial defined by

 BNm(2)mN−1=m∑l=0(−1)lσml(z)m−l, (16)

where is a normalization factor. Hence, the -tuple is completely specified by .

### 3.1 Derivation

The ()-tuple is constructed from the elementary symmetric polynomials corresponding to the solution , where indicates the ordered elements , zero superscripts denote multiplicity, and . Through Eq. 16, the ()-tuple maps to the degree polynomial . Inserting into Eq. 14 and appealing to the properties of the coefficients , as given by Eq. 15, yields a direct correspondence to . Hence, we derive the association

 b\tinyR∑g=0(b\tinyRg)¯¯¯¯¯¯Wm−bj−b\tinyR+2g∼(2z)b\tinyRBNm−b(2)(m−b)N−1. (17)

By Eq. 13, the solution state generates a pattern scaled by a factor of 1/2 with respect to the pattern corresponding to the solution state . Hence, a recurrence relation generating the correct pattern comprising any weighted average of over available values of will yield a valid solution state .

We define a ray as any connection on a uniformly spaced graph which passes through nodes on every segment level , . The sum of the recurrence weightings over any ray terminating at must be unity if the ray originates at the origin of a pattern flow at , and zero otherwise. The coefficients of the Gaussian polynomials (), denoted , possess the required properties. In Fig. 3, rays are shown summing to unity and zero, with a list of weightings satisfying these properties for all possible rays for the particular case of . Defining , the primitive form of the recurrence relation for is

 ¯¯¯¯¯¯WLj=N∑k=1(−1)k+1Gk∑l=0[Pk]lq[Pk(12)kN¯¯¯¯¯¯WL−kNj−12Gk+l−(12)(P−k)N¯¯¯¯¯¯WL−(P−k)Nj−12Gk+l]+(12)PN¯¯¯¯¯¯WL−PNj, (18)

where is the degree of for . We note that the Gaussian polynomial possesses a unique representation as a summation of the binomial powers , for , given by

 [Pk]q=12Gk∑g=0rP,kgq12Gk−g(1+q2)g, (19)

where the coefficients follow the generating function

 ∞∑k=0∞∑g=0(−1)k(2)grP,kg(t)k(z)g=(1−t)N∏k=1(1+(t)2−2tCk). (20)

Then, using Eq. 19, we may recast Eq. 18 in the form

 ¯¯¯¯¯¯WLj=N∑k=1(−1)k+112Gk∑g=0rP,kgg∑l=0(gl)[Pk(12)kN¯¯¯¯¯¯WL−kNj−g+2l−(12)(P−k)N¯¯¯¯¯¯WL−(P−k)Nj−g+2l]+(12)PN¯¯¯¯¯¯WL−PNj. (21)

Applying the association given in Eq. 17 to the terms in Eq. 21, the recurrence relation for follows as

 BNM=N∑k=1(−1)k+1[BNM−k−BNM−P+k]12Gk∑g=0rP,kg(2z)g+BNM−P. (22)

We continue by noting that the generating function for derived from the recurrence relation given by Eq. 22 is

 ∞∑k=0(t)kBNk=∑2Nk=0(t)kbNk1−∑Nk=1(−1)k+1[(t)k−(t)(P−k)]∑12Gkg=0rP,kg(2z)g−(t)P, (23)

where are coefficients determined by the seed states of . Appealing to Eq. 20, the generating function derived from the recurrence relation given by Eq. 22 is

 ∞∑k=0(t)kBNk=bN01−t+2N∑k=1bNk(1−zt)1+(t)2−2tCk(1−t)N∏k=1(1+(t)2−2tCk), (24)

where are coefficients determined by the seed states of . The normalization has been imposed in order to fix the forms of the numerators in the separated fractions.

Noting that the generating function for is , we conclude that . Consideration of the particular case indicates a correspondence between and is required in order to match the required solution pattern and normalization properties. A general correspondence between and is established by considering successive values of , with , for the limiting case , . This completes the derivation of the analytic expression for the FRKC stability polynomial given by Eq. 4.

## 4 Tests

In this section numerical studies of two-dimensional two-species Brusselator diffusion-reaction problems are presented which confirm that high-order FRKC stability polynomials meet all relevant linear order conditions and that the derived factorized numerical schemes are both stable and efficient. Split schemes, denoted FRKCs, are obtained by means of complex splitting techniques: linear diffusion operators are treated via FRKC methods, while nonlinear reaction terms are integrated using standard Runge-Kutta techniques. The performance of the second-order accurate unsplit FRKC2 scheme is compared to second-order RKC and CVODE2 codes . Finally, comparisons are presented of higher order split FRKCs schemes (at orders 4 and 6), with fourth-order ROCK4 , and fifth-order CVODE schemes.

### 4.1 High-order splitting

FRKC stability polynomials satisfy linear order conditions to an arbitrary order of accuracy. This property may be exploited in solving numerical problems for semi-linear stiff systems of equations through operator splitting methods (Castella et al., 2009; Hansen and Ostermann, 2009, 2010; Dörsek and Hansen, 2014). We note that in the literature, the linear and nonlinear terms of reaction-diffusion models have been decoupled under a variety of numerical integration techniques including: splitting methods (McLachlan and Quispel, 2002, and previous references), Implicit-Explicit Runge-Kutta-Chebyshev (IMEX RKC) methods (Ascher et al., 1997; Shampine et al., 2006), PIROCK (Abdulle and Vilmart, 2013), and Local Linearization Runge-Kutta (LLRK) methods (Cruz et al., 2006; De la Cruz et al., 2013).

High-order splitting has been shown to give rise to an order reduction effect in some reaction-diffusion cases (Warnez and Muite, 2013). For Dirichlet and Neumann boundary conditions, splitting techniques may result in order reduction at boundaries (Hundsdorfer and Verwer, 2003; Hansen and Ostermann, 2009). It has also been observed that the full order is recovered on the interior of the computational domain when it is taken sufficiently far from the influence of the boundaries (Lubich and Ostermann, 1995; Lubich and Makridakis, 2013). Boundary conditions for the separate operator updates are necessary to avoid order reduction effectively, however, as yet, no consistent treatment exists (Hundsdorfer and Verwer, 1995).

Assuming Eq. 2 is linearized and split in the form , the solution over a timestep requires an approximation to the operator . High-order approximations may be obtained through appropriate choice of partial steps where

 wn+1=eTkJBeTkJ−1A⋯eTk3BeTk2AeTk1Bwn. (25)

Formally, with support from numerical studies (Castella et al., 2009; Blanes et al., 2013), the splitting scheme given by Eq. 25 may be may be extended to the semi-linear parabolic form of Eq. 2 given by

 w′=Aw+fB(w) (26)

by replacing the exponential operator with a step of the nonlinear equation over the interval . For reference, the complex splitting schemes used in this work are provided in Table 4.

### 4.2 Brusselator

The Brusselator (Lefever and Nicolis, 1971; Hairer et al., 1993) is a stiff nonlinear diffusion-reaction problem describing chemical kinetics of a tri-molecular chemical reaction. The test case considered here is a two-dimensional hybrid of the one- and two-dimensional Brusselator problems presented by Hairer et al. (1993), and Hairer and Wanner (1996), with governing equations given by

 ∂v/∂t = ϵ(∂2v/∂x12+∂2v/∂x22)+A−(B+1)v+wv2, ∂w/∂t = ϵ(∂2w/∂x12+∂2w/∂x22)+Bv−v2w, (27)

and initial conditions , . The initial state is therefore a simple perturbation of the equilibrium solution. The problem is configured with parameters , , and , and the solution is obtained at , or , on the domain , , under periodic boundary conditions.

### 4.3 Linear order conditions

The semi-discrete form of Eq. 27 may be written , where describes the discretization of the Laplacian with respect to and , and contains the reaction terms. Linear diffusion terms are integrated using FRKC methods and nonlinear reaction terms via standard techniques. The linear order properties of the FRKC stability polynomials are confirmed by considering the convergence rates of the approximated solution to the exact solution at as a function of step size.

For all presented results, we use , and the approximation . The number of grid points is 400 in each of the two spatial variables. For these parameters, the FRKC stability polynomials achieve approximately (), (), and () of the optimal intervals for respectively. With respect to the corresponding standard explicit Runge-Kutta schemes, these values represent a speedup in efficiency by factors of approximately for , and for both and . All polynomials are damped with damping parameter , reducing the stability domains’ real extents by factors of . Finally, in order to meet the specified solution time, timesteps are scaled by 0.9846, 0.9001, 0.7563 for respectively. Quadruple precision is used in all calculations. Results are presented in Table 1 where the and errors are shown over a range of resolutions at each considered value of . Fig. 4 illustrates the dependence of the errors on the number of timesteps, , for species . With the exception of the final point for the sixth-order integration, where machine precision is exceeded, all solutions are converging in good agreement with the nominal orders of accuracy (i.e. ). Fitting the