Efficient exponential Runge--Kutta methods of high order: construction and implementation
aa r X i v : . [ m a t h . NA ] S e p BIT manuscript No. (will be inserted by the editor)
Efficient exponential Runge–Kutta methods of high order:construction and implementation
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 timeintegration of stiff semilinear parabolic PDEs. The current construction of stiffly accurateexponential Runge–Kutta methods, however, relies on a convergence result that requiresweakening many of the order conditions, resulting in schemes whose stages must be im-plemented 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–Kuttamethods, which, in contrast to the existing methods, have multiple stages that are indepen-dent of one another and share the same format, thereby allowing them to be implementedin parallel or simultaneously, and making the methods to behave like using with much lessstages. Moreover, all of their stages involve only one linear combination of the product of ϕ -functions (using the same argument) with vectors. Overall, these features make these newmethods to be much more efficient to implement when compared to the existing methods ofthe same orders. Numerical experiments on a one-dimensional semilinear parabolic prob-lem, a nonlinear Schrdinger equation, and a two-dimensional Gray–Scott model are givento confirm the accuracy and efficiency of the two newly constructed methods. Keywords
Exponential Runge–Kutta methods · Exponential integrators · Stiff PDEs · Efficient implementation
Mathematics Subject Classification (2000)
MSC 65L04 · MSC 65M06 · MSC 65N12
In this paper, we are concerned with the construction and implementation of new efficientexponential Runge–Kutta integrators for solving stiff parabolic PDEs. These PDEs, upontheir spatial discretizations, can be cast in the form of semilinear problems u ′ ( t ) = Au ( t ) + g ( t , u ( t )) = F ( t , u ( t )) , u ( t ) = u , (1.1) This work has been supported in part by National Science Foundation through award NSF DMS–2012022.Vu Thai LuanDepartment of Mathematics and Statistics, Mississippi State University410 Allen Hall, 175 President’s Circle, Mississippi State, MS, 39762E-mail: [email protected] Vu Thai Luan where the linear part Au usually causes stiffness. The nonlinearity g ( t , u ) is assumed tosatisfy a local Lipschitz condition in a strip along the exact solution.Exponential Runge–Kutta methods are a popular class of exponential integrators [9],which have shown a great promise as an alternative to standard time integration solvers forstiff systems and applications in recent years, see e.g. [8, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 22]. The main idea behind these methods is to solve the linear portion of (1.1)exactly and integrate the remaining nonlinear portion explicitly based on a representation ofthe exact solution using the variation-of-constants formula.A s -stage explicit exponential Runge–Kutta (expRK) method [8] applied to (1.1) can bereformulated (see [15, 17]) as U ni = u n + c i h ϕ ( c i hA ) F ( t n , u n ) + h i − ∑ j = a i j ( hA ) D n j , ≤ i ≤ s , (1.2a) u n + = u n + h ϕ ( hA ) F ( t n , u n ) + h s ∑ i = b i ( hA ) D ni , (1.2b)where D ni = g ( t n + c i h , U ni ) − g ( t n , u n ) , ≤ i ≤ s . (1.3)Here, U ni denote the internal stages that approximate u ( t n + c i h ) using the time step size h = t n + − t n > c i . By construction, the coefficients a i j ( z ) and b i ( z ) are usuallylinear combinations of the entire functions ϕ k ( z ) = Z e ( − θ ) z θ k − ( k − ) ! d θ , k ≥ ϕ k ( c i z ) .A common approach that has been used to determine the unknown matrix functions a i j ( hA ) and b i ( hA ) is to expand them as a i j ( hA ) = ∑ k ≥ α ( k ) i j ( hA ) k , b i ( hA ) = ∑ k ≥ β ( k ) i ( hA ) k (e.g. using classical Taylor series expansions) to obtain order conditions. Clearly, the bound-edness of the remainder terms of these expansions (and thus the error terms) are dependentof k A k . Due to stability reasons, such resulting methods might not be suitable for integrat-ing stiff PDEs, which A typically has a large norm or is even unbounded operator. Thesemethods are thus usually referred as classical (non-stiffly accurate) expRK methods. Unlikethis approach, in a seminal contribution [8], Hochbruck and Ostermann derived a new errorexpansion with the remainder terms that are bounded independently of the stiffness (i.e. notinvolving the powers of A ), leading to stiff order conditions, which give rise to the construc-tion of stiffly accurate expRK methods of orders up to four. Following this, in [14] Luanand Ostermann developed a systematic theory of deriving stiff order conditions for expRKmethods of arbitrary order, thereby allowing the construction of a fifth-order method in [15].In view of the existing stiffly accurate expRK methods in the literature, we observe thatthey were derived based on a convergence result that requires weakening many of the stifforder conditions (in order to minimize the number of required stages s and matrix func-tions 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 imple-mented by computing each of these stages sequentially. Also, the very last stages usuallyinvolve several different linear combinations of ϕ k ( c i hA ) -functions (using different nodes c i in their arguments) acting on different sets of vectors. This would introduce additionalcomputational effort for these stages. For more details, we refer to Sections 2 and 5. fficient exponential Runge–Kutta methods of high order 3 Motivated by the observations above, in this work we show a stronger convergence re-sult for expRK methods up to order five which requires weakening only one order con-dition (thereby could improve the stability and accuracy) and offers more degree of free-doms in solving order conditions. Using this result and inspired by our recent algorithm, phipm simul iom , proposed in [19] (which allows one to simultaneously compute multi-ple linear combinations of ϕ - functions acting on a same set of vectors), we construct newmethods of orders 4 and 5 which involve only one linear combination of ϕ - functions foreach stage and have multiple internal stages U ni that are independent of one another, therebyallowing them to be computed in parallel. Furthermore, one can derive these independentstages in a way that they share the same form of linear combination of ϕ k ( c i hA ) - functionsacting on the same set of vectors, allowing them to be implemented simultaneously (byone evaluation). While these independent states can be computed in parallel (as mentionedabove) by any algorithm which approximates the action of (the linear combination of) ϕ -functions, we note that the possibility to compute them simultaneously is a new feature thatcan be used with our algorithm phipm simul iom (other algorithms, e.g., that do not requirethe construction of Krylov subspaces, might not support computing these stages simultane-ously). Overall, this makes the new methods to behave like methods using much less numberof stages (even when compared to the existing methods of the same orders), meaning thatthey require much less number of evaluations for linear combinations of ϕ - functions, andare thus more efficient.The paper is organized as follows. In Section 2, we describe our motivation, proposenew ideas, and review the existing expRK methods in the literature with respect to theseideas. 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 toconstruct more efficient methods in Section 4. In particular, we are able to derive two newfamilies of fourth- and fifth- order stiffly accurate expRK methods called expRK4s6 (4th-order 6-stage but requires 4 independent stage evaluation only) and expRK5s10 (5th-order10-stage but requires 5 independent stage evaluation only), respectively. In Section 5, wepresent details implementation of these two new methods, as well as the existing stifflyaccurate expRK schemes of the same orders (for comparison). In the last section, numericalexamples including one and two-dimensional stiff PDEs are presented to demonstrate theaccuracy and efficiency of the two newly constructed expRK integrators. In this section, we start our motivation by taking a closer look at an efficient way for imple-menting expRK methods (1.2). Then, we propose some ideas to derive more efficient meth-ods with respect to this efficient implementation along with reviewing the current methods.2.1 An efficient way of implementationClearly, each stage ( U ni or u n + ) of (1.2) requires computing matrix functions of the form ϕ k ( c i hA ) v k (0 < c i ≤ v k is some vector (could be F ( t n , u n ) , D ni or a linear combi-nation of these). Thanks to recent developments [1, 4, 6, 21], one can efficiently compute alinear combination of ϕ -functions acting on a set of input vectors V , . . ., V q ϕ ( M ) V + ϕ ( M ) V + ϕ ( M ) V + · · · + ϕ q ( M ) V q , (2.1) Vu Thai Luan where M is some square matrix. This is crucial when implementing exponential integrators.Very recently, in [19], we were able to improve the implementations presented in [6, 21], re-sulting in the routine phipm simul iom . The underlying method in this algorithm is the useof an adaptive time-stepping technique combined with Krylov subspace methods, which al-lows us to simultaneously compute multiple linear combinations of type (2.1) using differentscaling factors ρ , · · · , ρ r of M , i.e., ϕ ( ρ M ) V + ϕ ( ρ M ) V + ϕ ( ρ M ) V + · · · + ϕ N ( ρ M ) V q , ... ϕ ( ρ r M ) V + ϕ ( ρ r M ) V + ϕ ( ρ r M ) V + · · · + ϕ N ( ρ r M ) V q . (2.2)Now taking M = hA and considering ρ k (1 ≤ k ≤ r ) as nodes c i used in expRK methodsimmediately suggests that one can compute the following ( s −
1) linear linear combinations ϕ ( c i hA ) V + ϕ ( c i hA ) V + . . . + ϕ N ( c i hA ) V q , ≤ i ≤ s (2.3)simultaneously by using only one evaluation (i.e., one call to phipm simul iom ). Note thatthis requires the use of a same set of vectors [ V , . . ., V q ] for all the linear combinationsin (2.3).Motivated by this, we see that if a s -stage expRK scheme (1.2) is constructed in such away that each internal stage U ni has the form U ni = u n + ϕ ( c i hA ) V i + ϕ ( c i hA ) V i + . . . + ϕ N ( c i hA ) V qi , (2.4)which includes only one linear combination of ϕ - functions using exactly node c i as an argu-ment in all ϕ k functions, then the scheme will contain a total of s such linear combinations( s − U ni and 1 for u n + as (1.2b) can be always written in the form of (2.4) with c i = s evaluations only. Furthermore, since the sets of vectors [ V i , V i , · · · , V qi ] in (2.4) are usually different for each U ni , (2.3) also suggests that the efficiency will be sig-nificantly increased if one could build such stages (or a group of) U ni of the form (2.4) thatshare the same format (i.e., having the same set of acting vectors V i ≡ V , . . ., V qi ≡ V q ) orthat are independent of one another. As this allows to compute such stages simultaneouslyby one evaluation or to implement them in parallel similarly to our construction of paral-lel exponential Rosenbrock methods [18]), it certainly reduces the total number of requiredevaluations and thus speedups the computing time.With respect to these observations, we now review the existing expRK schemes in theliterature. Since our focus is on stiff problems, we will discuss only on stiffly accurate expRKmethods, meaning that they satisfy the stiff order conditions (see Section 3 below).2.2 Existing schemes and remarksIn [8], expRK methods of orders up to four have been derived. For later reference, we namethe second-order, the third-order, and the fourth-order methods in that work as expRK2s2 , expRK3s3 , and expRK4s5 , respectively. In [15], we have constructed an expRK method oforder five called expRK5s8 . To discuss all of these schemes in terms of the implementation,we rewrite their internal stages U ni and u n + as linear combinations of ϕ - functions like (2.4)and display them as follows (Note that since the first-order method, the exponential Eulerscheme u n + = u n + ϕ ( hA ) hF ( t n , u n ) , has no internal stage, we do not consider it here). fficient exponential Runge–Kutta methods of high order 5 expRK2s2 : U n = u n + ϕ ( c hA ) c hF ( t n , u n ) , u n + = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) c hD n . (2.5) expRK3s3 (a representative with c = ): U n = u n + ϕ ( c hA ) c hF ( t n , u n ) , U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) c hD n , u n + = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) hD n . (2.6) expRK4s5 (the only existing fourth-order stiffly accurate expRK method constructed byHochbruck and Ostermann [8]): U n = u n + ϕ ( hA ) hF ( t n , u n ) , U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) hD n , U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( D n + D n ) , U n = u n + [ ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( D n + D n − D n )+ ϕ ( hA ) h ( − D n − D n + D n )] + [ ϕ ( hA ) h ( D n + D n − D n )+ ϕ ( hA ) h ( − D n − D n + D n )] , u n + = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( − D n + D n ) + ϕ ( hA ) h ( D n − D n ) . (2.7) expRK5s8 (the only existing fifth-order stiffly accurate expRK method constructed by Luanand Ostermann [15]): U n = u n + ϕ ( hA ) hF ( t n , u n ) , U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) hD n , U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) hD n , U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( − D n + D n ) + ϕ ( hA ) h ( D n − D n ) U n = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( D n − D n ) + ϕ ( hA ) h ( − D n + D n ) , U n = u n + [ ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( − D n + D n ) + ϕ ( hA ) h ( D n − D n )]+ [ ϕ ( hA ) h ( − D n + D n + D n ) + ϕ ( hA ) h ( D n − D n − D n )] , U n = u n + (cid:2) ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( − D n + D n + D n )+ ϕ ( hA ) h ( D n − D n − D n ) + ϕ ( hA ) h ( − D n + D n + D n ) (cid:3) + (cid:2) ϕ ( hA ) h ( − D n + D n + D n ) + ϕ ( hA ) h ( D n − D n − D n )+ ϕ ( hA ) h ( − D n + D n + D n ) (cid:3) + (cid:2) ϕ ( hA ) h ( − D n + D n + D n ) + ϕ ( hA ) h ( D n − D n − D n )+ ϕ ( hA ) h ( − D n + D n + D n ) (cid:3) , u n + = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h ( D n − D n + D n )+ ϕ ( hA ) h ( − D n + D n − D n ) + ϕ ( hA ) h ( D n − D n + D n ) . Vu Thai Luan
Remark 2.1
In view of the structures of these schemes, one can see that only the second- andthird-oder schemes ( expRK2s2 , expRK3s3 ) have all U ni in the form (2.4). While expRK2s2 requires one internal stage U n , expRK3s3 needs two internal stages with U n depends on U n , making these stages cannot be computed simultaneously. As for expRK4s5 , to the bestof our knowledge, this 5-stage scheme is the only existing fourth-order stiffly accurate ex-pRK method. As seen, among its internal stages the three internal stages U n , U n , and U n are of the form (2.4) but again their corresponding sets of vectors [ V ki ] are not thesame ( [ V k ] = [ hF ( t n , u n )] , [ V k ] = [ hF ( t n , u n ) , hD n ] , [ V k ] = [ hF ( t n , u n ) , h ( D n + D n )] ),and because of (1.3), they are not independent of one another ( U n and U n depend on theirpreceding stages). Therefore one needs 3 sequential evaluations for computing these threestages. Also, we note that the last internal stage U n depends on all of its preceding stagesand involves two different linear combinations of ϕ k - functions with different scaling factors c = and c =
1, namely, ∑ k ϕ k ( hA ) V k and ∑ k ϕ k ( hA ) W k (grouped in two different brack-ets [ ] ), which has to be implemented by 2 separate evaluations. The final stage u n + dependson U n and U n . As a result, this scheme must be implemented in a sequential way, whichrequires totally 6 evaluations for 6 different linear combinations. Similarly, to the best ofour knowledge, expRK5s8 is also the only existing fifth-order stiffly accurate expRK meth-ods. From the construction of this scheme [15], one needs 8 stages. Among them, the firstfive internal stages are of the form (2.4). We note, however, that the last two internal stages U n and U n involves 2 and 3 different linear combinations (grouped in different brackets [ ] ) of ϕ k - functions (with different scaling factors) acting on different sets of vectors. Andeach stage ( U ni or u n + ) depends on all the preceding stages (except for the first stage U n ).Thus, this scheme must be also implemented in a sequential way (it also does not have anygroup of internal stages that can be computed simultaneously). Clearly, it requires totally 11evaluations (11 different linear combinations of ϕ functions). Remark 2.2
The resulting structures of the expRK schemes discussed in Remark 2.1 canbe explained by taking a closer look at their constructions presented in [8, 15]. 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-ordermethods (e.g., expRK4s5 ) 4 out of 9 order conditions have to be relaxed and for fifth-ordermethods (e.g., 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 schemesmust be implemented by computing each stage sequentially) and the very last stages usuallyinvolve different linear combinations of ϕ k -functions (with several different nodes c i as scal-ing factors) acting on not the same set of vectors, which then require additional sequentialevaluations. For more details, see Section 4 below. Inspired by the motivation and remarks in Section 2, we next present a stronger convergenceresult which later allows a construction of new efficient methods of high order. For this, wefirst recall the stiff order conditions for expRK methods up to order 5 (see [15, 17]). fficient exponential Runge–Kutta methods of high order 7 e n + = ˆ u n + − u ( t n + ) denote the local error of (1.2), i.e., the difference between thenumerical solution ˆ u n + obtained by (1.2) after one step starting from the ‘initial condition’ u ( t n ) and the corresponding exact solution u ( t n + ) of (1.1) at t n + .To simplify the notation in this section we set f ( t ) = g ( t , u ( t )) as done in [8], and addi-tionally denote G k , n = D k g ( t n , u ( t n )) be the k -th partial Fr´echet derivative (with respect to u )evaluated at u ( t n ) . Our results in [17] (Sect. 4.2) or [14] (Sect. 5.1) showed that˜ e n + = h ψ ( hA ) f ′ ( t n ) + h ψ ( hA ) f ′′ ( t n ) + h ψ ( hA ) f ′′′ ( t n ) + h ψ ( hA ) f ( ) ( t n )+ R n + O ( h ) (3.1)with the remaining terms R n = h s ∑ i = b i G , n ψ , i f ′ ( t n ) + h s ∑ i = b i G , n ψ , i f ′′ ( t n ) + h s ∑ i = b i G , n i − ∑ j = a i j G , n ψ , j f ′ ( t n )+ h s ∑ i = b i c i G , n (cid:0) u ′ ( t n ) , ψ , i f ′ ( t n ) (cid:1) + h s ∑ i = b i G , n ψ , i f ′′′ ( t n )+ h s ∑ i = b i G , n i − ∑ j = a i j G , n ψ , j f ′′ ( t n ) + h s ∑ i = b i G , n i − ∑ j = a i j G , n j − ∑ k = a jk G , n ψ , k f ′ ( t n ) (3.2) + h s ∑ i = b i G , n i − ∑ j = a i j c j G , n (cid:0) u ′ ( t n ) , ψ , j f ′ ( t n ) (cid:1) + h s ∑ i = b i c i G , n (cid:0) u ′ ( t n ) , ψ , i f ′′ ( t n ) (cid:1) + h s ∑ i = b i c i G , n (cid:0) u ′ ( t n ) , i − ∑ j = a i j G , n ψ , j f ′ ( t n ) (cid:1) + h s ∑ i = b i G , n (cid:0) ψ , i f ′ ( t n ) , ψ , i f ′ ( t n ) (cid:1) + h s ∑ i = b i c i G , n (cid:0) u ′′ ( t n ) , ψ , i f ′ ( t n ) (cid:1) + h s ∑ i = b i c i G , n (cid:0) u ′ ( t n ) , u ′ ( t n ) , ψ , i f ′ ( t n ) (cid:1) . Here, (and from now on) we use the abbreviations a i j = a i j ( hA ) , b i = b i ( hA ) , ϕ j , i = ϕ j ( c i hA ) and ψ j ( hA ) = s ∑ i = b i c j − i ( j − ) ! − ϕ j ( hA ) , j ≥ ψ j , i = ψ j , i ( hA ) = i − ∑ k = a ik c j − k ( j − ) ! − c ji ϕ j , i . (3.3b)Requiring a local error truncation ˜ e n + = O ( h ) results in the stiff order conditions formethods of order up to 5, which are displayed in Table 3.1 below.3.2 A stronger convergence resultThe convergence analysis of exponential Runge–Kutta methods is usually performed in theframework of analytic semigroups on a Banach space X with the following assumptions (seee.g. [8, 15]): Vu Thai Luan
Table 3.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 ψ l and ψ k , l are defined in (3.3). No. Stiff order condition Order1 ψ ( Z ) = ⇐⇒ ∑ si = b i ( Z ) c i = ϕ ( Z ) ψ ( Z ) = ⇐⇒ ∑ si = b i ( Z ) c i = ϕ ( Z ) ∑ si = b i ( Z ) J ψ , i ( Z ) = ψ ( Z ) = ⇐⇒ ∑ si = b i ( Z ) c i = ϕ ( Z ) ∑ si = b i ( Z ) J ψ , i ( Z ) = ∑ si = b i ( Z ) J ∑ i − j = a i j ( Z ) J ψ , j ( Z ) = ∑ si = b i ( Z ) c i K ψ , i ( Z ) = ψ ( Z ) = ⇐⇒ ∑ si = b i ( Z ) c i = ϕ ( Z ) ∑ si = b i ( Z ) J ψ , i ( Z ) = ∑ si = b i ( Z ) J ∑ i − j = a i j ( Z ) J ψ , j ( Z ) = ∑ si = b i ( Z ) J ∑ i − j = a i j ( Z ) J ∑ j − k = a jk ( Z ) J ψ , k ( Z ) = ∑ si = b i ( Z ) J ∑ i − j = a i j ( Z ) c j K ψ , j ( Z ) = ∑ si = b i ( Z ) c i K ψ , i ( Z ) = ∑ si = b i ( Z ) c i K ∑ i − j = a i j ( Z ) J ψ , j ( Z ) = ∑ si = b i ( Z ) B (cid:0) ψ , i ( Z ) , ψ , i ( Z ) (cid:1) = ∑ si = b i ( Z ) c i L ψ , i ( Z ) = Assumption 1. The linear operator A is the infinitesimal generator of an analytic semi-group e tA on X . This implies that k e tA k X ← X ≤ C , t ≥ ϕ k ( hA ) , the coefficients a i j ( hA ) and b i ( hA ) of the method are boundedoperators. Furthermore, the following stability bound (see [8, Lemma 1]) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) hA n ∑ j = e jhA (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X ← X ≤ C (3.5)holds uniformly for all n ≥ h > < nh ≤ T − t . Assumption 2 (for high-order methods). The solution u : [ t , T ] → X of (1.1) is sufficientlysmooth with derivatives in X and g : [ t , T ] → X is sufficiently often Fr´echet differentiablein a strip along the exact solution. All occurring derivatives are assumed to be uniformlybounded. fficient exponential Runge–Kutta methods of high order 9
Let e n + = u n + − u ( t n + ) denote the global error at time t n + . In [15], we have shownthat e n satisfies the recursion e n = h n − ∑ j = e ( n − j ) hA K j ( e j ) e j + n − ∑ j = e jhA ˜ e n − j , (3.6)where K j ( e j ) are bounded operators on X .Motivated by Remark 2.2, we now give a stronger convergence result (compared to thoseresults given in [8, 15]) in the sense that it requires relaxing only one order condition. Theorem 3.1 (Convergence) Let the initial value problem (1.1) satisfy Assumptions 1–2.Consider for its numerical solution an explicit exponential Runge–Kutta method (1.2) thatfulfills the order conditions of Table 3.1 up to order p ( ≤ p ≤ ) in a strong form withthe exception that only one condition ψ p ( Z ) = holds in a weakened form, i.e., ψ p ( ) = .Then, the method is convergent of order p. In particular, the numerical solution u n satisfiesthe error bound k u n − u ( t n ) k ≤ Ch p (3.7) uniformly on compact time intervals t ≤ t n = t + nh ≤ T with a constant C that dependson T − t , but is independent of n and h.Proof The proof can be carried out in a very similar way as done in [15, Theorem 4.2].In view of (3.1) and (3.2) and employing the assumptions of Theorem 3.1 on the orderconditions, we have R n = e n + = h p (cid:0) ψ p ( hA ) − ψ p ( ) (cid:1) G p − , n + h p + S n , (3.8)where G p − , n is defined in Section 3.1 and S n involves the terms multiplying h p + and higherorder in (3.1) (clearly, k S n k ≤ C ). Inserting (3.8) (with index n − j − n ) into (3.6)and using the fact that there exists a bounded operator ˜ ψ p ( hA ) such that ψ p ( hA ) − ψ p ( ) = ˜ ψ p ( hA ) hA yields e n = h n − ∑ j = e ( n − j ) hA K j ( e j ) e j + h p n − ∑ j = hA e jhA ˜ ψ p ( hA ) G p − , n − j − + h p + S n − j − . (3.9)Using (3.4), (3.5) and an application of a discrete Gronwall lemma shows (3.7). ⊓⊔ 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 3.1 in the contextof Theorem 3.1. It turns out that for methods of high order this will require an increase inthe number of stages s . However, we will have more degree of freedoms for constructingour desired methods as seen in Section 4 below. In addition, by relaxing only one ordercondition, we expect methods resulted from Theorem 3.1 to have better stability (and thusmay be more accurate) when integrating stiff systems (see Section 6). In this section, we will derive methods which have the following features: (i) containingmultiple internal stages U ni that are independent of each other (henceforth called parallelstages ) and share the same format (thereby allowing them to be implemented in parallel);(ii) involving less number of evaluations of the form (2.4) when compared to the existingmethods of the same orders (thus behaving like methods that use fewer number of stages s ).We first start with methods of order p ≤
3. When solving order conditions for thesemethods (requiring at least s = s = expRK2s2 (order 2, 2-stage) and expRK3s3 (order 3, 3-stage)mentioned in Section 2. We omit the details. Therefore, we will focus on the derivation ofnew methods of higher orders, namely, orders 4 and 5.4.1 A family of fourth-order methods with parallel stagesDeriving methods of order 4 requires solving the set of 7 stiff order conditions 1–7 in Ta-ble 3.1. First, we discuss on the required number of stages s . It is shown in [8, Sect.5.3] that s = b i ( Z ) as b i ( ) ). In other words, with s = ψ ( ) = ∑ si = b i ( ) c i = ϕ ( ) = /
24. Therefore, we consider s = b c + b c + b c + b c + b c = ϕ , (4.1a) b c + b c + b c + b c + b c = ϕ , (4.1b) b ( ) c + b ( ) c + b ( ) c + b ( ) c + b ( ) c = ϕ ( ) = / , (4.1c)and conditions 3, 5, 7 and 6 are b J ψ , + b J ψ , + b J ψ , + b J ψ , + b J ψ , = , (4.2a) b J ψ , + b J ψ , + b J ψ , + b J ψ , + b J ψ , = , (4.2b) b c K ψ , + b c K ψ , + b c K ψ , + b c K ψ , + b c K ψ , = , (4.2c) b Ja J ψ , + b J ( a J ψ , + a J ψ , ) + b J ( a J ψ , + a J ψ , + a J ψ , ) (4.2d) + b J ( a J ψ , + a J ψ , + a J ψ , + a J ψ , ) = . We now solve these order conditions. We note from (3.3b) that ψ , i = i − ∑ j = a i j c j − c i ϕ , i , ψ , i = i − ∑ j = a i j c j − c i ϕ , i (4.3)and thus ψ , = − c ϕ , = , ψ , = − c ϕ , = c = ψ , or ψ , must be nonzero as well (if both are zero then a = c c ϕ , = c c ϕ , , which is impossible since c > { ϕ , ϕ } are linearly independent). This fficient exponential Runge–Kutta methods of high order 11 strongly suggests that b = b = J and K . Next, we further observe that if b = ψ , = ψ , = a = a = b = b = b =
0. Using this sufficient condition we caneasily solve (4.1) to get b = − c ϕ + ϕ c ( c − c ) , b = − c ϕ + ϕ c ( c − c ) for any choice of distinct nodes c , c >
0, satisfying the condition c = c − c − . (4.4)Since b , b =
0, we must enforce ψ , = ψ , = ψ , = ψ , = a c + a c + a c = c ϕ , , (4.5a) a c + a c + a c = c ϕ , , (4.5b)and a c + a c + a c + a c = c ϕ , , (4.6a) a c + a c + a c + a c = c ϕ , . (4.6b)To satisfy conditions (4.2d), we further enforce a = a = ψ , = c = c ) a = − c c ϕ , + c ϕ , c ( c − c ) = , a = − c c ϕ , + c ϕ , c ( c − c ) = , (4.7)and thus we also need ψ , = ψ , = ψ , = a = c c ϕ , , (4.8a) a c + a c = c ϕ , . (4.8b)After fulfilling all the required order conditions in (4.1)–(4.2), we see from (4.6) and (4.8b)that either a or a and one of the coefficients among a , a , a can be taken as freeparameters. We now use them to construct parallel stages. Guided by (4.7) and (4.8a), wechoose a = U n is independent of U n so that both these stages only depend on U n , and choose a = U n is independent of U n so that both these stages onlydepend on the two preceding stages U n , U n (since a = a = a = c c ϕ , , a = − c c ϕ , + c ϕ , c ( c − c ) , a = − c c ϕ , + c ϕ , c ( c − c ) . (4.9) Putting altogether and rearranging terms in U ni , u n + as linear combinations of ϕ functions,we obtain the following family of 4th-order 6-stage methods (with the pairs of parallel stages { U n , U n } and { U n , U n } ), which will be called expRK4s6 : U n = u n + ϕ ( c hA ) c hF ( t n , u n ) , (4.10a) U n , k = u n + ϕ ( c k hA ) c k hF ( t n , u n ) + ϕ ( c k hA ) c k c hD n , k = , U n , j = u n + ϕ ( c j hA ) c j hF ( t n , u n ) + ϕ ( c j hA ) c j c − c h (cid:0) − c c D n + c c D n (cid:1) + ϕ ( c j hA ) c j c − c h (cid:0) c D n − c D n (cid:1) , j = , u n + = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) c − c h (cid:0) − c c D n + c c D n (cid:1) + ϕ ( hA ) c − c h (cid:0) c D n − c D n (cid:1) . (4.10d)For the numerical experiments given in Section 6, we choose c = c = , c = , c = which gives c = due to (4.4). Remark 4.1 ( A comparison with expRK4s5 ). As seen, expRK4s6 is resulted from weakeningonly condition 4 of Table 3.1 instead of weakening four conditions 4–7 as in the derivationof expRK4s5 . While the 5-stage method expRK4s5 requires 6 sequential evaluations in eachstep (as mentioned in Section 2), the new fourth-order 6-stage method expRK4s6 requiresonly 4 sequential evaluations, making it to behave like a 4-stage method. This is due to thefact expRK4s6 has the pairs of parallel stages { U n , U n } and { U n , U n } and all U ni withinthese pairs have the same format, i.e., same (one) linear combination of ϕ k ( c i hA ) v k , allowingthem to be computed in parallel or simultaneously (see Section 5).4.2 A family of fifth-order methods with parallel stagesConstructing fifth-order exponential Runge-Kutta methods needs much more effort as onehas to solve 16 order conditions in Table 3.1. As mentioned in Section 2, the only existingmethod of order 5 in the literature is expRK5s8 (see [15]) which requires s = expRK4s5 , this method does not have any parallel stages and must be implementedin a sequential way. It also does not satisfy the assumption on the order conditions statedin Theorem 3.1. Indeed, it was constructed by fulfilling conditions 1–7 in the strong formand weakening conditions 8–16 (9 out of 16 order conditions) with b i ( ) in place of b i ( Z ) .This resulted in the last two internal stages U n and U n that involve several different linearcombinations of ϕ k ( c i hA ) v k (with different scalings c , c , c of hA ), for which additionalcomputational 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 deriveparallel stages schemes with U ni involving only one linear combination of ϕ k ( c i hA ) v k , wehave to increase s >
9. To make it easier for readers to follow, we consider s =
10 first andlater employ the similar procedure to show that it is not possible to fulfill condition 11 ofTable 3.1 in the strong form (and thus not satisfying Theorem 3.1) with s = a) The case s =
10: Similarly to the derivation presented in Subsection 4.1, using (4.3),it strongly suggests b = b = b = b = b = b = fficient exponential Runge–Kutta methods of high order 13
16, 13, and 15 in their strong form. Using this, these conditions now read as b J ψ , + b J ψ , + b J ψ , = , (4.11a) b J ψ , + b J ψ , + b J ψ , = , (4.11b) b J ψ , + b J ψ , + b J ψ , = , (4.11c) b c K ψ , + b c K ψ , + b c K ψ , = , (4.11d) b c L ψ , + b c L ψ , + b c L ψ , = , (4.11e) b c K ψ , + b c K ψ , + b c K ψ , = , (4.11f) b B ( ψ , , ψ , ) + b B ( ψ , , ψ , ) + b B ( ψ , , ψ , ) = , (4.11g)respectively. And conditions 1, 2, 4, and 8 (weakened form) become b c + b c + b c = ϕ , (4.12a) b c + b c + b c = ϕ , (4.12b) b c + b c + b c = ϕ , (4.12c) b ( ) c + b ( ) c + b ( ) c = ϕ ( ) = / . (4.12d)Solving (4.12) gives b = c c ϕ − ( c + c ) ϕ + ϕ c ( c − c )( c − c ) , (4.13a) b = c c ϕ − ( c + c ) ϕ + ϕ c ( c − c )( c − c ) , (4.13b) b = c c ϕ − ( c + c ) ϕ + ϕ c ( c − c )( c − c ) (4.13c)where c , c , and c are distinct and positive nodes satisfying the algebraic equation c + c + c − c c + c c + c c + c c c = . (4.14)Clearly, b , b , b = ψ , j = ψ , j = ψ , j = ( j = , , ) (4.15)to satisfy (4.11) in the strong sense with arbitrary square matrices J , K , L and B . Next, weconsider conditions 6 and 10 taken into account that b i = i = , · · · ,
7) and (4.15), whichcan be now simplified as ∑ j = ( b Ja j + b Ja j + b Ja j ) J ψ m , j = ( m = , ) , (4.16)respectively. In order to satisfy the strong form of (4.16) one needs a j = a j = a j = ( j = , , ) (4.17)(this is again due to (4.3)) and ψ , j = ψ , j = ( j = , , ) . (4.18) With (4.17), we note that U n , U n , U n are independent of the internal stages U n , U n , U n .Taking into all the requirements above, one can easily see that conditions 12 and 14 are nowautomatically fulfilled. Therefore, the only remaining condition to satisfy is condition 11.Before working with condition 11, we first solve (4.15) using (4.17). For this, we observethat several coefficients a i j can be considered as free parameters. To have U n , U n , U n areindependent of each other, we choose a = a , = a , = . (4.19)The resulting systems of linear equations from (4.15) is then solved with the unique solution a i j = c i c k c l ϕ , i − c i ( c k + c l ) ϕ , i + c i ϕ , i c j ( c j − c k )( c j − c l ) , i = , , j , k , l ∈ { , , } , j = k = l (4.20)(i.e., c , c , c > b i = i = , · · · , ∑ i = b i J ∑ j = a i j J (cid:0) a j J ψ , + a j J ψ , + a j J ψ , (cid:1) = . (4.21)Since b , b , b =
0, coefficients a i j in (4.20) ( i ∈ { , , } , j ∈ { , , } ) are also nonzero,and that ψ , =
0, we must enforce a j = ( j = , , ) , i.e., a = a = a = ψ , = ψ , = a = a = j =
5) in (4.21) or both as this does not agree with the requirement ψ , = ψ , = a i j . When solving(4.23) (see (4.8)), we choose a = U n is independent of U n . This gives a = c c ϕ , , a = c c ϕ , . (4.24)When solving (4.18) (using (4.22)), we choose a = a = a = U n , U n , U n areindependent of each other. This results in the following 6 coefficients: a i j = − c i c k ϕ , i + c i ϕ , i c j ( c j − c k ) , i = , , j , k ∈ { , } , j = k (4.25)(i.e., c , c > a i j and b i into U ni , u n + and rewriting these stagesas linear combinations of ϕ functions, we obtain the following family of 5th-order 10-stage fficient exponential Runge–Kutta methods of high order 15 methods (with the groups of parallel stages { U n , U n } , { U n , U n , U n } , and { U n , U n , U n } )which will be called expRK5s10 : U n = u n + ϕ ( c hA ) c hF ( t n , u n ) , U n ℓ = u n + ϕ ( c ℓ hA ) c ℓ hF ( t n , u n ) + ϕ ( c ℓ hA ) c ℓ c hD n , ℓ = , U nm = u n + ϕ ( c m hA ) c m hF ( t n , u n ) + ϕ ( c m hA ) c m h (cid:0) c c ( c − c ) D n + c c ( c − c ) D n (cid:1) + ϕ ( c m hA ) c m h (cid:0) c ( c − c ) D n − c ( c − c ) D n (cid:1) , m = , , U nq = u n + ϕ ( c q hA ) c q hF ( t n , u n ) + ϕ ( c q hA ) c q h (cid:0) α D n + α D n + α D n (cid:1) + ϕ ( c q hA ) c q h (cid:0) β D n − β D n − β D n (cid:1) + ϕ ( c q hA ) c q h (cid:0) γ D n + γ D n + γ D n (cid:1) , q = , , u n + = u n + ϕ ( hA ) hF ( t n , u n ) + ϕ ( hA ) h (cid:0) α D n + α D n + α D n (cid:1) − ϕ ( hA ) h (cid:0) β D n + β D n + β D n (cid:1) + ϕ ( hA ) h (cid:0) γ D n + γ D n + γ D n (cid:1) , where α i = c k c l c i ( c i − c k )( c i − c l ) , β i = ( c k + c l ) c i ( c i − c k )( c i − c l ) , γ i = c i ( c i − c k )( c i − c l ) (4.26)with i ∈ { , , } for k , l ∈ { , , } , and i ∈ { , , } for k , l ∈ { , , } (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 = c = c = , c = c = , c = , c = , c = , and c = Remark 4.2 ( A comparison with expRK5s8 ). Although the new fifth-order method expRK5s10 has 10 stages (compared to 8 stages of expRK5s8 displayed in Section 2), its special struc-ture offers much more efficient for implementation. In particular, all U ni in this schemeinvolve only one linear combination of ϕ k ( c i hA ) v k which can be computed by one evalua-tion for each; and more importantly, due to the same format of multiple stages within each ofthe three groups { U n , U n } , { U n , U n , U n } , and { U n , U n , U n } (same linear combinationwith different inputs c i ), they can be computed simultaneously or implemented in paral-lel (see Section 5). This makes expRK5s10 to behave like a 5-stage method only, therebyrequiring only 5 sequential evaluations in each step. Moreover, while expRK5s8 requiresweakening 9 out of 16 order conditions of Table 3.1, expRK5s10 requires only one condition(number 8) held in the weakened form. Note that by following the similar way of deriving expRK5s10 , we can derive a scheme that satisfies all the stiff order conditions in Table 3.1 inthe 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 = (which does not work) : Clearly, in this case we have less degree offreedoms than the case s =
10 when solving the order conditions in Table 3.1. Nonethe-less, one can still proceed in a similar way as done for s =
10. Again, it strongly suggests b = b = b = b = b = b , b , b = ψ , j = ψ , j = ψ , j = ( j = , , ) (4.27)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 ∑ j = ( b Ja j + b Ja j + b Ja j ) J ψ m , j = ( m = , ) . (4.28) Again, due to the fact that ψ , , ψ , = ψ , or ψ , must be nonzero, one needsto enforce a j = a j = a j = ( j = , ) in (4.28). Using this to solve (4.27) for j = ψ , = ψ , = ψ , =
0) gives a unique solution (with c , c , c > a , a , a =
0, which then determines U n . Next, one can solve (4.27) for j = , U n , U n that are independent of U n , as well as are independent of each other, by requiringthe three free parameters a = a = a =
0. As a result, one gets a j , a j , a j = ( j = , ) .This immediately suggests ψ , j = ψ , j = ( j = , , ) to completely fulfill (4.28) witharbitrary square matrix J . With all of these in place, conditions 12 and 14 are automaticallyfulfilled, and condition 11 is now reduced to ∑ i = b i J ∑ j = a i j J (cid:0) a j J ψ , + a j J ψ , (cid:1) = . (4.29)Clearly, since b , b , b = a j , a j , a j = ( j = , , ) , and ψ , =
0, (4.29) can be sat-isfied in the strong sense only if we have one of the following conditions: a j = a j = a j = ψ , = , ( j = , , ) . Unfortunately, either of these requirements is in contradictionwith ψ , j = ψ , j = ( j = , , ) which is needed for conditions 6 and 10 mentioned above.For example, solving ψ , = ψ , = a , a = In this section, we present details implementation of the old and new fourth- and fifth-orderexpRK schemes ( expRK4s5 , expRK5s8 , expRK4s6 , expRK5s10 ) mentioned above.As mentioned in Section 2.1, we will use the MATLAB routine phipm simul iom (de-scribed in details in [19]) to implement expRK methods. In particular, given the followinginputs: an array of scaling factors t = [ ρ , · · · , ρ r ] with 0 < ρ < ρ < · · · < ρ r ≤ t couldbe a positive scalar), an n -by- n matrix M , and a set of column vectors V = [ V , . . ., v q ] (each v i is an n -by-1 vector), a tolerance tol , an initial value m for the dimension of the Krylovsubspace, and an incomplete orthogonalization length of iom , a call to this function phipm simul iom ( t , M , V , tol , m , iom ) (5.1)simultaneously computes the following r linear combinations L ρ i , V = ϕ ( ρ i M ) v + ϕ ( ρ i M ) ρ i v + ϕ ( ρ i M ) ρ i v + · · · + ϕ q ( ρ i M ) ρ qi v q , ≤ i ≤ r . (5.2)Note that, by setting V j = ρ ji v j ( j = , · · · , q ), (5.2) becomes (2.2). In other words, all thelinear combinations in (2.2) ( if V j are given instead of v j ) can be then computed at the sametime with one call (5.1) by using scaled vectors v j = V j / ρ ji for the input V .In the following, we set M = hA , ρ i = c i , v = hF ( t n , u n ) , d i = hD ni . (5.3)to simplify notations in presenting details of implementation of the fourth- and fifth-ordermethods mentioned above. When calling (5.1), we use tol = − , m = imo = Implementation of expRK4s5 ( c = c = c = , c = expRK4s5 requires a sequential implementation of the following 6 different linear combi-nations of the form (5.2), corresponding to 6 calls to phipm simul iom : fficient exponential Runge–Kutta methods of high order 17 (i) Evaluate L c , V with t = c , V = [ , v ] to get U n = u n + L c , V .(ii) Evaluate L c , V with t = c , V = [ , v , d / c ] to get U n = u n + L c , V .(iii) Evaluate L c , V with t = c , V = [ , v , d + d ] to get U n = u n + L c , V .(iv) Evaluate L c , V with t = c , V = [ , v , d + d − d , ( − d − d + d ) / c ] and(v) Evaluate L c , V with t = c , V = [ , , ( d + d − d ) / , ( − d − d + d )] to get U n = u n + L c , V + L c , V .(vi) Evaluate L , V with t = , V = [ , v , − d + d , d − d ] to get u n + = u n + L , V .Since d i = hD ni which depends on U ni , these are the 6 (sequential) evaluations. Implementation of expRK4s6 ( c = c = , c = c = , c = ): As discussed in Remark 4.1, expRK4s6 can be implemented like a 4-stage method by evaluating the following 4 sequen-tial evaluations, corresponding to 4 calls to phipm simul iom :(i) Evaluate L c , V with t = c , V = [ , v ] to get U n = u n + L c , V .(ii) Evaluate L c , V and L c , V simultaneously using t = [ c , c ] , V = [ , v , d / c ] to get both U n = u n + L c , V and U n = u n + L c , V .(iii) Evaluate L c , V and L c , V simultaneously with t = [ c , c ] , V = [ , v , − c ( c − c ) c d + c ( c − c ) c d , ( c − c ) c d − ( c − c ) c d ] to get both U n = u n + L c , V and U n = u n + L c , V .(iv) Evaluate L , V with t = , V = [ , v , c − c ( − c c d + c c d ) , c − c ( c d − c d )] to get u n + = u n + L , V . Implementation of expRK5s8 ( c = c = c = , c = , c = , c = , c = expRK5s8 requires a sequential implementation of 11 different linear com-binations of the form (5.2), corresponding to the following 11 calls to phipm simul iom :(i) Evaluate L c , V with t = c , V = [ , v ] to get U n = u n + L c , V .(ii) Evaluate L c , V with t = c , V = [ , v , d / c ] to get U n = u n + L c , V .(iii) Evaluate L c , V with t = c , V = [ , v , d / c ] to get U n = u n + L c , V .(iv) Evaluate L c , V with t = c , V = [ , v , ( − d + d ) / c , ( d − d ) / c ] to get U n = u n + L c , V .(v) Evaluate L c , V with t = c , V = [ , v , ( d − d ) / c , ( − d + d ) / c ] to get U n = u n + L c , V .(vi) Evaluate L c , V with t = c , V = [ , v , ( − d + d ) / c , ( d − d n ) / c ] and(vii) Evaluate L c , V with t = c , V = [ , , ( − d + d + d ) / c , ( d − d − d ) / c ] to get U n = u n + L c , V + L c , V .(viii) Evaluate L c , V with t = c , V = [ , v , ( − d + d + d ) / c , ( d − d − d ) / c , ( − d + d + d ) / c ] and(ix) Evaluate L c , V with t = c , V = [ , , ( − d + d + d ) / c , ( d − d − d ) / c , ( − d + d + d ) / c ] and(x) Evaluate L c , V with t = c , V = [ , , ( − d + d + d ) / c , ( d − d − d ) / c , ( − d + d + d ) / c ] to get U n = u n + L c , V + L c , V + L c , V .(xi) Evaluate L , V with t = , V = [ , v , d − d + d , − d + d − d , d − d + d ] to get u n + = u n + L , V . Implementation of expRK5s10 ( c = c = c = , c = c = , c = , c = , c = ,and c = 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 , V with t = c , V = [ , v ] to get U n = u n + L c , V .(ii) Evaluate L c , V and L c , V simultaneously using t = [ c , c ] , V = [ , v , d / c ] to get both U n = u n + L c , V and U n = u n + L c , V .(iii) Evaluate L c , V , L c , V , and L c , V simultaneously using t = [ c , c , c ] , V = [ , v , c c ( c − c ) d + c c ( c − c ) d , c ( c − c ) d − c ( c − c ) d ] to get U n = u n + L c , V , U n = u n + L c , V , U n = u n + L c , V .(iv) Evaluate L c , V , L c , V , and L c , V simultaneously using t = [ c , c , c ] , V = [ , v , α d + α d + α d , β d − β d − β d , γ d + γ d + γ d ] to get U n = u n + L c , V , U n = u n + L c , V , U n = u n + L c , V .(v) Evaluate L , V with t = , V = [ , v , α d + α d + α d , β d + β d + β d , γ d + γ d + γ d ] to get u n + = u n + L , V (coefficients α i , β i , γ i are given in (4.26)). In this section, we demonstrate the efficiency of our newly derived fourth- and fifth-orderexpRK time integration methods ( expRK4s6 , expRK5s10 ). Specifically, we will comparetheir performance against the existing methods of the same orders ( expRK4s5 , expRK5s8 )on several examples of stiff PDEs. All the numerical simulations are performed in MATLABon a single workstation, using an iMac 3.6 GHz Intel Core i7, 32 GB 2400 MHz DDR4. Example 6.1 ( A one-dimensional semilinear parabolic problem [8]): We first verify the or-der of convergence for the new derived fourth- and fifth-order expRK schemes ( expRK4s6 , expRK5s10 ) by considering the following PDE for u ( x , t ) , x ∈ [ , ] , t ∈ [ , ] , and subject tohomogeneous Dirichlet boundary conditions, ∂ u ( x , t ) ∂ t − ∂ u ( x , t ) ∂ x = + u ( x , t ) + Φ ( x , t ) , (6.1)whose exact solution is known to be u ( x , t ) = x ( − x ) e t for a suitable choice of the sourcefunction Φ ( x , t ) . Spatial discretization : For this example, we use standard second order finite differences with200 grid points. This leads to a very stiff system of the form (1.1) (with k A k ∞ ≈ . × ).The resulting system is then integrated on the time interval [ , ] using constant stepsizes, corresponding to the number of time steps N = , , , ,
64. The time integrationerrors at the final time t = expRK4s6 and expRK5s10 fully achieve orders 4 and 5, respectively. When com-pared to the old integrators of the same orders expRK4s5 and expRK5s8 , we note that, giventhe same number of time steps, expRK4s6 is slightly more accurate but is much faster than expRK4s5 (see the right diagram). In a similar manner, expRK5s10 gives almost identicalglobal errors but is also much faster than expRK5s8 . Finally, we observe that, for this exam-ple, for a global error that is larger than 10 − , the new fourth-order method expRK4s6 is thefastest one, and for more stringent errors, expRK5s10 is the fastest integrator. fficient exponential Runge–Kutta methods of high order 19 number of time steps -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 e rr o r expRK4s5 (Hochbruck-Ostermann)expRK4s6expRK5s8expRK5s10slope 4slope 5 CPU time -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 e rr o r expRK4s5 (Hochbruck-Ostermann)expRK4s6expRK5s8expRK5s10 Fig. 6.1
Order plots (left) and total CPU times (right) of expRK4s5 , expRK4s6 , expRK5s8 , and expRK5s10 when applied to (6.1). The global errors at time t = Example 6.2 ( A nonlinear Schrdinger equation [2, 5] ): We consider the following one-dimensional nonlinear Schrdinger (NLS) equation with periodic boundary conditions i ∂Ψ ( x , t ) ∂ t = − ∂ Ψ ( x , t ) ∂ x + (cid:0) V ( x ) + λ | Ψ ( x , t ) | (cid:1) Ψ ( x , t ) , Ψ ( − π , t ) = Ψ ( π , t ) , t ≥ Ψ ( , t ) = Ψ ( x ) , x ∈ [ − π , π ] (6.2)where the potential function V ( x ) = + sin ( x ) , the initial condition Ψ ( x ) = e sin ( x ) , andthe constant λ = Spatial discretization : For this example, we use a discrete Fourier transform F with ND =
128 modes, leading to a mildly stiff system of the form (1.1) with A = diag ( − i k ) , k = − ND + , · · · , ND = − , · · · , g ( t , u ) = − i F (( V ( x ) + λ | F − ( u ) | ) F − ( u ) . (6.3)Next, we integrate this system on the time interval [ , ] with constant step sizes, correspond-ing to the number of time steps N = , , , Ψ ( x , t ) of (6.2) is unknown, a reliable reference solution is computed by the stiff solver ode15s with ATOL = RTOL = − . Again, the time integration errors are measured in a discretemaximum norm at the final time t = expRK4s5 , expRK4s6 , expRK5s8 , and expRK5s10 ) as func-tions of the number of time steps (left) and the total CPU time (right). The left digram clearlyindicates that the two new integrators expRK4s6 and expRK5s10 achieve their correspond-ing expected orders 4 and 5. While expRK5s10 is a little more accurate than expRK5s8 , expRK4s6 is much more accurate than expRK4s5 for a given same number of time steps,meaning that it can take much larger time steps while achieving the same accuracy. More-over, the right precision digram displays the efficiency plot indicating that both expRK4s6 and expRK5s10 are much faster than their counterparts expRK4s5 and expRK5s8 , respec-tively. More specifically, a similar story is observed: for lower accuracy requirements, say error ∼ − , the new fourth-order method expRK4s6 is the most efficient, whereas for error ∼ − or tighter the new fifth-order method expRK5s10 is the most efficient.
64 128 256 512 1024 number of time steps -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 e rr o r expRK4s5 (Hochbruck-Ostermann)expRK4s6expRK5s8expRK5s10slope 4slope 5 CPU time -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 e rr o r expRK4s5 (Hochbruck-Ostermann)expRK4s6expRK5s8expRK5s10 Fig. 6.2
Order plots (left) and total CPU times (right) of expRK4s5 , expRK4s6 , expRK5s8 , and expRK5s10 when applied to Example 6.2. The errors at time t = Example 6.3 ( A 2D Gray–Scott model [3, 7] ): Consider the following two-dimensionalreaction-diffusion equation–the Gray–Scott equation model, for u = u ( x , y , t ) , v = v ( x , y , t ) on the square Ω = [ , L ] , (here, we choose L = .
5) subject to periodic boundary conditions ∂ u ∂ t = d u ∆ u − uv + α ( − u ) , ∂ v ∂ t = d v ∆ v + uv − ( α + β ) v , (6.4)where ∆ is the Laplacian operator, the diffusion coefficients d u = . , d v = .
01, and thebifurcation parameters α = . , β = . u ( x , y , ) = − e − (cid:0) ( x − L ) +( y − L ) (cid:1) , v ( x , y , ) = e − (cid:0) ( x − L ) + ( y − L ) (cid:1) . Spatial discretization : For this example, we use standard second order finite differencesusing 150 grid points in each direction with mesh width ∆ x = ∆ y = L / [ , ] using constant step sizes. In theabsence of an analytical solution of (6.4), a high-accuracy reference solution is computedusing the expRK4s6 method with a sufficient small time step. Errors are measured in a dis-crete maximum norm at the final time t = N = , , , , expRK4s6 is muchmore accurate than expRK4s5 and expRK5s10 is slightly more accurate than expRK5s8 . fficient exponential Runge–Kutta methods of high order 21
16 32 64 128 256 512 1024 number of time steps -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 e rr o r expRK4s5 (Hochbruck-Ostermann)expRK4s6expRK5s8expRK5s10slope 4slope 5 Fig. 6.3
Order plots of expRK4s5 , expRK4s6 , expRK5s8 , and expRK5s10 when applied to Example 6.3.The errors at time t = In Figure 6.4, we display the efficiency plot for which the time step sizes were chosen foreach integrator to obtain about same error thresholds 10 − i , i = , · · · ,
11 (The correspondingnumber of time steps for each integrator are displayed in Table 6.1. As seen, given aboutthe same level of accuracy, the new methods use smaller steps than the old ones of thesame order, meaning that they can take larger step sizes). Again, expRK4s6 is much fasterthan expRK4s5 and it is interesting that this new fourth-order method turns out to be themost efficient (although for error thresholds tighter than 10 − the new fifth-order method expRK5s10 seems to become the most efficient). CPU time -12 -11 -10 -9 -8 -7 -6 -5 -4 e rr o r expRK4s5 (Hochbruck-Ostermann)expRK4s6expRK5s8expRK5s10 Fig. 6.4
Total CPU times of expRK4s5 , expRK4s6 , expRK5s8 , and expRK5s10 when applied to Exam-ple 6.3. The time step sizes were chosen in such a way that each integrator achieves about the same errorthresholds 10 − i , i = , ··· ,
11. The errors at time t = − − − − − − − expRK4s5
18 36 66 121 215 385 685 expRK4s6
10 19 28 46 122 230 420 expRK5s8 expRK5s10
Table 6.1
The number of time steps taken to achieve about the same error thresholds 10 − i , i = , ··· , The numerical results presented on the three examples above clearly confirm the advan-tage of constructing parallel stages expRK methods based on Theorem 3.1, leading to moreefficient and accurate methods expRK4s6 and expRK5s10 . Acknowledgements
The author would like to thank Reviewer 1 for the valuable comments and helpfulsuggestions. He would like also to thank the National Science Foundation, which supported this researchunder award NSF DMS–2012022.
References
1. Al-Mohy, A.H., Higham, N.J.: Computing the action of the matrix exponential, with anapplication to exponential integrators. SIAM J. Sci. Comput. , 488–511 (2011)2. Berland, H., Skaflestad, B.: Solving the nonlinear Schr¨odinger equation using exponen-tial integrators. Tech. rep. (2005)3. Berland, H., Skaflestad, B., Wright, W.M.: Expint—a matlab package for exponentialintegrators. ACM Transactions on Mathematical Software (TOMS) (1), 4–es (2007)4. Caliari, M., Kandolf, P., Ostermann, A., Rainer, S.: The Leja method revisited: Back-ward error analysis for the matrix exponential. SIAM J. Sci. Comp. (3), A1639–A1661 (2016)5. Cazenave, T.: An introduction to nonlinear Schr¨odinger equations, vol. 22. Universi-dade Federal do Rio de Janeiro, Centro de Ciˆencias Matem´aticas e da (1989)6. Gaudreault, S., Pudykiewicz, J.: An efficient exponential time integration method forthe numerical solution of the shallow water equations on the sphere. J. Comput. Phys. , 827–848 (2016)7. Gray, P., Scott, S.: Autocatalytic reactions in the isothermal, continuous stirred tankreactor: Oscillations and instabilities in the system A + B → B ; B → C . ChemicalEngineering Science (6), 1087–1097 (1984)8. Hochbruck, M., Ostermann, A.: Explicit exponential Runge–Kutta methods for semi-linear parabolic problems. SIAM J. Numer. Anal. , 1069–1090 (2005)9. Hochbruck, M., Ostermann, A.: Exponential integrators. Acta Numerica , 209–286(2010)10. Ju, L., Wang, Z.: Exponential time differencing Gauge method for incompressible vis-cous flows. Communications in Computational Physics (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-dependentPDEs. Applied Numerical Mathematics , 91–103 (2017)13. Luan, V.T., Michels, D.: Explicit exponential Rosenbrock methods and their applicationin visual computing. (Revised) (2020) fficient exponential Runge–Kutta methods of high order 23
14. Luan, V.T., Ostermann, A.: Exponential B-series: The stiff case. SIAM J. Numer. Anal. , 3431–3445 (2013)15. Luan, V.T., Ostermann, A.: Explicit exponential Runge–Kutta methods of high orderfor parabolic problems. J. Comput. Appl. Math. , 168–179 (2014)16. Luan, V.T., Ostermann, A.: Exponential Rosenbrock methods of order five–construction, analysis and numerical comparisons. J. Comput. Appl. Math. , 417–431 (2014)17. Luan, V.T., Ostermann, A.: Stiff order conditions for exponential Runge–Kutta methodsof order five. In: H.B. et al. (ed.) Modeling, Simulation and Optimization of ComplexProcesses - HPSC 2012, pp. 133–143. Springer (2014)18. Luan, V.T., Ostermann, A.: Parallel exponential Rosenbrock methods. Comput. Math.Appl. , 1137–1150 (2016)19. Luan, V.T., Pudykiewicz, J.A., Reynolds, D.R.: Further development of efficient andaccurate time integration schemes for meteorological models. J. Comput. Phys. ,817–837 (2019)20. Michels, D.L., Luan, V.T., Tokman, M.: A stiffly accurate integrator for elastodynamicproblems. ACM Transactions on Graphics (TOG) (4), 116 (2017)21. Niesen, J., Wright, W.M.: Algorithm 919: A Krylov subspace algorithm for evaluatingthe ϕ -functions appearing in exponential integrators. ACM Trans. Math. Soft. (TOMS) (3), 22 (2012)22. Pieper, K., Sockwell, K.C., Gunzburger, M.: Exponential time differencing for mimeticmultilayer ocean models. J. Comput. Phys.398