Factorized Runge-Kutta-Chebyshev Methods

TL;DR

This paper introduces the FRKC2 class of explicit, parallelizable Runge-Kutta-Chebyshev methods for large PDE systems, achieving high stability and accuracy with potential for extension to higher orders.

Contribution

It presents a new factorized scheme that is easy to implement, highly stable, and suitable for large-scale parallel computations, with discussions on extending to fourth-order accuracy.

Findings

Preserves 7 digits of accuracy at 16-digit precision.

Achieves acceleration factors over 6000 compared to standard methods.

Offers stability domains comparable or larger than existing schemes.

Abstract

The second-order extended stability Factorized Runge-Kutta-Chebyshev (FRKC2) class of explicit schemes for the integration of large systems of PDEs with diffusive terms is presented. FRKC2 schemes are straightforward to implement through ordered sequences of forward Euler steps with complex stepsizes, and easily parallelised for large scale problems on distributed architectures. Preserving 7 digits for accuracy at 16 digit precision, the schemes are theoretically capable of maintaining internal stability at acceleration factors in excess of 6000 with respect to standard explicit Runge-Kutta methods. The stability domains have approximately the same extents as those of RKC schemes, and are a third longer than those of RKL2 schemes. Extension of FRKC methods to fourth-order, by both complex splitting and Butcher composition techniques, is discussed. A publicly available implementation of the FRKC2 class of schemes may be obtained from maths.dit.ie/frkc

Figures

Tables

Table 1. Table 1 . FRKC4 order composition tableaux.

0
$c_2$	$a_1$
$c_3$	$a_1$	$a_2$
$c_4$	$a_1$	$a_2$	$a_3$
$⋮$	$⋮$	$⋮$	$⋮$	$\ddots$
$c_L$	$a_1$	$a_2$	$a_3$	$\cdots$	$a_{L-5}$
	$a_1$	$a_2$	$a_3$	$\cdots$	$a_{L-5}$	$a_{L-4}$

Equations

$$w^{\prime }=f(t,w).$$

$$w^{\prime }=Aw,$$

$$W^L=W^0+T\sum _{l=1}^L{a}_{l}f(W^{l-1})W^0=w^n,w^{n+1}=W^L.$$

$$R_M^N(TA)={\mathrm{e}}^{TA}+\mathcal{O}(TA{)}^{N+1}.$$

$$R_M^N{}^{(n)}(0)=1,n=1,\cdots ,N.$$

$$|R_M^N(z)|\le 1\forall z=T\lambda .$$

$$B_M^N(z)=d_0^N+2\sum _{k=1}^N{d}_k^N{C}_{kM}(z),$$

$$a_l=\frac{1}{M^2{\alpha }_M}\frac{1}{1-{\zeta }_l}.$$

$$a_l=\frac{1}{(1-\nu )M^2{\alpha }_M}\frac{1-{\mu }_l}{1-(1-2{\mu }_l){\zeta }_l},$$

$$\mathcal{Q}\equiv \max ({\mathcal{Q}}_{j,k}(x)),1\le j\le k\le L,x\in [-\beta ,0],$$

$${\mathcal{Q}}_{j,k}(x)=\prod _{l=j}^k|1+a_{l}x|.$$

$$\overline{\mathcal{Q}}\equiv {\Vert \max \left(\prod _{j=1}^l{v}_{j,k},\prod _{j=l+1}^{L/2}{v}_{j,k}\right)\Vert }_1.$$

$$W^{l+1}=W^l+T{a}_{l}f(W^l).$$

$${\hat{W}}^{l+1}={\hat{W}}^l+T\mathrm{Re}(a_l)\mathrm{Re}(f(W^l)).$$

$$\Vert err\Vert =\Vert \frac{|W^{l+1}-{\hat{W}}^{l+1}|}{TOL(1+\max (|W^{l+1}|,|{\hat{W}}^{l+1}|)}\Vert .$$

$$T^{n+1}=\left(\frac{SAFE}{\sqrt{\Vert er{r}^n\Vert }}\right)\left(\frac{T^n}{T^{n-1}}\right)\sqrt{\frac{\Vert er{r}^{n-1}\Vert }{\Vert er{r}^n\Vert }},$$

$$\partial v/\partial t$$

$$\partial w/\partial t$$

$$v(0,x)$$

$$w^{n+1}={\mathrm{e}}^{T_{k_J}B}{\mathrm{e}}^{T_{k_{J-1}}A}\cdots {\mathrm{e}}^{T_{k_3}B}{\mathrm{e}}^{T_{k_2}A}{\mathrm{e}}^{T_{k_1}B}{w}^n.$$

$$a\overline{a}(t)=\frac{1}{\alpha (t)}\sum \left(\sum _{i=0}^q\left(\genfrac{}{}{0ex}{}{q}{i}\right)\overline{a}(s_i(t))\prod _{z\in d_i(t)}a(z)\right),$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure0.pdf}}})$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 1.0pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure3.pdf}}}})$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 3.35728pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure8.pdf}}}})$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 1.0pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure15.pdf}}}})$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 1.0pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure22.pdf}}}})$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 3.35728pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure31.pdf}}}})$$

$$\displaystyle+2a(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 1.0pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure41.pdf}}}})\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 1.0pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure42.pdf}}}})+a(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 3.35728pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure43.pdf}}}})+\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 3.35728pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure44.pdf}}}}),$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 3.35728pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure45.pdf}}}})$$

$$\displaystyle a\overline{a}(\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{\lower 5.71457pt\hbox{{\leavevmode\includegraphics{}{}{TikzPictures/osullivan2016-figure54.pdf}}}})$$

Full text

A publicly available implementation of the FRKC2 class of schemes may be obtained from maths.dit.ie/frkc

1. Introduction

Factorized Runge–Kutta–Chebyshev (FRKC) methods are well suited to the numerical integration of problems where diffusion limits the efficiency of standard explicit techniques. In general, such systems of PDEs may be presented as semi-discrete ordinary differential equations of the form

[TABLE]

Here, the associated Jacobian is assumed to have negative eigenvalues lying close to the real axis, a good approximation for many systems of interest in astrophysical contexts.

The main use of extended stability Runge–Kutta (ESRK) methods is to fill the gap between unconditionally stable but operationally complex implicit methods, and simply implemented explicit schemes which suffer from stability constraints for stiff problems. ESRK methods are particularly useful for problems involving diffusion, where the work required by standard explicit techniques goes as the inverse square of the mesh spacing, while for extended stability methods it goes as the inverse mesh spacing. ESRK explicit schemes may be broadly divided into factorized and recursive types.

Factorized ESRK methods are particularly straightforward to implement at second-order, consisting solely of forward Euler steps. At orders above two, splitting methods or, alternatively, additional finishing stages are required for nonlinear problems. First considered by [23, 9, 8], factorized ESRK methods fell out of common usage for some time, until revived in 1996 as SuperTimeStepping [4] at first-order, and later extended to second-order applications in astrophysical simulations by means of Richardson extrapolation [21, 22]. DUMKA methods exist at third- and fourth-order [16].

The perennial problem with the factorized formalism has been that, when a very large number of stages is used, the internal amplification factors can easily drown numerical precision. Factorized methods have, as a result, largely taken a back seat to recursive methods which manage internal stability by mapping the three-term recurrence relations for orthogonal polynomials to the stability polynomials [25]. Recursive methods have been developed up to fourth-order [27, 24, 3, 2, 15, 17, 18].

In the following, a formulation of factorized methods is presented which has high internal stability, and is more straightforward to implement than recursive methods, and demonstrates comparable efficiency.

2. FRKC stability polynomials

Stability polynomials form the backbone of ESRK numerical schemes and encapsulate the linear stability properties. Linearising a system of semi-discrete ODEs

[TABLE]

the FRKC stability polynomial [20] is obtained by seeking a polynomial of degree $L$ which yields a forward Euler scheme of order $N$ for linear problems through its roots via

[TABLE]

The $M$-th stability polynomial of order $N$, with $L=MN$, determined via ${w}^{n+1}=R^{N}_{M}(TA){w}^{n}$, must match the first $N+1$ terms in the Taylor expansion of the evolution operator

[TABLE]

Equivalently, the linear order conditions may be expressed as constraints on the derivatives of the stability polynomial evaluated at zero:

[TABLE]

In addition, stability requires that the polynomial is bounded according to

[TABLE]

The objective is to determine a closed form for the polynomial such that the extent of the stability domain along the negative real axis $\beta$ is as great as possible.

It is shown in [20] that the FRKC stability polynomial of rank $N$, and degree $L$, is a sum of Chebyshev polynomials of the first kind given by

[TABLE]

where $C_{kM}$ is the Chebyshev polynomial of the first kind of degree $kM$, and the coefficients $d_{k}$ are determined through the linear order conditions given by Equation 5. The resultant scheme follows immediately from the roots of the polynomial, $\zeta_{l}$, with coefficients given by

[TABLE]

The dependency of $\beta$ on $L$ is presented in Figure 1 at second-order ($N=2$). While the optimal value for $\beta/L^{2}\beta_{\scriptsize{RK2}}$ is $0.41$ [26], where $\beta_{\scriptsize{RK2}}$ is the conventional second order Runge–Kutta stability limit, values of 0.330, 0.333, and 0.25 are obtained for FRKC2, RKC [25], and RKL2 [18] respectively.

2.1. Damping

Along the real axis on the interior of the stability domain of the stability polynomial, there are points which are marginally stable, as shown in Figure 2. This is remedied by introducing the damping parameter $\nu\equiv\nu_{0}/N$ via

[TABLE]

and again enforcing the order conditions given by Equation 5 via Newton-Raphson iteration over the parameters $\mu_{l}$, which consist of $N$ distinct values, each repeated $M$ times. As a result, the marginally stable maxima in $|R|$ along the real axis are scaled by $\sim(1-\nu_{0})$ at the expense of reducing the extent of the stability domain $\beta$ along the real axis by approximately $(1-\nu)$.

The damping process may also be used to make the scheme applicable to problems with small hyperbolic components, with Péclet numbers Pe\,\raisebox{-1.72218pt}{\stackrel{{\scriptstyle<}}{{\scriptstyle\sim}}}\,10. For the case $Pe\approx 10$, illustrated in Figure 2, there is a 25% loss in the extent of the stability domain observed.

2.2. Internal stability

Internal instability may be caused when the product of any of the possible sub-sequences of steps act to generate large values which drown out numerical precision. Following the idea of Lebedev [14, 13], but with a more effective approach, the timesteps are ordered to approximately minimise $\mathcal{Q}$, where

[TABLE]

is the maximum over the internal amplification factors defined by

[TABLE]

While the optimal value of $\mathcal{Q}$ is approximately $L^{2}$, for the purposes of constructing a stabilization algorithm, a practical upper bound of $10L^{2}$ is chosen. In order to maintain this bound, the estimated maximum amplification factor $\overline{\mathcal{Q}}$ is held to a minimum while $l$ runs from 1 to $L$, where

[TABLE]

The amplification factors $v_{j,\,k}$ are defined by $v_{j,\,k}=|1+a_{j}x_{k}|$, where $x_{k}\in[\overline{\beta},\,0]$ are $L$ uniformly distributed values over the reduced range $\overline{\beta}=(1-nC)\beta$, with $C=10^{-4}$. Initially, $n=1$, however, in a limited number of cases, the process is repeated with $n$ incremented until the required bound is satisfied. In this work, the mean value of $n$ for the second-order schemes was found to be 1.5. Figure 3 shows the realised values of $\mathcal{Q}$ obtained via the stabilization algorithm for second-order schemes. Preserving 7 digits for accuracy, a scheme consisting of $10^{4}$ stages is therefore theoretically viable in a numerical integration carried out to 16 digit precision.

2.3. Convex Monotone Property

The convex monotone property (CMP) is of relevance to problems with spatially varying diffusion coefficients [18]. Figure 4 shows solutions obtained for a problem considered in [18] describing two materials placed into contact with differing temperatures and diffusion coefficients. For Chebyshev polynomial-based schemes such as RKC and FRKC2, noise in the solutions associated with failure to meet the CMP is evident if the schemes are forced to take steps with a uniform number of stages. Figures 4(a) and 4(b) illustrate this for the cases $L=7$ and $L=8$ respectively. The RKL2 scheme maintains the CMP naturally with the result that the solution is smooth (Figure 4(c)). As shown in Figure 4(d), for FRKC2 with an identical number of function evaluations as used for Figure 4(b), but with an error-control procedure (discussed further in Section 3), the noise is also absent.

3. FRKC2 schemes

A system $w^{\prime}={f}(w)$ is assumed such that the associated Jacobian has an eigenvalue of maximum magnitude $|\lambda|_{\mbox{{{\scriptsize\rm{max}}}}}$. Then, given a numerical solution ${w}^{n}$ at some time index $n$, $L=2M$ stages ${W}^{l}$ (for $l=0,\,\cdots,\,L$) are evaluated such that ${W}^{0}={w}^{n}$ and a second-order accurate solution ${w}^{n+1}={W}^{L}$ is obtained a time $T$ later. The intermediate stages of the scheme are determined via the Euler steps

[TABLE]

Error control is straightforward since a first-order solution is available at no additional cost in function evaluations. This first-order solution ${\widehat{W}}^{L}$ is obtained by considering only the real parts of $a_{l}$ and ${f}({W}^{l})$. Setting ${\widehat{W}}^{0}={W}^{0}$,

[TABLE]

The error, scaled to a specified tolerance $TOL$, is estimated using

[TABLE]

If $\|err\|>1$, the step is rejected and retried with $T$ scaled by $SAFE/\sqrt{\|err\|}$. Otherwise, a predictive controller is used to determine the subsequent timestep calibrated to the required tolerance via

[TABLE]

using values of $T$ and $\|err\|$ from previous timesteps. $SAFE$ is a safety factor chosen with a value 0.8 here. The reader is referred to [11] for further details of error control procedures.

3.1. FRKC2 public code

A freely available C implementation of the second-order FRKC2 schemes may be accessed at maths.dit.ie/frkc. The files frkc2core.c and frkc2user.c provide the code for internal calculations required by the FRKC2 scheme and the code specific to the particular problem respectively. For a given value of $M$, up to 257, the extent of the stability domain along the real axis, $\overline{\beta}$, and the maximum realised internal amplification factor $\mathcal{Q}$ (see also Figure 3) are given on line $3M-2$ of frkc2arks.dat. The real and imaginary parts of $a_{l}$ are recorded on lines $3M-2$ and $3M$ respectively.

The default problem provided in frkc2user.c is a two-species reaction diffusion Brusselator

[TABLE]

which possesses a spectral radius of approximately $6400$ for a $200\times 200$ mesh. The initial state is a perturbation of the equilibrium solution [19]which is given by $v=1$, $w=3$. Figure 5 shows the number of steps required to attain a given error in the solution for FRKC2, ROCK2 [1], and RKC. It is evident that, for a given precision, there is little difference in the number of steps required by FRKC2 and ROCK2. At higher degrees of acceleration (fewer steps), the difference between the three schemes is negligible.

4. FRKC4 scheme

4.1. Complex splitting

Above second-order, nonlinear order conditions are present which require additional consideration. One approach, given a semi-linear parabolic (reaction-diffusion) equation of the form $w^{\prime}=Aw+f_{B}(w)$, is to split the nonlinear part $f_{B}(w)$, which is typically easily integrated, from the linear diffusion terms $Aw$. For orders above two, this requires complex timesteps [6, 5] and may be prescribed in the form

[TABLE]

Figure 6 shows the split FRKC4s scheme is competitive with ROCK4 [2]. However, in support of the splitting approach, it may be noted that ROCK4 suffers significantly from internal stability issues arising from the finishing stages required for the nonlinear order conditions (discussed further in Section 4.2) which effectively limits the scheme to about 150 stages.

4.2. Butcher composition

An alternative to complex splitting is Butcher composition. At fourth-order, as illustrated in Table 1, $L-4$ forward Euler steps are adopted from the FRKC stability polynomial with $N=4$ as the scheme $a$, and appropriate finishing stages for the scheme $\overline{a}$ are subsequently derived.

According to a theorem of Hairer & Wanner [10], given the B-series $a$, $b$, the composite scheme $ab$ is determined via

[TABLE]

where rooted trees $t$ are used to represent derivatives in Taylor series. The first summation in Equation 19 is over all $\alpha(t)$ different labelings of $t$, $s_{i}(t)$ is the subtree formed by the first $i$ indices, and $\{z\in d_{i}(t)\}$ is the difference set of subtrees formed by the remaining indices. The eight order conditions at fourth order are then given by

[TABLE]

Hence, for given scheme $a$, imposing the required order conditions on $a\overline{a}$ yields equations for $\overline{a}$, which are in turn easily solved for $\overline{a}$. The reader is referred to [10] for further details.

Figure 6 shows a comparison of the FRKC4 scheme based on composition methods with other schemes. The reference solution is provided by a fifth-order implicit preconditioned BDF solver from the CVODE numerical package [7]. The number of steps required for a given precision is comparable for ROCK4 and FRKC4 and somewhat greater than the split FRKC4s scheme. This difference may be ascribed to a loss of precision due to the accumulation of errors over the finishing stages [12].

5. Conclusions

FRKC extended stability Runge–Kutta (ESRK) schemes are shown to well-suited to the integration of large-scale problems governed by systems of PDEs where the explicit timescale is constrained by diffusion. An implementation of the FRKC2 second-order schemes, publicly available at maths.dit.ie/frkc, is presented. The fourth-order FRKC4 schemes are also presented with nonlinear order conditions addressed via both splitting and composition methods. These schemes are shown to be competitive with alternative ESRK methods.

5.1. Acknowledgments

I am grateful to organisers of Astronum2016 for the invitation to present this work in Monterey and to an anonymous referee for helpful comments.

6. References