Explicit exponential Runge–Kutta methods of high order for parabolic problems 1footnote 11footnote 1This work was supported by the FWF doctoral program ‘Computational Interdisciplinary Modelling’ W1227. The work of the first author was in addition supported by the Tiroler Wissenschaftsfond grant UNI-0404/1284.

# Explicit exponential Runge–Kutta methods of high order for parabolic problems 111This work was supported by the FWF doctoral program ‘Computational Interdisciplinary Modelling’ W1227. The work of the first author was in addition supported by the Tiroler Wissenschaftsfond grant UNI-0404/1284.

Vu Thai Luan and Alexander Ostermann Institut für Mathematik, Universität Innsbruck, Technikerstr. 19a,
A–6020 Innsbruck, Austria
###### Abstract

Exponential Runge–Kutta methods constitute efficient integrators for semilinear stiff problems. So far, however, explicit exponential Runge–Kutta methods are available in the literature up to order 4 only. The aim of this paper is to construct a fifth-order method. For this purpose, we make use of a novel approach to derive the stiff order conditions for high-order exponential methods. This allows us to obtain the conditions for a method of order 5 in an elegant way. After stating the conditions, we first show that there does not exist an explicit exponential Runge–Kutta method of order 5 with less than or equal to 6 stages. Then, we construct a fifth-order method with 8 stages and prove its convergence for semilinear parabolic problems. Finally, a numerical example is given that illustrates our convergence bound.

###### keywords:
exponential integrators, exponential Runge–Kutta methods, stiff order conditions, error bounds, semilinear parabolic problems, stiff problems
journal: J. of Computational and Applied Mathematics

## 1 Introduction

In this paper, we are concerned with the construction of high-order exponential Runge–Kutta methods for the time discretization of stiff semilinear problems

 u′(t)=F(u(t))=Au(t)+g(u(t)),u(t0)=u0 (1.1)

on the interval . Parabolic partial differential equations, written as abstract ordinary differential equations in some Banach space and their spatial discretizations are typical examples of such problems. Our main interest is the case of stiff problems where has a large norm or is an unbounded operator. On the other hand, the nonlinearity is assumed to satisfy a local Lipschitz condition with a moderate Lipschitz constant in a strip along the exact solution.

In recent years, exponential integrators turned out to be very competitive for stiff problems, see HLS98 ; KT05 . For a detailed overview of such integrators and their implementation, we refer to HO10 . The main idea behind these methods is to treat the linear part of problem (1.1) exactly and the nonlinearity in an explicit way. An early paper by Friedli F78 derived the class of exponential Runge–Kutta methods for non-stiff problems. He used classical Taylor series expansion for this purpose. For stiff problems, methods of orders up to four were constructed in HO05b . Very recently, in LO12b , we proposed a new and simple approach to derive stiff order conditions for exponential Runge–Kutta methods up to order five. With these order conditions at hand, we are going to construct a fifth-order method in this paper. We first show that there does not exist an explicit exponential Runge–Kutta method of order 5 with less than or equal to 6 stages. On the other hand, we are able to construct a family of fifth-order methods with 8 stages. These methods satisfy some of the order conditions only in a weakened form, however, this will be sufficient for obtaining convergence results for semilinear parabolic problems. We note that exponential Runge–Kutta methods may also be applied to the solution of other problems, see for example DP11 ; D09 .

The paper is organized as follows. In section 2, we recall exponential Runge–Kutta methods for our further analysis. Our abstract framework is given in section 3. In section 4, we recall the stiff order conditions derived in LO12b . A new convergence result is also given in this section. Section 5 is devoted to the construction of exponential Runge–Kutta methods of order 5. In section 6, a numerical experiment is presented to illustrate the order of the new method. The main results of the paper are Theorem 4.1, Theorem 5.1 and the new integrator .

Throughout the paper, will denote a generic constant that may have different values at different occurrences.

## 2 Exponential Runge–Kutta methods in a restated form

For solving (1.1), we consider the following class of -stage explicit exponential Runge–Kutta methods HO05b , reformulated as in LO12b

 Uni =un+cihnφ1(cihnA)F(un)+hni−1∑j=2aij(hnA)Dnj, 1≤i≤s, (2.1a) un+1 =un+hnφ1(hnA)F(un)+hns∑i=2bi(hnA)Dni (2.1b) with Dni=g(Uni)−g(un),2≤i≤s. (2.1c)

As usual, denotes the time step size and the are the nodes. The internal stages approximate the exact solution at . For the sake of completeness, one can define and . However, we note that these quantities do not enter the scheme anyway. The coefficients and are chosen as linear combinations of the entire functions and scaled versions thereof. These functions are given by

 φ0(z)=ez,φk(z)=∫10e(1−θ)zθk−1(k−1)!dθ,k≥1. (2.2)

They satisfy the recurrence relation

 φk+1(z)=φk(z)−φk(0)z,k≥0. (2.3)

The reformulated scheme (2.1) can be implemented more efficiently than its original form, and it offers many advantages for the error analysis (see LO12b ).

As shown in HO05b , exponential Runge–Kutta methods are invariant under the transformation of non-autonomous problems

 u′(t)=F(t,u(t))=Au(t)+g(t,u(t)) (2.4)

to (1.1) by adding . Scheme (2.1) can be applied to (2.4) by replacing in (2.1a), (2.1b) by and in (2.1c) by

 Dni=g(tn+cihn,Uni)−g(tn,un).

As a consequence of this invariance we will consider henceforth the autonomous case only. However, all stated results stay mutatis mutandis valid for the non-autonomous problem (2.4) as well.

## 3 Analytical framework

Our analysis will be based on an abstract framework of analytic semigroups on a Banach space with norm . Background information on semigroups can be found in the monograph PAZY83 .

Throughout the paper we consider the following assumptions.

Assumption 1. The linear operator is the infinitesimal generator of an analytic semigroup on .

This assumption implies that there exist constants and such that

 ∥etA∥X←X≤Meωt,t≥0. (3.1)

In particular, the expressions and consequently the coefficients and of the method are bounded operators, see (2.2) for . This property is crucial in our proofs.

Under Assumption 1, the following stability bound was proved in (HO05a, , Lemma 1): There exists a constant such that

 ∥∥ ∥∥hAn∑j=1ejhA∥∥ ∥∥X←X≤C. (3.2)

This bound holds uniformly for all and with . We will employ this bound later on.

For high-order convergence results, we require the following regularity assumption.

Assumption 2. We suppose that (1.1) possesses a sufficiently smooth solution with derivatives in and that is sufficiently often Fréchet differentiable in a strip along the exact solution. All occurring derivatives are assumed to be uniformly bounded.

Assumption 2 implies that is locally Lipschitz in a strip along the exact solution. It is well known that semilinear reaction-diffusion-advection equations can be put into this abstract framework, see H81 .

## 4 Local error, stiff order conditions and convergence results for fifth-order methods

In this section, we recall some notations and results from LO12b that will be used later. We further give a result which allows us to relax the stiff order conditions so that high-order convergence is still guaranteed.

### 4.1 Local error

For the error analysis of scheme (2.1), we consider one step with initial value on the exact solution, i.e.

 ˆUni =~un+cihnφ1(cihnA)F(~un)+hni−1∑j=2aij(hnA)ˆDnj, (4.1a) ^un+1 =~un+hnφ1(hnA)F(~un)+hns∑i=2bi(hnA)ˆDni (4.1b)

with

 ˆDni=g(ˆUni)−g(~un),ˆUni≈u(tn+cihn). (4.2)

Let denote the th derivative of the exact solution of (1.1), evaluated at time . For we use the common notation for simplicity. We further denote the th derivative of with respect to by .

Let denote the local error, i.e., the difference between the numerical solution after one step starting from and the corresponding exact solution of (1.1) at , and let

 ψj(hnA)=s∑i=2bi(hnA)cj−1i(j−1)!−φj(hnA),j≥2. (4.3)

For simplicity, we also use the abbreviations , and

 ψj,i=ψj,i(hnA)=i−1∑k=2aikcj−1k(j−1)!−cjiφj,i. (4.4)

In LO12b , we have shown that

 ~en+1 =h2nψ2(hnA)Ln+h3nψ3(hnA)Mn+h4nψ4(hnA)Nn (4.5) +h5nψ5(hnA)Pn+Rn+O(h6n)

with

 Ln =g′(~un)~u′n, (4.6) Mn =g′(~un)~u′′n+g′′(~un)(~u′n,~u′n), Nn =g′(~un)~u(3)n+3g′′(~un)(~u′n,~u′′n)+g(3)(~un)(~u′n,~u′n,~u′n), Pn =g′(~un)~u(4)n+3g′′(~un)(~u′′n,~u′′n)+4g′′(~un)(~u′n,~u(3)n) +6g(3)(~un)(~u′n,~u′n,~u′′n)+g(4)(~un)(~u′n,~u′n,~u′n,~u′n),

and the remaining terms

 Rn=h3ns∑i=2big′(~un)ψ2,iLn+h4ns∑i=2big′(~un)ψ3,iMn +h4ns∑i=2big′(~un)i−1∑j=2aijg′(~un)ψ2,jLn+h4ns∑i=2bicig′′(~un)(~u′n,ψ2,iLn) +h5ns∑i=2big′(~un)ψ4,iNn+h5ns∑i=2big′(~un)i−1∑j=2aijg′(~un)ψ3,jMn +h5ns∑i=2big′(~un)i−1∑j=2aijg′(~un)j−1∑k=2ajkg′(~un)ψ2,kLn (4.7) +h5ns∑i=2big′(~un)i−1∑j=2aijcjg′′(~un)(~u′n,ψ2,jLn)+h5ns∑i=2bicig′′(~un)(~u′n,ψ3,iMn) +h5ns∑i=2bicig′′(~un)(~u′n,i−1∑j=2aijg′(~un)ψ2,jLn)+h5ns∑i=2bi2!g′′(~un)(ψ2,iLn,ψ2,iLn) +h5ns∑i=2bic2i2!g′′(~un)(~u′′n,ψ2,iLn)+h5ns∑i=2bic2i2!g(3)(~un)(~u′n,~u′n,ψ2,iLn).

The remainder term in the asymptotic expansion (4.5) has the form

 s∑i=2bi(hA)Rni−Rn (4.8)

with

 Rni =h6n∫10(1−θ)44!g(5)(~un+θhnVi)(Vi,…,Vi5 times)dθ, Rn =h6n∫10e(1−θ)hnA∫10θ5(1−s)44!g(5)(~un+sθhnV)(V,…,V5 times)dsdθ.

Note that the arguments

 Vi =1hn(ˆUni−~un)=ciφ1(cihnA)F(~un)+i−1∑j=2aij(hnA)ˆDnj, V =1θhn(u(tn+θhn)−u(tn))

remain bounded as . Therefore, (4.8) is bounded by , where the constant only depends on values that are uniformly bounded by Assumptions 1 and 2. This justifies the notation for the remainder in (4.5).

### 4.2 Stiff order conditions for methods of order five

By zeroing the corresponding terms in (4.5), the stiff order conditions for methods of order five can easily be identified, see (LO12b, , Table 1). We note, however, that the first and the third condition of that table are automatically satisfied in our context since they were used to derive our reformulated scheme (2.1). As we will heavily need the order conditions for constructing a method of order five, we display them again in Table 1 below, where we have omitted the two redundant conditions.

### 4.3 Stability and convergence results

It was shown in (LO12b, , Theorem 4.1) that the numerical method (2.1) is stable and converges up to order five, if all conditions of Table 1 are fulfilled. However, a method satisfying all these conditions will have a plenty of stages. We therefore proceed differently. Motivated by the approach taken in HO05b , we relax the order conditions 8–16 in a way described below.

First, we consider the case of constant step size, i.e.  for all . Let denote the global error, i.e., the difference between the numerical and the exact solution of (1.1), and let and denote the differences between the two numerical solutions obtained by the schemes (2.1) and (4.1), respectively. It is easy to see that the global error at time satisfies

 en+1=^en+1+~en+1. (4.9)

Subtracting (4.1b) from (2.1b) gives

 ^en+1=ehAen+hTn (4.10)

with

 (4.11)

Denoting ,   and using the Taylor series expansion of , we get

 g(un)−g(~un) =Gnen+Rn, (4.12a) g(Uni)−g(ˆUni) =GniˆEni+Rni (4.12b)

with remainders and

 Rn =∫10(1−θ)g′′(~un+θen)(en,en)dθ, (4.13a) Rni =∫10(1−θ)g′′(ˆUni+θˆEni)(ˆEni,ˆEni)dθ. (4.13b)

Employing Assumption 2 shows that

 ∥Rn∥≤C∥en∥2,∥Rni∥≤C∥ˆEni∥2, (4.14)

as long as and remain in a sufficiently small neighborhood of 0.

Subtracting (4.12a) from (4.12b), we get

 Dni−ˆDni=GniˆEni+Rni−(Gnen+Rn). (4.15)

Inserting (4.12a) and (4.15) into (4.11), we now obtain

 Tn=b1(hA)(Gnen+Rn)+s∑i=2bi(hA)(GniˆEni+Rni) (4.16)

with  .

###### Lemma 4.1.

Under Assumptions 1 and 2, there exist bounded operators on such that

 Tn=Kn(en)en. (4.17)
###### Proof.

Subtracting (4.1a) from (2.1a), using and employing (4.12b), we obtain

 ˆEni=Uni−ˆUni=ecihAen+hi−1∑j=1aij(hA)(GnjˆEnj+Rnj). (4.18)

Solving recursion (4.18), using (4.13) with an induction argument and inserting the obtained result into (4.16) yields the formula (4.17). ∎

In view of (4.9), (4.10) and (4.17), we get

 en+1=ehAen+hKn(en)en+~en+1. (4.19)

Solving recursion (4.19) and using finally yields

 en=hn−1∑j=0e(n−j)hAKj(ej)ej+n−1∑j=0ejhA~en−j. (4.20)

We are now ready to give a convergence result for constant step sizes that only requires a weakened form of conditions 8–16 of Table 1.

###### Theorem 4.1.

Let the initial value problem (1.1) satisfy Assumptions 1–2. Consider for its numerical solution an explicit exponential Runge–Kutta method (2.1) that fulfills the order conditions 1–7 of Table 1 (up to order four). Further assume that condition 8 holds in a weakened form with and the remaining conditions of order five hold in a weaker form with in place of . Then, the method is convergent of order 5. In particular, the numerical solution satisfies the error bound

 ∥un−u(tn)∥≤Ch5 (4.21)

uniformly on compact time intervals   with a constant that depends on , but is independent of and .

###### Proof.

In view of (4.7) and (4.5), under the assumptions of Theorem 4.1, we obtain

 ~en+1=h5(ψ5(hA)−ψ5(0))Pn+h5s∑i=2(bi(hA)−bi(0))Qn+h6Sn, (4.22)

where denotes the terms multiplying in (4.7) and Inserting (4.22) (with index in place of ) into (4.20) yields

 en =hn−1∑j=0e(n−j)hAKj(ej)ej+h5n−1∑j=0ejhA(ψ5(hA)−ψ5(0))Pn−j−1 (4.23) +h5n−1∑j=0ejhAs∑i=2(bi(hA)−bi(0))Qn−j−1+h6Sn−j−1.

We can now employ the same techniques that were used in the proof of (HO05b, , Theorem 4.7). The bound (4.21) follows from the stability bound (3.2), the bound (HO05b, , Lemma 4.8) and an application of a discrete Gronwall lemma. ∎

Under further assumptions on the regularity of the solution or on the step size sequence, it is possible to generalize Theorem 4.1 to variable step sizes. We omit the details.

## 5 Construction of methods of order 5

In this section, we will construct exponential Runge–Kutta methods of order five based on the weakened form of the order conditions 8–16, as mentioned in Theorem 4.1. In fact, one needs to solve a system of 16 equations (corresponding to 16 order conditions of Table 1) for unknown coefficients The obtained coefficients of the method will be displayed in the following Butcher tableau:

 c2c3a32⋮⋮⋱csas2…as−1,sb2…bs−1bs

We are interested in finding methods with a minimal number of stages. In order to answer this question, we first state the following two technical lemmas.

###### Lemma 5.1.

Let be linearly independent (analytic) functions and assume that for all square matrices of the same format

 ℓ∑i=1Xi(Z)JYi(Z)=0 (5.1)

for given (analytic) functions . Then for all .

###### Proof.

Let Inserting these matrices into (5.1) shows that for all . Since the set is linearly independent, we get for all , which proves the lemma. ∎

###### Lemma 5.2.

Let be linearly independent analytic functions and assume that for all square matrices of the same format

 ℓ∑i=1Xi(Z)Jm∑j=1Yij(Z)JWj(Z)=0 (5.2)

for given (analytic) functions . Then, for all , we have

 m∑j=1Yij(μ)Wj(ν)=0. (5.3)

for all .

###### Proof.

We choose the matrices

 Z=⎡⎢⎣λ000μ000ν⎤⎥⎦,J=⎡⎢⎣010001100⎤⎥⎦

and insert them into (5.2). This shows that for all . Since the set is linearly independent, we obtain for all . ∎

With the results of Lemmas 5.1 and 5.2 at hand, we are now in the position to state the main result of this section.

###### Theorem 5.1.

There does not exist an explicit exponential Runge–Kutta method of order 5 with less than or equal to 6 stages.

###### Proof.

Note that scheme (2.1) reduces to a classical Runge–Kutta method for . A famous result by Butcher shows that an explicit Runge–Kutta method of order five with stages does not exist (see Butcher64 ). We deduce from this result that, for , no explicit exponential Runge–Kutta method of order 5 exists. Therefore, we now consider the case . Based on Theorem 4.1, one needs to check conditions 1–7 of Table 1. For , conditions 1, 2 and 4 are

 b2c2+b3c3+b4c4+b5c5+b6c6 =φ2, (5.4a) b2c22+b3c23+b4c24+b5c25+b6c26 =2φ3, (5.4b) b2c32+b3c33+b4c34+b5c35+b6c36 =6φ4, (5.4c)

and conditions 3, 5, 7 and 6 read

 b2Jψ2,2+b3Jψ2,3+b4Jψ2,4+b5Jψ2,5+b6Jψ2,6 =0, (5.5a) b2Jψ3,2+b3Jψ3,3+b4Jψ3,4+b5Jψ3,5+b6Jψ3,6 =0, (5.5b) b2c2Kψ2,2+b3c3Kψ2,3+b4c4Kψ2,4+b5c5Kψ2,5+b6c6Kψ2,6 =0, (5.5c)

and

 b3Ja32Jψ2,2 +b4J(a42Jψ2,2+a43Jψ2,3) (5.6) +b5J(a52Jψ2,2+a53Jψ2,3+a54Jψ2,4) +b6J(a62Jψ2,2+a63Jψ2,3+a64Jψ2,4+a65Jψ2,5)=0,

respectively. From (4.4), one derives

 ψ2,i=i−1∑j=2aijcj−c2iφ2,i,ψ3,i=i−1∑j=2aijc2j2!−c3iφ3,i. (5.7)

We note for later use that (otherwise and we are back to the case of 5 stages) and

 ψ2,2=−c22φ2,2≠0,ψ3,2=−c32φ3,2≠0. (5.8)