∎

¹¹institutetext: Vu Thai Luan ²²institutetext: Department of Mathematics and Statistics, Mississippi State University
410 Allen Hall, 175 President’s Circle, Mississippi State, MS, 39762
²²email: [email protected]

Efficient exponential Runge–Kutta methods of high order: construction and implementation ^†^†thanks: This work has been supported in part by National Science Foundation through award NSF DMS–2012022.

Vu Thai Luan

Dedicated to Professor Alexander Ostermann on the occasion of his 60th birthday.

Abstract

Exponential Runge–Kutta methods have shown to be competitive for the time integration of stiff semilinear parabolic PDEs. The current construction of stiffly accurate exponential Runge–Kutta methods, however, relies on a convergence result that requires weakening many of the order conditions, resulting in schemes whose stages must be implemented in a sequential way. In this work, after showing a stronger convergence result, we are able to derive two new families of fourth- and fifth-order exponential Runge–Kutta methods, which, in contrast to the existing methods, have multiple stages that are independent of one another and share the same format, thereby allowing them to be implemented in parallel or simultaneously, and making the methods to behave like using with much less stages. Moreover, all of their stages involve only one linear combination of the product of $\varphi$ -functions (using the same argument) with vectors. Overall, these features make these new methods to be much more efficient to implement when compared to the existing methods of the same orders. Numerical experiments on a one-dimensional semilinear parabolic problem, a nonlinear Schrödinger equation, and a two-dimensional Gray–Scott model are given to confirm the accuracy and efficiency of the two newly constructed methods.

Keywords:

Exponential Runge–Kutta methods Exponential integrators Stiff PDEs Efficient implementation

MSC:

MSC 65L04 MSC 65M06 MSC 65N12

^†^†journal: BIT

1 Introduction

In this paper, we are concerned with the construction and implementation of new efficient exponential Runge–Kutta integrators for solving stiff parabolic PDEs. These PDEs, upon their spatial discretizations, can be cast in the form of semilinear problems

u^{\prime}(t)=Au(t)+g(t,u(t))=F(t,u(t)),\qquad u(t_{0})=u_{0},

(1)

where the linear part $Au$ usually causes stiffness. The nonlinearity $g(t,u)$ is assumed to satisfy a local Lipschitz condition in a strip along the exact solution.

Exponential Runge–Kutta methods are a popular class of exponential integrators HO10 , which have shown a great promise as an alternative to standard time integration solvers for stiff systems and applications in recent years, see e.g. HO05 ; LO12b ; LO14b ; LO13 ; LO14a ; Luan2014 ; LO16 ; L17 ; Michels2017 ; Ju2017 ; LD18 ; Luan18 ; Pieper2019a . The main idea behind these methods is to solve the linear portion of (1) exactly and integrate the remaining nonlinear portion explicitly based on a representation of the exact solution using the variation-of-constants formula.

A $s$ -stage explicit exponential Runge–Kutta (expRK) method HO05 applied to (1) can be reformulated (see LO12b ; LO14b ) as


$\displaystyle U_{ni}$	$\displaystyle=u_{n}+c_{i}h\varphi_{1}(c_{i}hA)F(t_{n},u_{n})+h\sum_{j=2}^{i-1}a_{ij}(hA)D_{nj},\ 2\leq i\leq s,$	(2a)
$\displaystyle u_{n+1}$	$\displaystyle=u_{n}+h\varphi_{1}(hA)F(t_{n},u_{n})+h\sum_{i=2}^{s}b_{i}(hA)D_{ni},$	(2b)

where

D_{ni}=g(t_{n}+c_{i}h,U_{ni})-g(t_{n},u_{n}),\qquad 2\leq i\leq s.

(3)

Here, $U_{ni}$ denote the internal stages that approximate $u(t_{n}+c_{i}h)$ using the time step size $h=t_{n+1}-t_{n}>0$ and nodes $c_{i}$ . By construction, the coefficients $a_{ij}(z)$ and $b_{i}(z)$ are usually linear combinations of the entire functions

\varphi_{k}(z)=\int_{0}^{1}{\rm e}\hskip 1.0pt^{(1-\theta)z}\frac{\theta^{k-1}}{(k-1)!}\,\text{d}\theta,\quad k\geq 1

(4)

and their scaled versions $\varphi_{k}(c_{i}z)$ .

A common approach that has been used to determine the unknown matrix functions $a_{ij}(hA)$ and $b_{i}(hA)$ is to expand them as $a_{ij}(hA)=\sum_{k\geq 0}\alpha^{(k)}_{ij}(hA)^{k}$ , $b_{i}(hA)=\sum_{k\geq 0}\beta^{(k)}_{i}(hA)^{k}$ (e.g. using classical Taylor series expansions) to obtain order conditions. Clearly, the boundedness of the remainder terms of these expansions (and thus the error terms) are dependent of $\|A\|$ . Due to stability reasons, such resulting methods might not be suitable for integrating stiff PDEs, which $A$ typically has a large norm or is even unbounded operator. These methods are thus usually referred as classical (non-stiffly accurate) expRK methods. Unlike this approach, in a seminal contribution HO05 , Hochbruck and Ostermann derived a new error expansion with the remainder terms that are bounded independently of the stiffness (i.e. not involving the powers of $A$ ), leading to stiff order conditions, which give rise to the construction of stiffly accurate expRK methods of orders up to four. Following this, in LO13 Luan and Ostermann developed a systematic theory of deriving stiff order conditions for expRK methods of arbitrary order, thereby allowing the construction of a fifth-order method in LO14b .

In view of the existing stiffly accurate expRK methods in the literature, we observe that they were derived based on a convergence result that requires weakening many of the stiff order conditions (in order to minimize the number of required stages $s$ and matrix functions used in each internal stages $U_{ni}$ ). As a result, their structures contain internal stages $U_{ni}$ that are dependent of the preceding stages, implying that such methods must be implemented by computing each of these stages sequentially. Also, the very last stages usually involve several different linear combinations of $\varphi_{k}(c_{i}hA)$ -functions (using different nodes $c_{i}$ in their arguments) acting on different sets of vectors. This would introduce additional computational effort for these stages. For more details, we refer to Sections 2 and 5.

Motivated by the observations above, in this work we show a stronger convergence result for expRK methods up to order five which requires weakening only one order condition (thereby could improve the stability and accuracy) and offers more degree of freedoms in solving order conditions. Using this result and inspired by our recent algorithm, $\mathtt{phipm\_simul\_iom}$ , proposed in Luan18 (which allows one to simultaneously compute multiple linear combinations of $\varphi$ - functions acting on a same set of vectors), we construct new methods of orders 4 and 5 which involve only one linear combination of $\varphi$ - functions for each stage and have multiple internal stages $U_{ni}$ that are independent of one another, thereby allowing them to be computed in parallel. Furthermore, one can derive these independent stages in a way that they share the same form of linear combination of $\varphi_{k}(c_{i}hA)$ - functions acting on the same set of vectors, allowing them to be implemented simultaneously (by one evaluation). While these independent states can be computed in parallel (as mentioned above) by any algorithm which approximates the action of (the linear combination of) $\varphi$ - functions, we note that the possibility to compute them simultaneously is a new feature that can be used with our algorithm $\mathtt{phipm\_simul\_iom}$ (other algorithms, e.g., that do not require the construction of Krylov subspaces, might not support computing these stages simultaneously). Overall, this makes the new methods to behave like methods using much less number of stages (even when compared to the existing methods of the same orders), meaning that they require much less number of evaluations for linear combinations of $\varphi$ - functions, and are thus more efficient.

The paper is organized as follows. In Section 2, we describe our motivation, propose new ideas, and review the existing expRK methods in the literature with respect to these ideas. Following this, in Section 3 we prove a stronger convergence result (Theorem 3.1) for expRK methods, which requires relaxing only one order condition. This allows us to construct more efficient methods in Section 4. In particular, we are able to derive two new families of fourth- and fifth- order stiffly accurate expRK methods called $\mathtt{expRK4s6}$ (4th-order 6-stage but requires 4 independent stage evaluation only) and $\mathtt{expRK5s10}$ (5th-order 10-stage but requires 5 independent stage evaluation only), respectively. In Section 5, we present details implementation of these two new methods, as well as the existing stiffly accurate expRK schemes of the same orders (for comparison). In the last section, numerical examples including one and two-dimensional stiff PDEs are presented to demonstrate the accuracy and efficiency of the two newly constructed expRK integrators.

2 Motivation and existing methods

In this section, we start our motivation by taking a closer look at an efficient way for implementing expRK methods (2). Then, we propose some ideas to derive more efficient methods with respect to this efficient implementation along with reviewing the current methods.

2.1 An efficient way of implementation

Clearly, each stage ( $U_{ni}$ or $u_{n+1}$ ) of (2) requires computing matrix functions of the form $\varphi_{k}(c_{i}hA)v_{k}$ ( $0<c_{i}\leq 1$ ), where $v_{k}$ is some vector (could be $F(t_{n},u_{n}),D_{ni}$ or a linear combination of these). Thanks to recent developments AH11 ; NW12 ; CKOS16 ; GaudreaultPudykiewicz16 , one can efficiently compute a linear combination of $\varphi$ -functions acting on a set of input vectors $V_{0},\ldots,V_{q}$

\varphi_{0}(M)V_{0}+\varphi_{1}(M)V_{1}+\varphi_{2}(M)V_{2}+\cdots+\varphi_{q}(M)V_{q},

(5)

where $M$ is some square matrix. This is crucial when implementing exponential integrators. Very recently, in Luan18 , we were able to improve the implementations presented in NW12 ; GaudreaultPudykiewicz16 , resulting in the routine $\mathtt{phipm\_simul\_iom}$ . The underlying method in this algorithm is the use of an adaptive time-stepping technique combined with Krylov subspace methods, which allows us to simultaneously compute multiple linear combinations of type (5) using different scaling factors $\rho_{1},\cdots,\rho_{r}$ of $M$ , i.e.,

$\displaystyle\varphi_{0}(\rho_{1}M)V_{0}+\varphi_{1}(\rho_{1}M)V_{1}$	$\displaystyle+\varphi_{2}(\rho_{1}M)V_{2}+\cdots+\varphi_{N}(\rho_{1}M)V_{q},$	(6)
	$\displaystyle\vdots$
$\displaystyle\varphi_{0}(\rho_{r}M)V_{0}+\varphi_{1}(\rho_{r}M)V_{1}$	$\displaystyle+\varphi_{2}(\rho_{r}M)V_{2}+\cdots+\varphi_{N}(\rho_{r}M)V_{q}.$

Now taking $M=hA$ and considering $\rho_{k}$ ( $1\leq k\leq r$ ) as nodes $c_{i}$ used in expRK methods immediately suggests that one can compute the following ( $s-1$ ) linear linear combinations

\varphi_{1}(c_{i}hA)V_{1}+\varphi_{2}(c_{i}hA)V_{2}+\ldots+\varphi_{N}(c_{i}hA)V_{q},\quad 2\leq i\leq s

(7)

simultaneously by using only one evaluation (i.e., one call to $\mathtt{phipm\_simul\_iom}$ ). Note that this requires the use of a same set of vectors $[V_{1},\ldots,V_{q}]$ for all the linear combinations in (7).

Motivated by this, we see that if a $s$ -stage expRK scheme (2) is constructed in such a way that each internal stage $U_{ni}$ has the form

U_{ni}=u_{n}+\varphi_{1}(c_{i}hA)V_{1i}+\varphi_{2}(c_{i}hA)V_{2i}+\ldots+\varphi_{N}(c_{i}hA)V_{qi},

(8)

which includes only one linear combination of $\varphi$ - functions using exactly node $c_{i}$ as an argument in all $\varphi_{k}$ functions, then the scheme will contain a total of $s$ such linear combinations ( $s-1$ for $U_{ni}$ and 1 for $u_{n+1}$ as (2b) can be always written in the form of (8) with $c_{i}=1$ ), thereby requiring $s$ evaluations only. Furthermore, since the sets of vectors $[V_{1i},V_{2i},\cdots,V_{qi}]$ in (8) are usually different for each $U_{ni}$ , (7) also suggests that the efficiency will be significantly increased if one could build such stages (or a group of) $U_{ni}$ of the form (8) that share the same format (i.e., having the same set of acting vectors $V_{1i}\equiv V_{1},\ldots,V_{qi}\equiv V_{q}$ ) or that are independent of one another. As this allows to compute such stages simultaneously by one evaluation or to implement them in parallel similarly to our construction of parallel exponential Rosenbrock methods LO16 ), it certainly reduces the total number of required evaluations and thus speedups the computing time.

With respect to these observations, we now review the existing expRK schemes in the literature. Since our focus is on stiff problems, we will discuss only on stiffly accurate expRK methods, meaning that they satisfy the stiff order conditions (see Section 3 below).

2.2 Existing schemes and remarks

In HO05 , expRK methods of orders up to four have been derived. For later reference, we name the second-order, the third-order, and the fourth-order methods in that work as $\mathtt{expRK2s2}$ , $\mathtt{expRK3s3}$ , and $\mathtt{expRK4s5}$ , respectively. In LO14b , we have constructed an expRK method of order five called $\mathtt{expRK5s8}$ . To discuss all of these schemes in terms of the implementation, we rewrite their internal stages $U_{ni}$ and $u_{n+1}$ as linear combinations of $\varphi$ - functions like (8) and display them as follows (Note that since the first-order method, the exponential Euler scheme $u_{n+1}=u_{n}+\varphi_{1}(hA)hF(t_{n},u_{n}),$ has no internal stage, we do not consider it here).

$\mathtt{expRK2s2}$ :

	$\displaystyle U_{n2}=u_{n}$	$\displaystyle+\varphi_{1}(c_{2}hA)c_{2}hF(t_{n},u_{n}),$		(9)
	$\displaystyle u_{n+1}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)\tfrac{1}{c_{2}}hD_{n2}.$		(9)

$\mathtt{expRK3s3}$ (a representative with $c_{2}\neq\tfrac{2}{3}$ ):

$\displaystyle U_{n2}=u_{n}$	$\displaystyle+\varphi_{1}(c_{2}hA)c_{2}hF(t_{n},u_{n}),$	(10)
$\displaystyle U_{n3}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{2}{3}hA)\tfrac{2}{3}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{2}{3}hA)\tfrac{4}{9c_{2}}hD_{n2},$
$\displaystyle u_{n+1}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)\tfrac{3}{2}hD_{n2}.$

$\mathtt{expRK4s5}$ (the only existing fourth-order stiffly accurate expRK method constructed by Hochbruck and Ostermann HO05 ):

$\displaystyle U_{n2}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{2}hA)\tfrac{1}{2}hF(t_{n},u_{n}),$	(11)
$\displaystyle U_{n3}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{2}hA)\tfrac{1}{2}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{1}{2}hA)hD_{n2},$
$\displaystyle U_{n4}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)h(D_{n2}+D_{n3}),$
$\displaystyle U_{n5}=u_{n}$	$\displaystyle+[\varphi_{1}(\tfrac{1}{2}hA)\tfrac{1}{2}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{1}{2}hA)\tfrac{1}{4}h(2D_{n2}+2D_{n3}-D_{n4})$
	$\displaystyle+\varphi_{3}(\tfrac{1}{2}hA)\tfrac{1}{2}h(-D_{n2}-D_{n3}+D_{n4})]+[\varphi_{2}(hA)\tfrac{1}{4}h(D_{n2}+D_{n3}-D_{n4})$
	$\displaystyle+\varphi_{3}(hA)h(-D_{n2}-D_{n3}+D_{n4})],$
$\displaystyle u_{n+1}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)h(-D_{n4}+4D_{n5})+\varphi_{3}(hA)h(4D_{n4}-8D_{n5}).$

$\mathtt{expRK5s8}$ (the only existing fifth-order stiffly accurate expRK method constructed by Luan and Ostermann LO14b ):

	$\displaystyle U_{n2}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{2}hA)\tfrac{1}{2}hF(t_{n},u_{n}),$
	$\displaystyle U_{n3}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{2}hA)\tfrac{1}{2}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{1}{2}hA)\tfrac{1}{2}hD_{n2},$
	$\displaystyle U_{n4}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{4}hA)\tfrac{1}{4}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{1}{4}hA)\tfrac{1}{8}hD_{n3},$
	$\displaystyle U_{n5}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{2}hA)\tfrac{1}{2}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{1}{2}hA)\tfrac{1}{2}h(-D_{n3}+4D_{n4})+\varphi_{3}(\tfrac{1}{2}hA)h(2D_{n3}-4D_{n4})$
	$\displaystyle U_{n6}=u_{n}$	$\displaystyle+\varphi_{1}(\tfrac{1}{5}hA)\tfrac{1}{5}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{1}{5}hA)\tfrac{1}{25}h(8D_{n4}-2D_{n5})+\varphi_{3}(\tfrac{1}{5}hA)\tfrac{1}{125}h(-32D_{n4}+16D_{n5}),$
	$\displaystyle U_{n7}=u_{n}$	$\displaystyle+[\varphi_{1}(\tfrac{2}{3}hA)\tfrac{2}{3}hF(t_{n},u_{n})+\varphi_{2}(\tfrac{2}{3}hA)h(\tfrac{-16}{27}D_{n5}+\tfrac{100}{27}D_{n6})+\varphi_{3}(\tfrac{2}{3}hA)h(\tfrac{320}{81}D_{n5}-\tfrac{800}{81}D_{n6})]$
		$\displaystyle+[\varphi_{2}(\tfrac{1}{5}hA)h(\tfrac{-20}{81}D_{n4}+\tfrac{5}{243}D_{n5}+\tfrac{125}{486}D_{n6})+\varphi_{3}(\tfrac{1}{5}hA)h(\tfrac{16}{81}D_{n4}-\tfrac{4}{243}D_{n5}-\tfrac{50}{243}D_{n6})],$
	$\displaystyle U_{n8}=u_{n}$	$\displaystyle+\big{[}\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)h(\tfrac{-16}{3}D_{n5}+\tfrac{250}{21}D_{n6}+\tfrac{27}{14}D_{n7})$
		$\displaystyle+\varphi_{3}(hA)h(\tfrac{208}{3}D_{n5}-\tfrac{250}{3}D_{n6}-27D_{n7})+\varphi_{4}(hA)h(-240D_{n5}+\tfrac{1500}{7}D_{n6}+\tfrac{810}{7}D_{n7})\big{]}$
		$\displaystyle+\big{[}\varphi_{2}(\tfrac{1}{5}hA)h(\tfrac{-4}{7}D_{n5}+\tfrac{25}{49}D_{n6}+\tfrac{27}{98}D_{n7})+\varphi_{3}(\tfrac{1}{5}hA)h(\tfrac{8}{5}D_{n5}-\tfrac{10}{7}D_{n6}-\tfrac{27}{35}D_{n7})$
		$\displaystyle+\varphi_{4}(\tfrac{1}{5}hA)h(\tfrac{-48}{35}D_{n5}+\tfrac{60}{49}D_{n6}+\tfrac{162}{245}D_{n7})\big{]}$
		$\displaystyle+\big{[}\varphi_{2}(\tfrac{2}{3}hA)h(\tfrac{-288}{35}D_{n5}+\tfrac{360}{49}D_{n6}+\tfrac{972}{245}D_{n7})+\varphi_{3}(\tfrac{2}{3}hA)h(\tfrac{384}{5}D_{n5}-\tfrac{480}{7}D_{n6}-\tfrac{1296}{35}D_{n7})$
		$\displaystyle+\varphi_{4}(\tfrac{2}{3}hA)h(\tfrac{-1536}{7}D_{n5}+\tfrac{9600}{49}D_{n6}+\tfrac{5184}{49}D_{n7})\big{]},$
	$\displaystyle u_{n+1}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)h(\tfrac{125}{14}D_{n6}-\tfrac{27}{14}D_{n7}+\tfrac{1}{2}D_{n8})$
		$\displaystyle+\varphi_{3}(hA)h(\tfrac{-625}{14}D_{n6}+\tfrac{162}{7}D_{n7}-\tfrac{13}{2}D_{n8})+\varphi_{4}(hA)h(\tfrac{1125}{14}D_{n6}-\tfrac{405}{7}D_{n7}+\tfrac{45}{2}D_{n8}).$

Remark 1

In view of the structures of these schemes, one can see that only the second- and third-oder schemes ( $\mathtt{expRK2s2}$ , $\mathtt{expRK3s3}$ ) have all $U_{ni}$ in the form (8). While $\mathtt{expRK2s2}$ requires one internal stage $U_{n2}$ , $\mathtt{expRK3s3}$ needs two internal stages with $U_{n3}$ depends on $U_{n2}$ , making these stages cannot be computed simultaneously. As for $\mathtt{expRK4s5}$ , to the best of our knowledge, this 5-stage scheme is the only existing fourth-order stiffly accurate expRK method. As seen, among its internal stages the three internal stages $U_{n2}$ , $U_{n3}$ , and $U_{n4}$ are of the form (8) but again their corresponding sets of vectors $[V_{ki}]$ are not the same ( $[V_{k2}]=[\tfrac{1}{2}hF(t_{n},u_{n})],[V_{k3}]=[\tfrac{1}{2}hF(t_{n},u_{n}),hD_{n2}],[V_{k4}]=[hF(t_{n},u_{n}),h(D_{n2}+D_{n3})]$ ), and because of (3), they are not independent of one another ( $U_{n4}$ and $U_{n3}$ depend on their preceding stages). Therefore one needs 3 sequential evaluations for computing these three stages. Also, we note that the last internal stage $U_{n5}$ depends on all of its preceding stages and involves two different linear combinations of $\varphi_{k}$ - functions with different scaling factors $c_{5}=\tfrac{1}{2}$ and $c_{4}=1$ , namely, $\sum_{k}\varphi_{k}(\tfrac{1}{2}hA)V_{k}$ and $\sum_{k}\varphi_{k}(hA)W_{k}$ (grouped in two different brackets $[\ ]$ ), which has to be implemented by 2 separate evaluations. The final stage $u_{n+1}$ depends on $U_{n4}$ and $U_{n5}$ . As a result, this scheme must be implemented in a sequential way, which requires totally 6 evaluations for 6 different linear combinations. Similarly, to the best of our knowledge, $\mathtt{expRK5s8}$ is also the only existing fifth-order stiffly accurate expRK methods. From the construction of this scheme LO14b , one needs 8 stages. Among them, the first five internal stages are of the form (8). We note, however, that the last two internal stages $U_{n7}$ and $U_{n8}$ involves 2 and 3 different linear combinations (grouped in different brackets $[\ ]$ ) of $\varphi_{k}$ - functions (with different scaling factors) acting on different sets of vectors. And each stage ( $U_{ni}$ or $u_{n+1}$ ) depends on all the preceding stages (except for the first stage $U_{n2}$ ). Thus, this scheme must be also implemented in a sequential way (it also does not have any group of internal stages that can be computed simultaneously). Clearly, it requires totally 11 evaluations (11 different linear combinations of $\varphi$ functions).

Remark 2

The resulting structures of the expRK schemes discussed in Remark 1 can be explained by taking a closer look at their constructions presented in HO05 ; LO14b . Namely, these methods have been analyzed and derived by using a weakened convergence result, i.e., weakening of many order conditions in order to minimize the number of required stages $s$ and the number of matrix functions in each internal stage $U_{ni}$ . Specifically, for fourth-order methods (e.g., $\mathtt{expRK4s5}$ ) 4 out of 9 order conditions have to be relaxed and for fifth-order methods (e.g., $\mathtt{expRK5s8}$ ) 9 out of 16 order conditions have to be relaxed. As a trade off, each stage of these methods depends on the preceding stages (thus the resulting schemes must be implemented by computing each stage sequentially) and the very last stages usually involve different linear combinations of $\varphi_{k}$ -functions (with several different nodes $c_{i}$ as scaling factors) acting on not the same set of vectors, which then require additional sequential evaluations. For more details, see Section 4 below.

3 Stiff order conditions and convergence analysis

Inspired by the motivation and remarks in Section 2, we next present a stronger convergence result which later allows a construction of new efficient methods of high order. For this, we first recall the stiff order conditions for expRK methods up to order 5 (see LO12b ; LO14b ).

3.1 Stiff order conditions for methods up to order 5

Let $\tilde{e}_{n+1}=\hat{u}_{n+1}-u(t_{n+1})$ denote the local error of (2), i.e., the difference between the numerical solution $\hat{u}_{n+1}$ obtained by (2) after one step starting from the ‘initial condition’ $u(t_{n})$ and the corresponding exact solution $u(t_{n+1})$ of (1) at $t_{n+1}$ .

To simplify the notation in this section we set $f(t)=g(t,u(t))$ as done in HO05 , and additionally denote $G_{k,n}=D^{k}g(t_{n},u(t_{n}))$ be the $k$ -th partial Fréchet derivative (with respect to $u$ ) evaluated at $u(t_{n})$ . Our results in LO12b (Sect. 4.2) or LO13 (Sect. 5.1) showed that

	$\displaystyle\tilde{e}_{n+1}$	$\displaystyle=h^{2}\psi_{2}(hA)\,f^{\prime}(t_{n})+h^{3}\psi_{3}(hA)\,f^{\prime\prime}(t_{n})+h^{4}\psi_{4}(hA)\,f^{\prime\prime\prime}(t_{n})+h^{5}\psi_{5}(hA)\,f^{(4)}(t_{n})$		(12)
		$\displaystyle\quad+\mathbf{R}_{n}+\mathcal{O}(h^{6})$		(12)

with the remaining terms

$\displaystyle\mathbf{R}_{n}$	$\displaystyle=h^{3}\sum_{i=2}^{s}b_{i}G_{1,n}\psi_{2,i}\,f^{\prime}(t_{n})+h^{4}\sum_{i=2}^{s}b_{i}G_{1,n}\psi_{3,i}\,f^{\prime\prime}(t_{n})+h^{4}\sum_{i=2}^{s}b_{i}G_{1,n}\sum_{j=2}^{i-1}a_{ij}G_{1,n}\psi_{2,j}\,f^{\prime}(t_{n})$
	$\displaystyle+h^{4}\sum_{i=2}^{s}b_{i}c_{i}G_{2,n}\big{(}u^{\prime}(t_{n}),\psi_{2,i}\,f^{\prime}(t_{n})\big{)}+h^{5}\sum_{i=2}^{s}b_{i}G_{1,n}\psi_{4,i}\,f^{\prime\prime\prime}(t_{n})$
	$\displaystyle+h^{5}\sum_{i=2}^{s}b_{i}G_{1,n}\sum_{j=2}^{i-1}a_{ij}G_{1,n}\psi_{3,j}\,f^{\prime\prime}(t_{n})+h^{5}\sum_{i=2}^{s}b_{i}G_{1,n}\sum_{j=2}^{i-1}a_{ij}G_{1,n}\sum_{k=2}^{j-1}a_{jk}G_{1,n}\psi_{2,k}\,f^{\prime}(t_{n})$	(13)
	$\displaystyle+h^{5}\sum_{i=2}^{s}b_{i}G_{1,n}\sum_{j=2}^{i-1}a_{ij}c_{j}G_{2,n}\big{(}u^{\prime}(t_{n}),\psi_{2,j}\,f^{\prime}(t_{n})\big{)}+h^{5}\sum_{i=2}^{s}b_{i}c_{i}G_{2,n}\big{(}u^{\prime}(t_{n}),\psi_{3,i}\,f^{\prime\prime}(t_{n})\big{)}$
	$\displaystyle+h^{5}\sum_{i=2}^{s}b_{i}c_{i}G_{2,n}\big{(}u^{\prime}(t_{n}),\sum_{j=2}^{i-1}a_{ij}G_{1,n}\psi_{2,j}\,f^{\prime}(t_{n})\big{)}+h^{5}\sum_{i=2}^{s}\frac{b_{i}}{2!}G_{2,n}\big{(}\psi_{2,i}\,f^{\prime}(t_{n}),\psi_{2,i}\,f^{\prime}(t_{n})\big{)}$
	$\displaystyle+h^{5}\sum_{i=2}^{s}b_{i}\frac{c^{2}_{i}}{2!}G_{2,n}\big{(}u^{\prime\prime}(t_{n}),\psi_{2,i}\,f^{\prime}(t_{n})\big{)}+h^{5}\sum_{i=2}^{s}b_{i}\frac{c^{2}_{i}}{2!}G_{3,n}\big{(}u^{\prime}(t_{n}),u^{\prime}(t_{n}),\psi_{2,i}\,f^{\prime}(t_{n})\big{)}.$

Here, (and from now on) we use the abbreviations $a_{ij}=a_{ij}(hA)$ , $b_{i}=b_{i}(hA),\varphi_{j,i}=\varphi_{j}(c_{i}hA)$ and


$\displaystyle\psi_{j}(hA)$	$\displaystyle=\sum_{i=2}^{s}b_{i}\frac{c^{j-1}_{i}}{(j-1)!}-\varphi_{j}(hA),\quad j\geq 2$	(14a)
$\displaystyle\psi_{j,i}$	$\displaystyle=\psi_{j,i}(hA)=\sum_{k=2}^{i-1}a_{ik}\frac{c^{j-1}_{k}}{(j-1)!}-c^{j}_{i}\varphi_{j,i}.$	(14b)

Requiring a local error truncation $\tilde{e}_{n+1}=\mathcal{O}(h^{6})$ results in the stiff order conditions for methods of order up to 5, which are displayed in Table 1 below.

Table 1: Stiff order conditions for explicit exponential Runge–Kutta methods up to order 5. The variables

Z

J

K

L

denote arbitrary square matrices, and

B

an arbitrary bilinear mapping of appropriate dimensions. The functions

\psi_{l}

and

\psi_{k,l}

are defined in (14).

No.	Stiff order condition	Order
1	$\psi_{2}(Z)=0\Longleftrightarrow\sum_{i=2}^{s}b_{i}(Z)c_{i}=\varphi_{2}(Z)$	2
2	$\psi_{3}(Z)=0\Longleftrightarrow\sum_{i=2}^{s}b_{i}(Z)\frac{c_{i}^{2}}{2!}=\varphi_{3}(Z)$	3
3	$\sum_{i=2}^{s}b_{i}(Z)J\psi_{2,i}(Z)=0$	3
4	$\psi_{4}(Z)=0\Longleftrightarrow\sum_{i=2}^{s}b_{i}(Z)\frac{c_{i}^{3}}{3!}=\varphi_{4}(Z)$	4
5	$\sum_{i=2}^{s}b_{i}(Z)J\psi_{3,i}(Z)=0$	4
6	$\sum_{i=2}^{s}b_{i}(Z)J\sum_{j=2}^{i-1}a_{ij}(Z)J\psi_{2,j}(Z)=0$	4
7	$\sum_{i=2}^{s}b_{i}(Z)c_{i}K\psi_{2,i}(Z)=0$	4
8	$\psi_{5}(Z)=0\Longleftrightarrow\sum_{i=2}^{s}b_{i}(Z)\frac{c_{i}^{4}}{4!}=\varphi_{5}(Z)$	5
9	$\sum_{i=2}^{s}b_{i}(Z)J\psi_{4,i}(Z)=0$	5
10	$\sum_{i=2}^{s}b_{i}(Z)J\sum_{j=2}^{i-1}a_{ij}(Z)J\psi_{3,j}(Z)=0$	5
11	$\sum_{i=2}^{s}b_{i}(Z)J\sum_{j=2}^{i-1}a_{ij}(Z)J\sum_{k=2}^{j-1}a_{jk}(Z)J\psi_{2,k}(Z)=0\quad$	5
12	$\sum_{i=2}^{s}b_{i}(Z)J\sum_{j=2}^{i-1}a_{ij}(Z)c_{j}K\psi_{2,j}(Z)=0$	5
13	$\sum_{i=2}^{s}b_{i}(Z)c_{i}K\psi_{3,i}(Z)=0$	5
14	$\sum_{i=2}^{s}b_{i}(Z)c_{i}K\sum_{j=2}^{i-1}a_{ij}(Z)J\psi_{2,j}(Z)=0$	5
15	$\sum_{i=2}^{s}b_{i}(Z)B\big{(}\psi_{2,i}(Z),\psi_{2,i}(Z)\big{)}=0$	5
16	$\sum_{i=2}^{s}b_{i}(Z)c^{2}_{i}L\psi_{2,i}(Z)=0$	5

3.2 A stronger convergence result

The convergence analysis of exponential Runge–Kutta methods is usually performed in the framework of analytic semigroups on a Banach space $X$ with the following assumptions (see e.g. HO05 ; LO14b ):

Assumption 1. The linear operator $A$ is the infinitesimal generator of an analytic semigroup ${\rm e}\hskip 1.0pt^{tA}$ on $X$ . This implies that

\|{\rm e}\hskip 1.0pt^{tA}\|_{X\leftarrow X}\leq C,\quad t\geq 0

(15)

and consequently $\varphi_{k}(hA)$ , the coefficients $a_{ij}(hA)$ and $b_{i}(hA)$ of the method are bounded operators. Furthermore, the following stability bound (see (HO05, , Lemma 1))

\left\|hA\sum_{j=1}^{n}{\rm e}\hskip 1.0pt^{jhA}\right\|_{X\leftarrow X}\leq C

(16)

holds uniformly for all $n\geq 1$ and $h>0$ with $0<nh\leq T-t_{0}$ .

Assumption 2 (for high-order methods). The solution $u:[t_{0},T]\to X$ of (1) is sufficiently smooth with derivatives in $X$ and $g:[t_{0},T]\to X$ is sufficiently often Fréchet differentiable in a strip along the exact solution. All occurring derivatives are assumed to be uniformly bounded.

Let $e_{n+1}=u_{n+1}-u(t_{n+1})$ denote the global error at time $t_{n+1}$ . In LO14b , we have shown that $e_{n}$ satisfies the recursion

e_{n}=h\sum_{j=0}^{n-1}{\rm e}\hskip 1.0pt^{(n-j)hA}\mathcal{K}_{j}(e_{j})e_{j}+\sum_{j=0}^{n-1}{\rm e}\hskip 1.0pt^{jhA}\tilde{e}_{n-j},

(17)

where $\mathcal{K}_{j}(e_{j})$ are bounded operators on $X$ .

Motivated by Remark 2, we now give a stronger convergence result (compared to those results given in HO05 ; LO14b ) in the sense that it requires relaxing only one order condition.

Theorem 3.1

(Convergence) Let the initial value problem (1) satisfy Assumptions 1–2. Consider for its numerical solution an explicit exponential Runge–Kutta method (2) that fulfills the order conditions of Table 1 up to order $p$ ( $2\leq p\leq 5$ ) in a strong form with the exception that only one condition $\psi_{p}(Z)=0$ holds in a weakened form, i.e., $\psi_{p}(0)=0$ . Then, the method is convergent of order $p$ . In particular, the numerical solution $u_{n}$ satisfies the error bound

\|u_{n}-u(t_{n})\|\leq Ch^{p}

(18)

uniformly on compact time intervals $t_{0}\leq t_{n}=t_{0}+nh\leq T$ with a constant $C$ that depends on $T-t_{0}$ , but is independent of $n$ and $h$ .

Proof

The proof can be carried out in a very similar way as done in (LO14b, , Theorem 4.2). In view of (12) and (13) and employing the assumptions of Theorem 3.1 on the order conditions, we have $\mathbf{R}_{n}=0$ and thus

\tilde{e}_{n+1}=h^{p}\big{(}\psi_{p}(hA)-\psi_{p}(0)\big{)}G_{p-1,n}+h^{p+1}\mathbf{S}_{n},

(19)

where $G_{p-1,n}$ is defined in Section 3.1 and $\mathbf{S}_{n}$ involves the terms multiplying $h^{p+1}$ and higher order in (12) (clearly, $\|\mathbf{S}_{n}\|\leq C$ ). Inserting (19) (with index $n-j-1$ in place of $n$ ) into (17) and using the fact that there exists a bounded operator $\tilde{\psi}_{p}(hA)$ such that $\psi_{p}(hA)-\psi_{p}(0)=\tilde{\psi}_{p}(hA)hA$ yields

e_{n}=h\sum_{j=0}^{n-1}{\rm e}\hskip 1.0pt^{(n-j)hA}\mathcal{K}_{j}(e_{j})e_{j}+h^{p}\sum_{j=0}^{n-1}hA{\rm e}\hskip 1.0pt^{jhA}\tilde{\psi}_{p}(hA)G_{p-1,n-j-1}+h^{p+1}\mathbf{S}_{n-j-1}.

(20)

Using (15), (16) and an application of a discrete Gronwall lemma shows (18).∎

With the result of Theorem 3.1 in hand, we are now ready to derive more efficient methods. In particular, we will solve the system of stiff order conditions of Table 1 in the context of Theorem 3.1. It turns out that for methods of high order this will require an increase in the number of stages $s$ . However, we will have more degree of freedoms for constructing our desired methods as seen in Section 4 below. In addition, by relaxing only one order condition, we expect methods resulted from Theorem 3.1 to have better stability (and thus may be more accurate) when integrating stiff systems (see Section 6).

4 Derivation of new efficient exponential Runge–Kutta methods

In this section, we will derive methods which have the following features: (i) containing multiple internal stages $U_{ni}$ that are independent of each other (henceforth called parallel stages) and share the same format (thereby allowing them to be implemented in parallel); (ii) involving less number of evaluations of the form (8) when compared to the existing methods of the same orders (thus behaving like methods that use fewer number of stages $s$ ).

We first start with methods of order $p\leq 3$ . When solving order conditions for these methods (requiring at least $s=2$ and $s=3$ for second- and third-order methods, respectively), one can easily show that it is not possible to fulfill the desired feature (ii), particularly when comparing with $\mathtt{expRK2s2}$ (order 2, 2-stage) and $\mathtt{expRK3s3}$ (order 3, 3-stage) mentioned in Section 2. We omit the details. Therefore, we will focus on the derivation of new methods of higher orders, namely, orders 4 and 5.

4.1 A family of fourth-order methods with parallel stages

Deriving methods of order 4 requires solving the set of 7 stiff order conditions 1–7 in Table 1. First, we discuss on the required number of stages $s$ . It is shown in (HO05, , Sect.5.3) that $s=5$ is the minimal number of stages required to construct a family of fourth-order methods which satisfies conditions 1–3 in the strong sense and conditions 4–7 in the weakened form (relaxing $b_{i}(Z)$ as $b_{i}(0)$ ). In other words, with $s=5$ it is not possible to fulfill the order conditions in the context of Theorem 3.1, which requires only condition 4 holds in a weakened form $\psi_{4}(0)=0$ or equivalently $\sum_{i=2}^{s}b_{i}(0)\frac{c_{i}^{3}}{3!}=\varphi_{4}(0)=1/24$ . Therefore, we consider $s=6$ . In this case, conditions 1, 2, and the weakened condition 4 are


$\displaystyle b_{2}c_{2}+b_{3}c_{3}+b_{4}c_{4}+b_{5}c_{5}+b_{6}c_{6}$	$\displaystyle=\varphi_{2},$	(21a)
$\displaystyle b_{2}c^{2}_{2}+b_{3}c^{2}_{3}+b_{4}c^{2}_{4}+b_{5}c^{2}_{5}+b_{6}c^{2}_{6}$	$\displaystyle=2\varphi_{3},$	(21b)
$\displaystyle b_{2}(0)c^{3}_{2}+b_{3}(0)c^{3}_{3}+b_{4}(0)c^{3}_{4}+b_{5}(0)c^{3}_{5}+b_{6}(0)c^{3}_{6}$	$\displaystyle=6\varphi_{4}(0)=1/4,$	(21c)

and conditions 3, 5, 7 and 6 are


	$\displaystyle b_{2}J\psi_{2,2}+b_{3}J\psi_{2,3}+b_{4}J\psi_{2,4}+b_{5}J\psi_{2,5}+b_{6}J\psi_{2,6}=0,$		(22a)
	$\displaystyle b_{2}J\psi_{3,2}+b_{3}J\psi_{3,3}+b_{4}J\psi_{3,4}+b_{5}J\psi_{3,5}+b_{6}J\psi_{3,6}=0,$		(22b)
	$\displaystyle b_{2}c_{2}K\psi_{2,2}+b_{3}c_{3}K\psi_{2,3}+b_{4}c_{4}K\psi_{2,4}+b_{5}c_{5}K\psi_{2,5}+b_{6}c_{6}K\psi_{2,6}=0,$		(22c)
	$\displaystyle b_{3}Ja_{32}J\psi_{2,2}+b_{4}J(a_{42}J\psi_{2,2}+a_{43}J\psi_{2,3})+b_{5}J(a_{52}J\psi_{2,2}+a_{53}J\psi_{2,3}+a_{54}J\psi_{2,4})$		(22d)
	$\displaystyle+b_{6}J(a_{62}J\psi_{2,2}+a_{63}J\psi_{2,3}+a_{64}J\psi_{2,4}+a_{65}J\psi_{2,5})=0.$

We now solve these order conditions. We note from (14b) that

\psi_{2,i}=\sum_{j=2}^{i-1}a_{ij}c_{j}-c^{2}_{i}\varphi_{2,i},\quad\psi_{3,i}=\sum_{j=2}^{i-1}a_{ij}\frac{c^{2}_{j}}{2!}-c^{3}_{i}\varphi_{3,i}

(23)

and thus $\psi_{2,2}=-c^{2}_{2}\varphi_{2,2}\neq 0,\quad\psi_{3,2}=-c^{3}_{2}\varphi_{3,2}\neq 0$ (since $c_{2}\neq 0$ ). Using (23), one can infer that either $\psi_{2,3}$ or $\psi_{3,3}$ must be nonzero as well (if both are zero then $a_{32}=\dfrac{c^{2}_{3}}{c_{2}}\varphi_{2,3}=\dfrac{2c^{3}_{3}}{c^{2}_{2}}\varphi_{3,3}$ , which is impossible since $c_{3}>0$ and $\{\varphi_{2},\varphi_{3}\}$ are linearly independent). This strongly suggests that $b_{2}=b_{3}=0$ in order to later fulfill (22) in the strong sense with arbitrary square matrices $J$ and $K$ . Next, we further observe that if $b_{4}\neq 0$ one may need both $\psi_{2,4}=\psi_{3,4}=0$ (which solves $a_{42}\neq 0$ , $a_{43}\neq 0$ ). However, this makes the second term in (22d) to be nonzero which is then very difficult to satisfy (22d) in the strong form. Putting together, it requires that $b_{2}=b_{3}=b_{4}=0$ . Using this sufficient condition we can easily solve (21) to get

b_{5}=\dfrac{-c_{6}\varphi_{2}+2\varphi_{3}}{c_{5}(c_{5}-c_{6})},\quad b_{6}=\dfrac{-c_{5}\varphi_{2}+2\varphi_{3}}{c_{6}(c_{6}-c_{5})}

for any choice of distinct nodes $c_{5},c_{6}>0$ , satisfying the condition

c_{5}=\dfrac{4c_{6}-3}{6c_{6}-4}.

(24)

Since $b_{5},b_{6}\neq 0$ , we must enforce $\psi_{2,5}=\psi_{3,5}=0$ and $\psi_{2,6}=\psi_{3,6}=0$ to satisfy conditions (22a)–(22c). Using (23), this leads to the following 2 systems of two linear equations


$\displaystyle a_{52}c_{2}+a_{53}c_{3}+a_{54}c_{4}$	$\displaystyle=c^{2}_{5}\varphi_{2,5},$	(25a)
$\displaystyle a_{52}c^{2}_{2}+a_{53}c^{2}_{3}+a_{54}c^{2}_{4}$	$\displaystyle=2c^{3}_{5}\varphi_{3,5},$	(25b)

and


$\displaystyle a_{62}c_{2}+a_{63}c_{3}+a_{64}c_{4}+a_{65}c_{5}$	$\displaystyle=c^{2}_{6}\varphi_{2,6},$	(26a)
$\displaystyle a_{62}c^{2}_{2}+a_{63}c^{2}_{3}+a_{64}c^{2}_{4}+a_{65}c^{2}_{5}$	$\displaystyle=2c^{3}_{6}\varphi_{3,6}.$	(26b)

To satisfy conditions (22d), we further enforce $a_{52}=a_{62}=0$ (since $\psi_{2,2}\neq 0$ ), which immediately solves (25) for coefficients (with $c_{3}\neq c_{4}$ )

a_{53}=\dfrac{-c_{4}c^{2}_{5}\varphi_{2,5}+2c^{3}_{5}\varphi_{3,5}}{c_{3}(c_{3}-c_{4})}\neq 0,\quad a_{54}=\dfrac{-c_{3}c^{2}_{5}\varphi_{2,5}+2c^{3}_{5}\varphi_{3,5}}{c_{4}(c_{4}-c_{3})}\neq 0,

(27)

and thus we also need $\psi_{2,3}=\psi_{2,4}=0$ (since $\psi_{2,5}=0$ ), which gives


$\displaystyle a_{32}$	$\displaystyle=\dfrac{c^{2}_{3}}{c_{2}}\varphi_{2,3},$	(28a)
$\displaystyle a_{42}c_{2}+a_{43}c_{3}$	$\displaystyle=c^{2}_{4}\varphi_{2,4}.$	(28b)

After fulfilling all the required order conditions in (21)–(22), we see from (26) and (28b) that either $a_{42}$ or $a_{43}$ and one of the coefficients among $a_{63},\ a_{64},\ a_{65}$ can be taken as free parameters. We now use them to construct parallel stages. Guided by (27) and (28a), we choose $a_{43}=0$ to make $U_{n4}$ is independent of $U_{n3}$ so that both these stages only depend on $U_{n2}$ , and choose $a_{65}=0$ to make $U_{n6}$ is independent of $U_{n5}$ so that both these stages only depend on the two preceding stages $U_{n3},U_{n4}$ (since $a_{52}=a_{62}=0$ ). From this we determine the remaining coefficients

a_{42}=\dfrac{c^{2}_{4}}{c_{2}}\varphi_{2,4},\quad a_{63}=\dfrac{-c_{4}c^{2}_{6}\varphi_{2,6}+2c^{3}_{6}\varphi_{3,6}}{c_{3}(c_{3}-c_{4})},\quad a_{64}=\dfrac{-c_{3}c^{2}_{6}\varphi_{2,6}+2c^{3}_{6}\varphi_{3,6}}{c_{4}(c_{4}-c_{3})}.

(29)

Putting altogether and rearranging terms in $U_{ni},\ u_{n+1}$ as linear combinations of $\varphi$ functions, we obtain the following family of 4th-order 6-stage methods (with the pairs of parallel stages $\{U_{n3},U_{n4}\}$ and $\{U_{n5},U_{n6}\}$ ), which will be called $\mathtt{expRK4s6}$ :


$\displaystyle U_{n2}=u_{n}$	$\displaystyle+\varphi_{1}(c_{2}hA)c_{2}hF(t_{n},u_{n}),$	(30a)
$\displaystyle U_{n,k}=u_{n}$	$\displaystyle+\varphi_{1}(c_{k}hA)c_{k}hF(t_{n},u_{n})+\varphi_{2}(c_{k}hA)\tfrac{c^{2}_{k}}{c_{2}}hD_{n2},\quad\hskip 42.67912ptk=3,4$	(30b)
$\displaystyle U_{n,j}=u_{n}$	$\displaystyle+\varphi_{1}(c_{j}hA)c_{j}hF(t_{n},u_{n})+\varphi_{2}(c_{j}hA)\tfrac{c^{2}_{j}}{c_{3}-c_{4}}h\big{(}\tfrac{-c_{4}}{c_{3}}D_{n3}+\tfrac{c_{3}}{c_{4}}D_{n4}\big{)}$
	$\displaystyle+\varphi_{3}(c_{j}hA)\tfrac{2c^{3}_{j}}{c_{3}-c_{4}}h\big{(}\tfrac{1}{c_{3}}D_{n3}-\tfrac{1}{c_{4}}D_{n4}\big{)},\quad\quad\ \hskip 56.9055ptj=5,6$	(30c)
$\displaystyle u_{n+1}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)\tfrac{1}{c_{5}-c_{6}}h\big{(}\tfrac{-c_{6}}{c_{5}}D_{n5}+\tfrac{c_{5}}{c_{6}}D_{n6}\big{)}$
	$\displaystyle+\varphi_{3}(hA)\tfrac{2}{c_{5}-c_{6}}h\big{(}\tfrac{1}{c_{5}}D_{n5}-\tfrac{1}{c_{6}}D_{n6}\big{)}.$	(30d)

For the numerical experiments given in Section 6, we choose $c_{2}=c_{3}=\frac{1}{2},c_{4}=\frac{1}{3}$ , $c_{6}=\frac{1}{3}$ which gives $c_{5}=\frac{5}{6}$ due to (24).

Remark 3

(A comparison with $\mathtt{expRK4s5}$ ). As seen, $\mathtt{expRK4s6}$ is resulted from weakening only condition 4 of Table 1 instead of weakening four conditions 4–7 as in the derivation of $\mathtt{expRK4s5}$ . While the 5-stage method $\mathtt{expRK4s5}$ requires 6 sequential evaluations in each step (as mentioned in Section 2), the new fourth-order 6-stage method $\mathtt{expRK4s6}$ requires only 4 sequential evaluations, making it to behave like a 4-stage method. This is due to the fact $\mathtt{expRK4s6}$ has the pairs of parallel stages $\{U_{n3},U_{n4}\}$ and $\{U_{n5},U_{n6}\}$ and all $U_{ni}$ within these pairs have the same format, i.e., same (one) linear combination of $\varphi_{k}(c_{i}hA)v_{k}$ , allowing them to be computed in parallel or simultaneously (see Section 5).

4.2 A family of fifth-order methods with parallel stages

Constructing fifth-order exponential Runge-Kutta methods needs much more effort as one has to solve 16 order conditions in Table 1. As mentioned in Section 2, the only existing method of order 5 in the literature is $\mathtt{expRK5s8}$ (see LO14b ) which requires $s=8$ stages. Like $\mathtt{expRK4s5}$ , this method does not have any parallel stages and must be implemented in a sequential way. It also does not satisfy the assumption on the order conditions stated in Theorem 3.1. Indeed, it was constructed by fulfilling conditions 1–7 in the strong form and weakening conditions 8–16 (9 out of 16 order conditions) with $b_{i}(0)$ in place of $b_{i}(Z)$ . This resulted in the last two internal stages $U_{n7}$ and $U_{n8}$ that involve several different linear combinations of $\varphi_{k}(c_{i}hA)v_{k}$ (with different scalings $c_{6},c_{7},c_{8}$ of $hA$ ), for which additional computational efforts are required to compute those stages (as shown in Section 2).

Therefore, to derive a method based on Theorem 3.1 which later allows us to derive parallel stages schemes with $U_{ni}$ involving only one linear combination of $\varphi_{k}(c_{i}hA)v_{k}$ , we have to increase $s\geqslant 9$ . To make it easier for readers to follow, we consider $s=10$ first and later employ the similar procedure to show that it is not possible to fulfill condition 11 of Table 1 in the strong form (and thus not satisfying Theorem 3.1) with $s=9$ .

a) The case $s=10$ : Similarly to the derivation presented in Subsection 4.1, using (23), it strongly suggests $b_{2}=b_{3}=b_{4}=b_{5}=b_{6}=b_{7}=0$ in order to solve conditions 3, 5, 9, 7, 16, 13, and 15 in their strong form. Using this, these conditions now read as


	$\displaystyle b_{8}J\psi_{2,8}+b_{9}J\psi_{2,9}+b_{10}J\psi_{2,10}=0,$		(31a)
	$\displaystyle b_{8}J\psi_{3,8}+b_{9}J\psi_{3,9}+b_{10}J\psi_{3,10}=0,$		(31b)
	$\displaystyle b_{8}J\psi_{4,8}+b_{9}J\psi_{4,9}+b_{10}J\psi_{4,10}=0,$		(31c)
	$\displaystyle b_{8}c_{8}K\psi_{2,8}+b_{9}c_{9}K\psi_{2,9}+b_{10}c_{10}K\psi_{2,10}=0,$		(31d)
	$\displaystyle b_{8}c^{2}_{8}L\psi_{2,8}+b_{9}c^{2}_{9}L\psi_{2,9}+b_{10}c^{2}_{10}L\psi_{2,10}=0,$		(31e)
	$\displaystyle b_{8}c_{8}K\psi_{3,8}+b_{9}c_{9}K\psi_{3,9}+b_{10}c_{10}K\psi_{3,10}=0,$		(31f)
	$\displaystyle b_{8}B(\psi_{2,8},\psi_{2,8})+b_{9}B(\psi_{2,9},\psi_{2,9})+b_{10}B(\psi_{2,10},\psi_{2,10})=0,$		(31g)

respectively. And conditions 1, 2, 4, and 8 (weakened form) become


$\displaystyle b_{8}c_{8}+b_{9}c_{9}+b_{10}c_{10}$	$\displaystyle=\varphi_{2},$	(32a)
$\displaystyle b_{8}c^{2}_{8}+b_{9}c^{2}_{9}+b_{10}c^{2}_{10}$	$\displaystyle=2\varphi_{3},$	(32b)
$\displaystyle b_{8}c^{3}_{8}+b_{9}c^{3}_{9}+b_{10}c^{3}_{10}$	$\displaystyle=6\varphi_{4},$	(32c)
$\displaystyle b_{8}(0)c^{4}_{8}+b_{9}(0)c^{4}_{9}+b_{10}(0)c^{4}_{10}$	$\displaystyle=24\varphi_{5}(0)=1/5.$	(32d)

Solving (32) gives


$\displaystyle b_{8}$	$\displaystyle=\dfrac{c_{9}c_{10}\varphi_{2}-2(c_{9}+c_{10})\varphi_{3}+6\varphi_{4}}{c_{8}(c_{8}-c_{9})(c_{8}-c_{10})},$	(33a)
$\displaystyle b_{9}$	$\displaystyle=\dfrac{c_{8}c_{10}\varphi_{2}-2(c_{8}+c_{10})\varphi_{3}+6\varphi_{4}}{c_{9}(c_{9}-c_{8})(c_{9}-c_{10})},$	(33b)
$\displaystyle b_{10}$	$\displaystyle=\dfrac{c_{8}c_{9}\varphi_{2}-2(c_{8}+c_{9})\varphi_{3}+6\varphi_{4}}{c_{10}(c_{10}-c_{8})(c_{10}-c_{9})}$	(33c)

where $c_{8},c_{9}$ , and $c_{10}$ are distinct and positive nodes satisfying the algebraic equation

\dfrac{c_{8}+c_{9}+c_{10}}{4}-\dfrac{c_{8}c_{9}+c_{8}c_{10}+c_{9}c_{10}}{3}+\dfrac{c_{8}c_{9}c_{10}}{2}=\dfrac{1}{5}.

(34)

Clearly, $b_{8},b_{9},b_{10}\neq 0$ so one has to enforce

\psi_{2,j}=\psi_{3,j}=\psi_{4,j}=0\ (j=8,9,10)

(35)

to satisfy (31) in the strong sense with arbitrary square matrices $J,K,L$ and $B$ . Next, we consider conditions 6 and 10 taken into account that $b_{i}=0$ ( $i=2,\cdots,7$ ) and (35), which can be now simplified as

\sum_{j=2}^{7}(b_{8}Ja_{8j}+b_{9}Ja_{9j}+b_{10}Ja_{10j})J\psi_{m,j}=0\ \ (m=2,3),

(36)

respectively. In order to satisfy the strong form of (36) one needs

a_{8j}=a_{9j}=a_{10j}=0\ (j=2,3,4)

(37)

(this is again due to (23)) and

\psi_{2,j}=\psi_{3,j}=0\ (j=5,6,7).

(38)

With (37), we note that $U_{n8},U_{n9},U_{n10}$ are independent of the internal stages $U_{n2},U_{n3},U_{n4}$ . Taking into all the requirements above, one can easily see that conditions 12 and 14 are now automatically fulfilled. Therefore, the only remaining condition to satisfy is condition 11.

Before working with condition 11, we first solve (35) using (37). For this, we observe that several coefficients $a_{ij}$ can be considered as free parameters. To have $U_{n8},U_{n9},U_{n10}$ are independent of each other, we choose

a_{98}=a_{10,8}=a_{10,9}=0.

(39)

The resulting systems of linear equations from (35) is then solved with the unique solution

a_{ij}=\dfrac{c^{2}_{i}c_{k}c_{l}\varphi_{2,i}-2c^{3}_{i}(c_{k}+c_{l})\varphi_{3,i}+6c^{4}_{i}\varphi_{4,i}}{c_{j}(c_{j}-c_{k})(c_{j}-c_{l})},\quad i=8,9,10;\ j,k,l\in\{5,6,7\},\ j\neq k\neq l

(40)

(i.e., $c_{5},c_{6},c_{7}>0$ are distinct nodes).
We now use $b_{i}=0$ ( $i=2,\cdots,7$ ), (35), (37), (38), and (39) to simplify condition 11 as

\sum_{i=8}^{10}b_{i}J\sum_{j=5}^{7}a_{ij}J\big{(}a_{j2}J\psi_{2,2}+a_{j3}J\psi_{2,3}+a_{j4}J\psi_{2,4}\big{)}=0.

(41)

Since $b_{8},b_{9},b_{10}\neq 0$ , coefficients $a_{ij}$ in (40) ( $i\in\{8,9,10\},j\in\{5,6,7\}$ ) are also nonzero, and that $\psi_{2,2}\neq 0$ , we must enforce

a_{j2}=0\ (j=5,6,7),\ \text{i.e.,}\ a_{52}=a_{62}=a_{72}=0

(42)

and require that

\psi_{2,3}=\psi_{2,4}=0

(43)

in order to satisfy (41) in the strong sense. Note, because of (42), one could not require $a_{53}=0$ or $a_{54}=0$ ( $j=5$ ) in (41) or both as this does not agree with the requirement $\psi_{2,5}=\psi_{3,5}=0$ in (38) (in other words, the linear system of equations displayed in (25) represented for this requirement has no solution). This justifies the requirement (43).

Finally, we solve (43) and (38) for the remaining coefficients $a_{ij}$ . When solving (43) (see (28)), we choose $a_{43}=0$ to have $U_{n4}$ is independent of $U_{n3}$ . This gives

a_{32}=\dfrac{c^{2}_{3}}{c_{2}}\varphi_{2,3},\quad a_{42}=\dfrac{c^{2}_{4}}{c_{2}}\varphi_{2,4}.

(44)

When solving (38) (using (42)), we choose $a_{65}=a_{75}=a_{76}=0$ to have $U_{n5},U_{n6},U_{n7}$ are independent of each other. This results in the following 6 coefficients:

a_{ij}=\dfrac{-c^{2}_{i}c_{k}\varphi_{2,i}+2c^{3}_{i}\varphi_{3,i}}{c_{j}(c_{j}-c_{k})},\quad i=5,6,7;\ j,k\in\{3,4\},\ j\neq k

(45)

(i.e., $c_{3},c_{4}>0$ are distinct nodes).

Inserting all the obtained coefficients $a_{ij}$ and $b_{i}$ into $U_{ni},\ u_{n+1}$ and rewriting these stages as linear combinations of $\varphi$ functions, we obtain the following family of 5th-order 10-stage methods (with the groups of parallel stages $\{U_{n3},U_{n4}\}$ , $\{U_{n5},U_{n6},U_{n7}\}$ , and $\{U_{n8},U_{n9},U_{n10}\}$ ) which will be called $\mathtt{expRK5s10}$ :

	$\displaystyle U_{n2}=u_{n}$	$\displaystyle+\varphi_{1}(c_{2}hA)c_{2}hF(t_{n},u_{n}),$
	$\displaystyle U_{n\ell}=u_{n}$	$\displaystyle+\varphi_{1}(c_{\ell}hA)c_{\ell}hF(t_{n},u_{n})+\varphi_{2}(c_{\ell}hA)\tfrac{c^{2}_{\ell}}{c_{2}}hD_{n2},\hskip 93.89418pt\ell=3,4$
	$\displaystyle U_{nm}=u_{n}$	$\displaystyle+\varphi_{1}(c_{m}hA)c_{m}hF(t_{n},u_{n})+\varphi_{2}(c_{m}hA)c^{2}_{m}h\big{(}\tfrac{c_{4}}{c_{3}(c_{4}-c_{3})}D_{n3}+\tfrac{c_{3}}{c_{4}(c_{3}-c_{4})}D_{n4}\big{)}$
		$\displaystyle+\varphi_{3}(c_{m}hA)c^{3}_{m}h\big{(}\tfrac{2}{c_{3}(c_{3}-c_{4})}D_{n3}-\tfrac{2}{c_{4}(c_{3}-c_{4})}D_{n4}\big{)},\quad\hskip 71.13188ptm=5,6,7$
	$\displaystyle U_{nq}=u_{n}$	$\displaystyle+\varphi_{1}(c_{q}hA)c_{q}hF(t_{n},u_{n})+\varphi_{2}(c_{q}hA)c^{2}_{q}h\big{(}\alpha_{5}D_{n5}+\alpha_{6}D_{n6}+\alpha_{7}D_{n7}\big{)}$
		$\displaystyle+\varphi_{3}(c_{q}hA)c^{3}_{q}h\big{(}\beta_{5}D_{n5}-\beta_{6}D_{n6}-\beta_{7}D_{n7}\big{)}$
		$\displaystyle+\varphi_{4}(c_{q}hA)c^{4}_{q}h\big{(}\gamma_{5}D_{n5}+\gamma_{6}D_{n6}+\gamma_{7}D_{n7}\big{)},\hskip 96.73918ptq=8,9,10$
	$\displaystyle u_{n+1}=u_{n}$	$\displaystyle+\varphi_{1}(hA)hF(t_{n},u_{n})+\varphi_{2}(hA)h\big{(}\alpha_{8}D_{n8}+\alpha_{9}D_{n9}+\alpha_{10}D_{n10}\big{)}$
		$\displaystyle-\varphi_{3}(hA)h\big{(}\beta_{8}D_{n8}+\beta_{9}D_{n9}+\beta_{10}D_{n10}\big{)}+\varphi_{4}(hA)h\big{(}\gamma_{8}D_{n8}+\gamma_{9}D_{n9}+\gamma_{10}D_{n10}\big{)},$

where

\alpha_{i}=\dfrac{c_{k}c_{l}}{c_{i}(c_{i}-c_{k})(c_{i}-c_{l})},\quad\beta_{i}=\dfrac{2(c_{k}+c_{l})}{c_{i}(c_{i}-c_{k})(c_{i}-c_{l})},\quad\gamma_{i}=\dfrac{6}{c_{i}(c_{i}-c_{k})(c_{i}-c_{l})}

(46)

with $i\in\{5,6,7\}$ for $\ k,l\in\{5,6,7\}$ , and $i\in\{8,9,10\}$ for $\ k,l\in\{8,9,10\}$ (note that $i,k,l$ are distinct indices and that $c_{i},c_{k},c_{l}$ are distinct (positive) nodes).
For our numerical experiments, we choose $c_{2}=c_{3}=c_{5}=\tfrac{1}{2}$ , $c_{4}=c_{6}=\tfrac{1}{3}$ , $c_{7}=\tfrac{1}{4}$ , $c_{8}=\tfrac{3}{10}$ , $c_{9}=\tfrac{3}{4}$ , and $c_{10}=1$ (satisfying (34)).

Remark 4

(A comparison with $\mathtt{expRK5s8}$ ). Although the new fifth-order method $\mathtt{expRK5s10}$ has 10 stages (compared to 8 stages of $\mathtt{expRK5s8}$ displayed in Section 2), its special structure offers much more efficient for implementation. In particular, all $U_{ni}$ in this scheme involve only one linear combination of $\varphi_{k}(c_{i}hA)v_{k}$ which can be computed by one evaluation for each; and more importantly, due to the same format of multiple stages within each of the three groups $\{U_{n3},U_{n4}\}$ , $\{U_{n5},U_{n6},U_{n7}\}$ , and $\{U_{n8},U_{n9},U_{n10}\}$ (same linear combination with different inputs $c_{i}$ ), they can be computed simultaneously or implemented in parallel (see Section 5). This makes $\mathtt{expRK5s10}$ to behave like a 5-stage method only, thereby requiring only 5 sequential evaluations in each step. Moreover, while $\mathtt{expRK5s8}$ requires weakening 9 out of 16 order conditions of Table 1, $\mathtt{expRK5s10}$ requires only one condition (number 8) held in the weakened form. Note that by following the similar way of deriving $\mathtt{expRK5s10}$ , we can derive a scheme that satisfies all the stiff order conditions in Table 1 in the strong sense with $s=11$ . Such a scheme, however, still behaves like a 5-stage method. Therefore, we do not discuss further on this case.

b) The case $s=9$ (which does not work): Clearly, in this case we have less degree of freedoms than the case $s=10$ when solving the order conditions in Table 1. Nonetheless, one can still proceed in a similar way as done for $s=10$ . Again, it strongly suggests $b_{2}=b_{3}=b_{4}=b_{5}=b_{6}=0$ (which solves for $b_{7},b_{8},b_{9}\neq 0$ from conditions 1, 2, 4) and

\psi_{2,j}=\psi_{3,j}=\psi_{4,j}=0\ (j=7,8,9)

(47)

in order to satisfy conditions 1, 2, 3, 4, 5, 7, 9, 13, 15, 16 in the strong form. With this, conditions 6 and 10 now become

\sum_{j=2}^{6}(b_{7}Ja_{7j}+b_{8}Ja_{8j}+b_{9}Ja_{9j})J\psi_{m,j}=0\ \ (m=2,3).

(48)

Again, due to the fact that $\psi_{2,2},\psi_{3,2}\neq 0$ and either $\psi_{2,3}$ or $\psi_{3,3}$ must be nonzero, one needs to enforce $a_{7j}=a_{8j}=a_{9j}=0\ (j=2,3)$ in (48). Using this to solve (47) for $j=7$ ( $\psi_{2,7}=\psi_{3,7}=\psi_{4,7}=0$ ) gives a unique solution (with $c_{4},c_{5},c_{6}>0$ and are distinct) for $a_{74},a_{75},a_{76}\neq 0$ , which then determines $U_{n7}$ . Next, one can solve (47) for $j=8,9$ to obtain $U_{n8},U_{n9}$ that are independent of $U_{n7}$ , as well as are independent of each other, by requiring the three free parameters $a_{87}=a_{97}=a_{98}=0$ . As a result, one gets $a_{7j},a_{8j},a_{9j}\neq 0\ (j=5,6)$ . This immediately suggests $\psi_{2,j}=\psi_{3,j}=0\ (j=4,5,6)$ to completely fulfill (48) with arbitrary square matrix $J$ . With all of these in place, conditions 12 and 14 are automatically fulfilled, and condition 11 is now reduced to

\sum_{i=7}^{9}b_{i}J\sum_{j=4}^{6}a_{ij}J\big{(}a_{j2}J\psi_{2,2}+a_{j3}J\psi_{2,3}\big{)}=0.

(49)

Clearly, since $b_{7},b_{8},b_{9}\neq 0$ , $a_{7j},a_{8j},a_{9j}\neq 0\ (j=4,5,6)$ , and $\psi_{2,2}\neq 0$ , (49) can be satisfied in the strong sense only if we have one of the following conditions: $a_{j2}=a_{j3}=0$ or $a_{j2}=\psi_{2,3}=0,\ (j=4,5,6)$ . Unfortunately, either of these requirements is in contradiction with $\psi_{2,j}=\psi_{3,j}=0\ (j=4,5,6)$ which is needed for conditions 6 and 10 mentioned above. For example, solving $\psi_{2,4}=\psi_{3,4}=0$ results in $a_{42},a_{43}\neq 0$ .

5 Details implementation of fourth- and fifth-order schemes

In this section, we present details implementation of the old and new fourth- and fifth-order expRK schemes ( $\mathtt{expRK4s5}$ , $\mathtt{expRK5s8}$ , $\mathtt{expRK4s6}$ , $\mathtt{expRK5s10}$ ) mentioned above.

As mentioned in Section 2.1, we will use the MATLAB routine phipm_simul_iom (described in details in Luan18 ) to implement expRK methods. In particular, given the following inputs: an array of scaling factors $\mathtt{t}=[\rho_{1},\cdots,\rho_{r}]$ with $0<\rho_{1}<\rho_{2}<\cdots<\rho_{r}\leq 1$ ( $\mathtt{t}$ could be a positive scalar), an $n$ -by- $n$ matrix $M$ , and a set of column vectors $\mathtt{V}=[V_{0},\ldots,v_{q}]$ (each $v_{i}$ is an $n$ -by- $1$ vector), a tolerance $\mathtt{tol}$ , an initial value $m$ for the dimension of the Krylov subspace, and an incomplete orthogonalization length of $\mathtt{iom}$ , a call to this function

\mathtt{phipm\_simul\_iom(t,M,V,tol,m,iom)}

(50)

simultaneously computes the following $r$ linear combinations

L_{\rho_{i},\mathtt{V}}=\varphi_{0}(\rho_{i}M)v_{0}+\varphi_{1}(\rho_{i}M)\rho_{i}v_{1}+\varphi_{2}(\rho_{i}M)\rho_{i}^{2}v_{2}+\cdots+\varphi_{q}(\rho_{i}M)\rho_{i}^{q}v_{q},\ 1\leq i\leq r.

(51)

Note that, by setting $V_{j}=\rho_{i}^{j}v_{j}$ ( $j=0,\cdots,q$ ), (51) becomes (6). In other words, all the linear combinations in (6) ( if $V_{j}$ are given instead of $v_{j}$ ) can be then computed at the same time with one call (50) by using scaled vectors $v_{j}=V_{j}/\rho_{i}^{j}$ for the input $\mathtt{V}$ .

In the following, we set

M=hA,\quad\rho_{i}=c_{i},\quad v=hF(t_{n},u_{n}),\quad d_{i}=hD_{ni}.

(52)

to simplify notations in presenting details of implementation of the fourth- and fifth-order methods mentioned above. When calling (50), we use $\mathtt{tol}=10^{-12}$ , $\mathtt{m}=1$ (default value), and $\mathtt{imo}=2$ (as in Luan18 ).

Implementation of $\mathtt{expRK4s5}$ ( $c_{2}=c_{3}=c_{5}=\tfrac{1}{2},c_{4}=1$ ): As discussed in Remark 1, $\mathtt{expRK4s5}$ requires a sequential implementation of the following 6 different linear combinations of the form (51), corresponding to 6 calls to phipm_simul_iom:

(i)

Evaluate $L_{c_{2},\mathtt{V}}$ with $\mathtt{t}=c_{2},\mathtt{V}=[0,v]$ to get $U_{n2}=u_{n}+L_{c_{2},\mathtt{V}}$ .
(ii)

Evaluate $L_{c_{3},\mathtt{V}}$ with $\mathtt{t}=c_{3},\mathtt{V}=[0,v,d_{2}/c^{2}_{3}]$ to get $U_{n3}=u_{n}+L_{c_{3},\mathtt{V}}$ .
(iii)

Evaluate $L_{c_{4},\mathtt{V}}$ with $\mathtt{t}=c_{4},\mathtt{V}=[0,v,d_{2}+d_{3}]$ to get $U_{n4}=u_{n}+L_{c_{4},\mathtt{V}}$ .
(iv)

Evaluate $L_{c_{5},\mathtt{V_{1}}}$ with $\mathtt{t}=c_{5},\mathtt{V}_{1}=[0,v,2d_{2}+2d_{3}-d_{4},(-d_{2}-d_{3}+d_{4})/c^{2}_{5}]$ and
(v)

Evaluate $L_{c_{4},\mathtt{V_{2}}}$ with $\mathtt{t}=c_{4},\mathtt{V}_{2}=[0,0,(d_{2}+d_{3}-d_{4})/4,(-d_{2}-d_{3}+d_{4})]$
to get $U_{n5}=u_{n}+L_{c_{5},\mathtt{V}_{1}}+L_{c_{4},\mathtt{V_{2}}}$ .
(vi)

Evaluate $L_{1,\mathtt{V}}$ with $\mathtt{t}=1,\mathtt{V}=[0,v,-d_{4}+5d_{5},4d_{4}-8d_{5}]$ to get $u_{n+1}=u_{n}+L_{1,\mathtt{V}}$ .

Since $d_{i}=hD_{ni}$ which depends on $U_{ni}$ , these are the 6 (sequential) evaluations.

Implementation of $\mathtt{expRK4s6}$ ( $c_{2}=c_{3}=\frac{1}{2},c_{4}=c_{6}=\frac{1}{3},c_{5}=\frac{5}{6}$ ): As discussed in Remark 3, $\mathtt{expRK4s6}$ can be implemented like a 4-stage method by evaluating the following 4 sequential evaluations, corresponding to 4 calls to phipm_simul_iom:

(i)

Evaluate $L_{c_{2},\mathtt{V}}$ with $\mathtt{t}=c_{2},\mathtt{V}=[0,v]$ to get $U_{n2}=u_{n}+L_{c_{2},\mathtt{V}}$ .
(ii)

Evaluate $L_{c_{4},\mathtt{V}}$ and $L_{c_{3},\mathtt{V}}$ simultaneously using $\mathtt{t}=[c_{4},c_{3}],\mathtt{V}=[0,v,d_{2}/c_{2}]$ to get both $U_{n3}=u_{n}+L_{c_{3},\mathtt{V}}$ and $U_{n4}=u_{n}+L_{c_{4},\mathtt{V}}$ .
(iii)

Evaluate $L_{c_{5},\mathtt{V}}$ and $L_{c_{6},\mathtt{V}}$ simultaneously with $\mathtt{t}=[c_{6},c_{5}],\\ \mathtt{V}=[0,v,\tfrac{-c_{4}}{(c_{3}-c_{4})c_{3}}d_{3}+\tfrac{c_{3}}{(c_{3}-c_{4})c_{4}}d_{4},\tfrac{1}{(c_{3}-c_{4})c_{3}}d_{3}-\tfrac{1}{(c_{3}-c_{4})c_{4}}d_{4}]$ to get both $U_{n5}=u_{n}+L_{c_{5},\mathtt{V}}$ and $U_{n6}=u_{n}+L_{c_{6},\mathtt{V}}$ .
(iv)

Evaluate $L_{1,\mathtt{V}}$ with $\mathtt{t}=1,\mathtt{V}=[0,v,\tfrac{1}{c_{5}-c_{6}}(\tfrac{-c_{6}}{c_{5}}d_{5}+\tfrac{c_{5}}{c_{6}}d_{6}),\tfrac{2}{c_{5}-c_{6}}(\tfrac{1}{c_{5}}d_{5}-\tfrac{1}{c_{6}}d_{6})]$ to get $u_{n+1}=u_{n}+L_{1,\mathtt{V}}$ .

Implementation of $\mathtt{expRK5s8}$ ( $c_{2}=c_{3}=c_{5}=\frac{1}{2},c_{4}=\frac{1}{4},c_{6}=\frac{1}{5},c_{7}=\frac{2}{3},c_{8}=1$ ): As discussed in Remark 1, $\mathtt{expRK5s8}$ requires a sequential implementation of 11 different linear combinations of the form (51), corresponding to the following 11 calls to phipm_simul_iom:

(i)

Evaluate $L_{c_{2},\mathtt{V}}$ with $\mathtt{t}=c_{2},\mathtt{V}=[0,v]$ to get $U_{n2}=u_{n}+L_{c_{2},\mathtt{V}}$ .
(ii)

Evaluate $L_{c_{3},\mathtt{V}}$ with $\mathtt{t}=c_{3},\mathtt{V}=[0,v,d_{2}/c^{2}_{3}]$ to get $U_{n3}=u_{n}+L_{c_{3},\mathtt{V}}$ .
(iii)

Evaluate $L_{c_{4},\mathtt{V}}$ with $\mathtt{t}=c_{4},\mathtt{V}=[0,v,d_{3}/c^{2}_{4}]$ to get $U_{n4}=u_{n}+L_{c_{4},\mathtt{V}}$ .
(iv)

Evaluate $L_{c_{5},\mathtt{V}}$ with $\mathtt{t}=c_{5},\mathtt{V}=[0,v,(-d_{3}+4d_{4})/c^{2}_{5},(2d_{3}-4d_{4})/c^{3}_{5}]$ to get $U_{n5}=u_{n}+L_{c_{5},\mathtt{V}}$ .
(v)

Evaluate $L_{c_{6},\mathtt{V}}$ with $\mathtt{t}=c_{6},\mathtt{V}=[0,v,(8d_{4}-2d_{5})/c^{2}_{6},(-32d_{4}+16d_{5})/c^{3}_{6}]$ to get $U_{n6}=u_{n}+L_{c_{6},\mathtt{V}}$ .
(vi)

Evaluate $L_{c_{7},\mathtt{V_{1}}}$ with $\mathtt{t}=c_{7},\mathtt{V}_{1}=[0,v,(\tfrac{-16}{27}d_{5}+\tfrac{100}{27}d_{6})/c^{2}_{7},(\tfrac{320}{81}d_{5}-\tfrac{800}{81}d_{n6})/c^{3}_{7}]$ and
(vii)

Evaluate $L_{c_{6},\mathtt{V_{2}}}$ with $\mathtt{t}=c_{6},\mathtt{V}_{2}=[0,0,(\tfrac{-20}{81}d_{4}+\tfrac{5}{243}d_{5}+\tfrac{125}{486}d_{6})/c^{2}_{6},(\tfrac{16}{81}d_{4}-\tfrac{4}{243}d_{5}-\tfrac{50}{243}d_{6})/c^{3}_{6}]$ to get $U_{n7}=u_{n}+L_{c_{7},\mathtt{V}_{1}}+L_{c_{6},\mathtt{V_{2}}}$ .
(viii)

Evaluate $L_{c_{8},\mathtt{V_{1}}}$ with $\mathtt{t}=c_{8},\mathtt{V}_{1}=[0,v,(\tfrac{-16}{3}d_{5}+\tfrac{250}{21}d_{6}+\tfrac{27}{14}d_{7})/c^{2}_{8},(\tfrac{208}{3}d_{5}-\tfrac{250}{3}d_{6}-27d_{7})/c^{3}_{8},(-240d_{5}+\tfrac{1500}{7}d_{6}+\tfrac{810}{7}d_{7})/c^{4}_{8}]$ and
(ix)

Evaluate $L_{c_{6},\mathtt{V_{2}}}$ with $\mathtt{t}=c_{6},\\ \mathtt{V}_{2}=[0,0,(\tfrac{-4}{7}d_{5}+\tfrac{25}{49}d_{6}+\tfrac{27}{98}d_{7})/c^{2}_{6},(\tfrac{8}{5}d_{5}-\tfrac{10}{7}d_{6}-\tfrac{27}{35}d_{7})/c^{3}_{6},(\tfrac{-48}{35}d_{5}+\tfrac{60}{49}d_{6}+\tfrac{162}{245}d_{7})/c^{4}_{6}]$ and
(x)

Evaluate $L_{c_{7},\mathtt{V_{3}}}$ with $\mathtt{t}=c_{7},\mathtt{V}_{3}=[0,0,(\tfrac{-288}{35}d_{5}+\tfrac{360}{49}d_{6}+\tfrac{972}{245}d_{7})/c^{2}_{7},(\tfrac{384}{5}d_{5}-\tfrac{480}{7}d_{6}-\tfrac{1296}{35}d_{7})/c^{3}_{7},(\tfrac{-1536}{7}d_{5}+\tfrac{9600}{49}d_{6}+\tfrac{5184}{49}d_{7})/c^{4}_{7}]$
to get $U_{n8}=u_{n}+L_{c_{8},\mathtt{V}_{1}}+L_{c_{6},\mathtt{V}_{2}}+L_{c_{7},\mathtt{V_{3}}}$ .
(xi)

Evaluate $L_{1,\mathtt{V}}$ with $\mathtt{t}=1,\mathtt{V}=[0,v,\tfrac{125}{14}d_{6}-\tfrac{27}{14}d_{7}+\tfrac{1}{2}d_{8},\tfrac{-625}{14}d_{6}+\tfrac{162}{7}d_{7}-\tfrac{13}{2}d_{8},\tfrac{1125}{14}d_{6}-\tfrac{405}{7}d_{7}+\tfrac{45}{2}d_{8}]$ to get $u_{n+1}=u_{n}+L_{1,\mathtt{V}}$ .

Implementation of $\mathtt{expRK5s10}$ ( $c_{2}=c_{3}=c_{5}=\tfrac{1}{2}$ , $c_{4}=c_{6}=\tfrac{1}{3}$ , $c_{7}=\tfrac{1}{4}$ , $c_{8}=\tfrac{3}{10}$ , $c_{9}=\tfrac{3}{4}$ , and $c_{10}=1$ ): As discussed in Remark 4, $\mathtt{expRK5s10}$ can be implemented like a 5-stage method by evaluating the following 5 sequential evaluations, corresponding to 5 calls to phipm_simul_iom:

(i)

Evaluate $L_{c_{2},\mathtt{V}}$ with $\mathtt{t}=c_{2},\mathtt{V}=[0,v]$ to get $U_{n2}=u_{n}+L_{c_{2},\mathtt{V}}$ .
(ii)

Evaluate $L_{c_{4},\mathtt{V}}$ and $L_{c_{3},\mathtt{V}}$ simultaneously using $\mathtt{t}=[c_{4},c_{3}],\mathtt{V}=[0,v,d_{2}/c_{2}]$ to get both $U_{n3}=u_{n}+L_{c_{3},\mathtt{V}}$ and $U_{n4}=u_{n}+L_{c_{4},\mathtt{V}}$ .
(iii)

Evaluate $L_{c_{5},\mathtt{V}}$ , $L_{c_{6},\mathtt{V}}$ , and $L_{c_{7},\mathtt{V}}$ simultaneously using $\mathtt{t}=[c_{7},c_{6},c_{5}]$ ,
$\mathtt{V}=[0,v,\tfrac{c_{4}}{c_{3}(c_{4}-c_{3})}d_{3}+\tfrac{c_{3}}{c_{4}(c_{3}-c_{4})}d_{4},\tfrac{2}{c_{3}(c_{3}-c_{4})}d_{3}-\tfrac{2}{c_{4}(c_{3}-c_{4})}d_{4}]$
to get $U_{n5}=u_{n}+L_{c_{5},\mathtt{V}}$ , $U_{n6}=u_{n}+L_{c_{6},\mathtt{V}}$ , $U_{n7}=u_{n}+L_{c_{7},\mathtt{V}}$ .
(iv)

Evaluate $L_{c_{8},\mathtt{V}}$ , $L_{c_{9},\mathtt{V}}$ , and $L_{c_{10},\mathtt{V}}$ simultaneously using $\mathtt{t}=[c_{9},c_{10},c_{8}]$ ,
$\mathtt{V}=[0,v,\alpha_{5}d_{5}+\alpha_{6}d_{6}+\alpha_{7}d_{7},\beta_{5}d_{5}-\beta_{6}d_{6}-\beta_{7}d_{7},\gamma_{5}d_{5}+\gamma_{6}d_{6}+\gamma_{7}d_{7}]$
to get $U_{n8}=u_{n}+L_{c_{8},\mathtt{V}}$ , $U_{n9}=u_{n}+L_{c_{9},\mathtt{V}}$ , $U_{n10}=u_{n}+L_{c_{10},\mathtt{V}}$ .
(v)

Evaluate $L_{1,\mathtt{V}}$ with $\mathtt{t}=1,\mathtt{V}=[0,v,\alpha_{8}d_{8}+\alpha_{9}d_{9}+\alpha_{10}d_{10},\beta_{8}d_{8}+\beta_{9}d_{9}+\beta_{10}d_{10},\gamma_{8}d_{8}+\gamma_{9}d_{9}+\gamma_{10}d_{10}]$ to get $u_{n+1}=u_{n}+L_{1,\mathtt{V}}$ (coefficients $\alpha_{i},\beta_{i},\gamma_{i}$ are given in (46)).

6 Numerical experiments

In this section, we demonstrate the efficiency of our newly derived fourth- and fifth-order expRK time integration methods ( $\mathtt{expRK4s6}$ , $\mathtt{expRK5s10}$ ). Specifically, we will compare their performance against the existing methods of the same orders ( $\mathtt{expRK4s5}$ , $\mathtt{expRK5s8}$ ) on several examples of stiff PDEs. All the numerical simulations are performed in MATLAB on a single workstation, using an iMac 3.6 GHz Intel Core i7, 32 GB 2400 MHz DDR4.

Example 1

(A one-dimensional semilinear parabolic problem HO05 ): We first verify the order of convergence for the new derived fourth- and fifth-order expRK schemes ( $\mathtt{expRK4s6}$ , $\mathtt{expRK5s10}$ ) by considering the following PDE for $u(x,t)$ , $x\in[0,1],t\in[0,1]$ , and subject to homogeneous Dirichlet boundary conditions,

\frac{\partial u(x,t)}{\partial t}-\frac{\partial^{2}u(x,t)}{\partial x^{2}}=\frac{1}{1+u^{2}(x,t)}+\Phi(x,t),

(53)

whose exact solution is known to be $u(x,t)=x(1-x){\rm e}\hskip 1.0pt^{t}$ for a suitable choice of the source function $\Phi(x,t)$ .

Spatial discretization: For this example, we use standard second order finite differences with $200$ grid points. This leads to a very stiff system of the form (1) (with $\|A\|_{\infty}\approx 1.6\times 10^{5}$ ).

The resulting system is then integrated on the time interval $[0,1]$ using constant step sizes, corresponding to the number of time steps $N=4,8,16,32,64$ . The time integration errors at the final time $t=1$ are measured in the maximum norm.

In Figure 1, we plot orders for all the employed integrators in the left diagram and the total CPU time versus the global errors in the right diagram. The left diagram clearly shows a perfect agreement with our convergence result in Theorem 3.1, meaning that the two new integrators $\mathtt{expRK4s6}$ and $\mathtt{expRK5s10}$ fully achieve orders 4 and 5, respectively. When compared to the old integrators of the same orders $\mathtt{expRK4s5}$ and $\mathtt{expRK5s8}$ , we note that, given the same number of time steps, $\mathtt{expRK4s6}$ is slightly more accurate but is much faster than $\mathtt{expRK4s5}$ (see the right diagram). In a similar manner, $\mathtt{expRK5s10}$ gives almost identical global errors but is also much faster than $\mathtt{expRK5s8}$ . Finally, we observe that, for this example, for a global error that is larger than $10^{-6}$ , the new fourth-order method $\mathtt{expRK4s6}$ is the fastest one, and for more stringent errors, $\mathtt{expRK5s10}$ is the fastest integrator.

Refer to caption — Figure 1: Order plots (left) and total CPU times (right) of $\mathtt{expRK4s5}$ , $\mathtt{expRK4s6}$ , $\mathtt{expRK5s8}$ , and $\mathtt{expRK5s10}$ when applied to (53). The global errors at time $t=1$ are plotted as functions of the number of time steps (left) and the total CPU time in second (right). For comparison, straight lines with slopes 4 and 5 are added.

Example 2

(A nonlinear Schrödinger equation Cazenave1989 ; Berland2005 ): We consider the following one-dimensional nonlinear Schrödinger (NLS) equation with periodic boundary conditions

$\displaystyle\mathtt{i}\frac{\partial\Psi(x,t)}{\partial t}$	$\displaystyle=-\frac{\partial^{2}\Psi(x,t)}{\partial x^{2}}+\big{(}V(x)+\lambda\|\Psi(x,t)\|^{2}\big{)}\Psi(x,t),$	(54)
$\displaystyle\Psi(-\pi,t)$	$\displaystyle=\Psi(\pi,t),\quad t\geq 0$
$\displaystyle\Psi(0,t)$	$\displaystyle=\Psi_{0}(x),\quad x\in[-\pi,\pi]$

where the potential function $V(x)=\dfrac{1}{1+\sin^{2}(x)}$ , the initial condition $\Psi_{0}(x)={\rm e}\hskip 1.0pt^{\sin(2x)}$ , and the constant $\lambda=1$ (see Berland2005 ).

Spatial discretization: For this example, we use a discrete Fourier transform $\mathcal{F}$ with $ND=128$ modes, leading to a mildly stiff system of the form (1) with

	$\displaystyle A$	$\displaystyle=\text{diag}(-\mathtt{i}k^{2}),\ k=-\frac{ND}{2}+1,\cdots,\frac{ND}{2}=-63,\cdots,64$		(55)
	$\displaystyle g(t,u)$	$\displaystyle=-\mathtt{i}\mathcal{F}((V(x)+\lambda\|\mathcal{F}^{-1}(u)\|^{2})\mathcal{F}^{-1}(u).$		(55)

Next, we integrate this system on the time interval $[0,3]$ with constant step sizes, corresponding to the number of time steps $N=64,128$ , $256,512,1024$ . Since the exact solution $\Psi(x,t)$ of (54) is unknown, a reliable reference solution is computed by the stiff solver $\mathtt{ode15s}$ with $ATOL=RTOL=10^{-14}$ . Again, the time integration errors are measured in a discrete maximum norm at the final time $t=3$ .

As seen from the two double-logarithmic diagrams in Figure 2, we plot the accuracy of the four employed integrators ( $\mathtt{expRK4s5}$ , $\mathtt{expRK4s6}$ , $\mathtt{expRK5s8}$ , and $\mathtt{expRK5s10}$ ) as functions of the number of time steps (left) and the total CPU time (right). The left digram clearly indicates that the two new integrators $\mathtt{expRK4s6}$ and $\mathtt{expRK5s10}$ achieve their corresponding expected orders 4 and 5. While $\mathtt{expRK5s10}$ is a little more accurate than $\mathtt{expRK5s8}$ , $\mathtt{expRK4s6}$ is much more accurate than $\mathtt{expRK4s5}$ for a given same number of time steps, meaning that it can take much larger time steps while achieving the same accuracy. Moreover, the right precision digram displays the efficiency plot indicating that both $\mathtt{expRK4s6}$ and $\mathtt{expRK5s10}$ are much faster than their counterparts $\mathtt{expRK4s5}$ and $\mathtt{expRK5s8}$ , respectively. More specifically, a similar story is observed: for lower accuracy requirements, say error $\sim 10^{-7}$ , the new fourth-order method $\mathtt{expRK4s6}$ is the most efficient, whereas for error $\sim 10^{-8}$ or tighter the new fifth-order method $\mathtt{expRK5s10}$ is the most efficient.

Example 3

(A 2D Gray–Scott model Gray1984 ; Berland2007 ): Consider the following two-dimensional reaction-diffusion equation–the Gray–Scott equation model, for $u=u(x,y,t),\ v=v(x,y,t)$ on the square $\Omega=[0,L]^{2}$ , (here, we choose $L=1.5$ ) subject to periodic boundary conditions

	$\displaystyle\frac{\partial u}{\partial t}$	$\displaystyle=d_{u}\Delta u-uv^{2}+\alpha(1-u),$		(56)
	$\displaystyle\frac{\partial v}{\partial t}$	$\displaystyle=d_{v}\Delta v+uv^{2}-(\alpha+\beta)v,$		(56)

where $\Delta$ is the Laplacian operator, the diffusion coefficients $d_{u}=0.02,\ d_{v}=0.01$ , and the bifurcation parameters $\alpha=0.065,\ \beta=0.035$ . The initial conditions are Gaussian pulses

u(x,y,0)=1-e^{-150\big{(}(x-L)^{2}+(y-L)^{2}\big{)}},\ v(x,y,0)=e^{-150\big{(}(x-L)^{2}+2(y-L)^{2}\big{)}}.

Spatial discretization: For this example, we use standard second order finite differences using 150 grid points in each direction with mesh width $\Delta x=\Delta y=L/150$ . This gives a stiff system of the form (1).

The system is then solved on the time interval $[0,2]$ using constant step sizes. In the absence of an analytical solution of (56), a high-accuracy reference solution is computed using the $\mathtt{expRK4s6}$ method with a sufficient small time step. Errors are measured in a discrete maximum norm at the final time $t=2$ .

In Figure 3, using the same number of time steps $N=32,64,128$ , $256,512,1024$ , we again display the order plots of the taken integrators. One can see that $\mathtt{expRK4s6}$ is much more accurate than $\mathtt{expRK4s5}$ and $\mathtt{expRK5s10}$ is slightly more accurate than $\mathtt{expRK5s8}$ .

In Figure 4, we display the efficiency plot for which the time step sizes were chosen for each integrator to obtain about same error thresholds $10^{-i},\ i=5,\cdots,11$ (The corresponding number of time steps for each integrator are displayed in Table 2. As seen, given about the same level of accuracy, the new methods use smaller steps than the old ones of the same order, meaning that they can take larger step sizes). Again, $\mathtt{expRK4s6}$ is much faster than $\mathtt{expRK4s5}$ and it is interesting that this new fourth-order method turns out to be the most efficient (although for error thresholds tighter than $10^{-11}$ the new fifth-order method $\mathtt{expRK5s10}$ seems to become the most efficient).

Method	Error threshold vs. Number of time steps
Method	$10^{-5}$	$10^{-6}$	$10^{-7}$	$10^{-8}$	$10^{-9}$	$10^{-10}$	$10^{-11}$
$\mathtt{expRK4s5}$	18	36	66	121	215	385	685
$\mathtt{expRK4s6}$	10	19	28	46	122	230	420
$\mathtt{expRK5s8}$	7	18	33	57	92	149	238
$\mathtt{expRK5s10}$	8	17	30	51	82	130	208

Table 2: The number of time steps taken to achieve about the same error thresholds

10^{-i},\ i=5,\cdots,11

The numerical results presented on the three examples above clearly confirm the advantage of constructing parallel stages expRK methods based on Theorem 3.1, leading to more efficient and accurate methods $\mathtt{expRK4s6}$ and $\mathtt{expRK5s10}$ .

Acknowledgements.

The author would like to thank Reviewer 1 for the valuable comments and helpful suggestions. He would like also to thank the National Science Foundation, which supported this research under award NSF DMS–2012022.

References

(1) Al-Mohy, A.H., Higham, N.J.: Computing the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput. 33, 488–511 (2011)
(2) Berland, H., Skaflestad, B.: Solving the nonlinear Schrödinger equation using exponential integrators. Tech. rep. (2005)
(3) Berland, H., Skaflestad, B., Wright, W.M.: Expint—a matlab package for exponential integrators. ACM Transactions on Mathematical Software (TOMS) 33(1), 4–es (2007)
(4) Caliari, M., Kandolf, P., Ostermann, A., Rainer, S.: The Leja method revisited: Backward error analysis for the matrix exponential. SIAM J. Sci. Comp. 38(3), A1639–A1661 (2016)
(5) Cazenave, T.: An introduction to nonlinear Schrödinger equations, vol. 22. Universidade Federal do Rio de Janeiro, Centro de Ciências Matemáticas e da … (1989)
(6) Gaudreault, S., Pudykiewicz, J.: An efficient exponential time integration method for the numerical solution of the shallow water equations on the sphere. J. Comput. Phys. 322, 827–848 (2016)
(7) Gray, P., Scott, S.: Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Oscillations and instabilities in the system $A+2B\rightarrow 3B;B\rightarrow C$ . Chemical Engineering Science 39(6), 1087–1097 (1984)
(8) Hochbruck, M., Ostermann, A.: Explicit exponential Runge–Kutta methods for semilinear parabolic problems. SIAM J. Numer. Anal. 43, 1069–1090 (2005)
(9) Hochbruck, M., Ostermann, A.: Exponential integrators. Acta Numerica 19, 209–286 (2010)
(10) Ju, L., Wang, Z.: Exponential time differencing Gauge method for incompressible viscous flows. Communications in Computational Physics 22(2), 517–541 (2017)
(11) Luan, V.T.: High-order exponential integrators. Ph.D. thesis, University of Innsbruck (2014)
(12) Luan, V.T.: Fourth-order two-stage explicit exponential integrators for time-dependent PDEs. Applied Numerical Mathematics 112, 91–103 (2017)
(13) Luan, V.T., Michels, D.: Explicit exponential Rosenbrock methods and their application in visual computing. (Revised) (2020)
(14) Luan, V.T., Ostermann, A.: Exponential B-series: The stiff case. SIAM J. Numer. Anal. 51, 3431–3445 (2013)
(15) Luan, V.T., Ostermann, A.: Explicit exponential Runge–Kutta methods of high order for parabolic problems. J. Comput. Appl. Math. 256, 168–179 (2014)
(16) Luan, V.T., Ostermann, A.: Exponential Rosenbrock methods of order five–construction, analysis and numerical comparisons. J. Comput. Appl. Math. 255, 417–431 (2014)
(17) Luan, V.T., Ostermann, A.: Stiff order conditions for exponential Runge–Kutta methods of order five. In: H.B. et al. (ed.) Modeling, Simulation and Optimization of Complex Processes - HPSC 2012, pp. 133–143. Springer (2014)
(18) Luan, V.T., Ostermann, A.: Parallel exponential Rosenbrock methods. Comput. Math. Appl. 71, 1137–1150 (2016)
(19) Luan, V.T., Pudykiewicz, J.A., Reynolds, D.R.: Further development of efficient and accurate time integration schemes for meteorological models. J. Comput. Phys. 376, 817–837 (2019)
(20) Michels, D.L., Luan, V.T., Tokman, M.: A stiffly accurate integrator for elastodynamic problems. ACM Transactions on Graphics (TOG) 36(4), 116 (2017)
(21) Niesen, J., Wright, W.M.: Algorithm 919: A Krylov subspace algorithm for evaluating the $\varphi$ -functions appearing in exponential integrators. ACM Trans. Math. Soft. (TOMS) 38(3), 22 (2012)
(22) Pieper, K., Sockwell, K.C., Gunzburger, M.: Exponential time differencing for mimetic multilayer ocean models. J. Comput. Phys. 398, 817–837 (2019)

Efficient exponential Runge–Kutta methods of high order: construction and implementation ††thanks: This work has been supported in part by National Science Foundation through award NSF DMS–2012022.