On the implementation of exponential methods for semilinear parabolic equations

# On the implementation of exponential methods for semilinear parabolic equations

María López-Fernández 222Departamento de Matemáticas, Universidad Autónoma de Madrid, Madrid, Spain.  E-mail: maria.lopez@uam.es. Supported by DGI-MCYT under projects MTM 2005-00714 and MTM 2007-63257, cofinanced by FEDER funds, and the SIMUMAT project S-0505/ESP/0158 of the Council of Education of the Regional Government of Madrid (Spain).
July 28, 2008.
July 27, 2008
###### Abstract

The time integration of semilinear parabolic problems by exponential methods of different kinds is considered. A new algorithm for the implementation of these methods is proposed. The algorithm evaluates the operators required by the exponential methods by means of a quadrature formula that converges like , with the number of quadrature nodes. The algorithm allows also the evaluation of the associated scalar mappings and in this case the quadrature converges like . The technique is based on the numerical inversion of sectorial Laplace transforms. Several numerical illustrations are provided to test the algorithm.

Key words. exponential methods, numerical inverse Laplace transform, semilinear parabolic equations.

AMS subject classifications. 65M70, 65R10.

## 1 Introduction

The good numerical results obtained from the application of exponential methods to the time integration of stiff semilinear problems, have motivated much interest on these kind of methods during the last years, see for instance [1, 2, 3, 9, 10, 11]. The problems under consideration can be written in the abstract format

 u′(t)=Au(t)+f(t,u(t)),u(0)=u0,0≤t≤T, \hb@xt@.01(1.1)

where is a linear operator representing the highest order differential terms and is a lower order nonlinear operator. The solution to the initial value problem (LABEL:ivp) is then given by the variation of constants formula and most of the exponential methods considered in the literature are constructed from this representation of the solution.

Let us consider for instance the family of multistep exponential methods developed in [2] and briefly reviewed in Section LABEL:subsec_msm. Given a stepsize , and approximations , , , the -step method approximates the solution of (LABEL:ivp) at by

 un+k=ϕ0(k,hA)un+hk−1∑j=0ϕj+1(k,hA)Δjfn, \hb@xt@.01(1.2)

where , denotes the standard forward difference operator, and for and ,

 ϕ0(k,λ)=ekλ,ϕj(k,λ)=∫k0e(k−σ)λ(σj−1)dσ,1≤j≤k. \hb@xt@.01(1.3)

As we can see in (LABEL:multistep_intro), these methods require the evaluation of , , for defined in (LABEL:phis_multistep). This is in fact the main difficulty in the implementation of the methods in (LABEL:multistep_intro) and, in general, of exponential methods, since they typically require the evaluation of vector-valued mappings , with the time step in the discretization and either

 ϕ(λ)=emλ,λ∈C, \hb@xt@.01(1.4)

or

 ϕ(λ)=∫m0e(m−σ)λp(σ)dσ, \hb@xt@.01(1.5)

with an integer and a polynomial. The values of and in (LABEL:Phigeneral) depend on the method. For instance, in the case of the methods in (LABEL:multistep_intro), it is clear from (LABEL:phis_multistep) that , the number of steps of the method, and

 p(σ)=(σj−1)=σ(σ−1)…(σ−j+2)(j−1)!,1≤j≤k.

For the explicit exponential Runge–Kutta methods constructed in [10] it is and

 p(σ)=σk−1(k−1)!,k≥1,

as we can see in Section LABEL:subsec_rkm, (LABEL:phis_rk).

In the present paper we propose a way to evaluate the operators in (LABEL:Phiexponencial) and (LABEL:Phigeneral) when in (LABEL:ivp) is the infinitesimal generator of an analytic semigroup in a Banach space . Thus, we will assume that is sectorial, i.e., is a densely defined and closed linear operator on and there exist constants , , and an angle , such that the resolvent fulfils

 ∥(zI−A)−1∥≤M|z−γ|,  for |arg(z−γ)|<π−δ. \hb@xt@.01(1.6)

As a side product we also obtain an accurate way to evaluate the mappings in (LABEL:Phigeneral) at the scalar level, which is itself a well-known problem in numerical analysis. This is exemplified in [11] with the mapping

 φ(λ)=eλ−1λ, \hb@xt@.01(1.7)

which is required for instance by the exponential Runge–Kutta methods of [10]. The evaluation of for small by using formula (LABEL:vfi_1) suffers from cancellation error. On the other hand, the use of a truncated Taylor expansion only works well for small enough . This implies that for some intermediate values of it is not very much clear how to choose the proper formula and moreover both of them could lose accuracy. For inside a sector , these difficulties can be overcome by writing in the format (LABEL:Phigeneral), with and . By doing so, we will be able to evaluate , independently of the size of , by using essentially the same technique developed in principle to evaluate the vector-valued mapping .

In the recent literature, several alternatives have been proposed to evaluate the required operators or alternatively their action on given vector, based on a Krylov approach [7, 8], on Padé approximants [1] or on a suitable contour integral representation of the mappings by means of the Cauchy integral formula [11]. In particular, in [11], the goal is both the evaluation of the mappings at the scalar level, assuming that a diagonalization of is available, but also at the operator level, since it also allows the evaluation of the operators working with the full matrix . However, despite the good computational results reported in [11], the way of selecting the parameters involved in the quadrature formulas is not very much clear and they depend very strongly on the equation considered and the spatial discretization parameters. The algorithm we propose is derived by using Laplace transformation formulas and is finally based on a suitable contour integral representation of the mappings , too. However, the quadrature formulas we use borrow their parameters from the method in [14] for the numerical inversion of sectorial Laplace transform, where a rigorous analysis of the error is performed together with an optimization process to choose the different parameters involved in the approximation.

In order to apply the quadrature formulas developed in [14], we derive a representation of the operators in (LABEL:Phiexponencial) and (LABEL:Phigeneral) as the inverses of suitable Laplace transforms at certain values of the original variable . These Laplace transforms have all the form

 Φ(z,hA)=R(z)(zI−hA)−1, \hb@xt@.01(1.8)

with a scalar rational mapping of . Due to (LABEL:sectorial), the mappings turn out to be sectorial in the variable , i.e., there exist constants and , possibly different from the constants in (LABEL:sectorial), such that

 Φ(z,hA) is analytic for z in the sector |arg(z−γ)|<π−δ and there∥Φ(z,hA)∥≤M|z−γ|ν,for some% ν≥1. \hb@xt@.01(1.9)

In this way, we reduce the problem of computing to the inversion of a sectorial Laplace transform of the form (LABEL:LTPhi_general). We then use the method for the numerical inversion of sectorial Laplace transforms developed in [14] and reviewed in Section LABEL:sec_LT. The inversion method consists on a quadrature formula to discretize the inversion formula (LABEL:invLT), so that we finally approximate

 ϕ(hA)≈K∑ℓ=−KwℓemzℓΦ(zℓ,hA)=K∑ℓ=−KwℓemzℓR(zℓ)(zℓI−hA)−1, \hb@xt@.01(1.10)

with the quadrature weights and nodes given in (LABEL:pesos_nodos). The convergence results in [14] assure an error estimate in the approximation (LABEL:quadgeneral) at least like . Further, by following this approach, our selection of parameters in the implementation of exponential methods depends only on in (LABEL:sectorial), being independent of and . In the particular case that we want to evaluate at a scalar in the sector , the approximation becomes simply

 ϕ(λ)≈K∑ℓ=−KwℓemzℓΦ(zℓ,λ)=K∑ℓ=−KwℓemzℓR(zℓ)zℓ−λ \hb@xt@.01(1.11)

and the convergence rate improves to an .

The approximation in (LABEL:quadgeneral) allows the computation of all the operators required by an exponential method before the time-stepping begins, so that only the products matrixvector need to be implemented at every time step. However, formula (LABEL:quadgeneral) requires the computation and storage of all the inverses , . Even if is a sparse matrix and these inversions can be performed efficiently, the storage of the resulting full matrices and the subsequent products matrixvector can become prohibitive for large problems. In this situation, (LABEL:quadgeneral) should be combined with a data sparse procedure to approximate the resolvent operators, as it is proposed in [4, 5].

Another way to avoid the computation and storage of all the resolvents in (LABEL:quadgeneral) could be the application of the formula to evaluate the action on a given vector , instead of the full operator . In this case, it is not necessary the computation of the full inverses , but the resolution of the linear systems

 (zℓI−hA)x=v.

However, solving all these linear systems for all the quadrature nodes in (LABEL:quadgeneral) at every time step can become quite expensive and, in this sense, the Krylov approach developed in [7] appears as a more competitive alternative. In this situation, only formula (LABEL:quadscalar) can be useful, in order to evaluate the mappings at the eigenvalues of the Hessenberg matrices involved in the Krylov approximation.

By using (LABEL:quadgeneral) we are in fact computing an approximation to the numerical solution of (LABEL:ivp) provided by an exponential method. Thus, the global error after applying (LABEL:quadgeneral) to the time integration of (LABEL:ivp) can be split into the error in the time integration by the “pure” exponential method and the deviation from the numerical solution introduced by the approximation (LABEL:quadgeneral) of the operators . The error in the time integration for the exponential integrators considered in the present paper is analyzed in [2] and [10] (see Section LABEL:sec_expmethods), and the quadrature error is analyzed in [14] (see Section LABEL:sec_LT, Theorem LABEL:thm:err). In order to visualize the effect of this approximation, we show in Section LABEL:sec_experiments the performance of our implementation for several problems with known exact solution and moderate size after the spatial discretization. In the error plots provided we can observe that the error coincides with the expected error for the exact exponential integrators up to high accuracy for quite moderate values of in (LABEL:quadgeneral), i.e., the error induced by the quadrature (LABEL:quadgeneral) is negligible compared to the error in the time integration.

Finally we notice that the matrix exponential and also certain rational approximations to it originating from Runge–Kutta schemes have already been successfully approximated by using this approach [12, 13, 14].

The paper is organized as follows. In Section LABEL:sec_expmethods we consider the class of explicit multistep exponential methods proposed in [2] and the explicit exponential Runge–Kutta methods in [10]. Section LABEL:sec_LT is a review of the method for the numerical inversion of sectorial Laplace transforms presented in [14]. In Section LABEL:sec_eval we deduce a representation for the operators required in the implementation of these integrators in terms of suitable Laplace transforms and apply the inversion method to the implementation of exponential methods. In Section LABEL:sec_evalscalar we consider with some detail the evaluation of the associated mappings at the scalar level and present some numerical results. We finally test our algorithm with several examples in Section LABEL:sec_experiments, where we implemented (LABEL:quadgeneral) by using the full matrices.

## 2 Exponential methods

In this section we review some of the exponential methods in the recent literature. In this way, we consider the class of explicit multistep exponential integrators developed in [2] and the exponential Runge–Kutta methods of [10], which are explicit as well.

### 2.1 Multistep exponential methods

In [2], a class of explicit exponential methods is constructed for the time integration of problems of the form (LABEL:ivp) with the infinitesimal generator of a -semigroup , , of linear and bounded operators in a Banach space .

As we already stated in the Introduction, we will restrict ourselves to the case of in (LABEL:ivp) sectorial. Then, for , the fractional powers are well defined for in (LABEL:sectorial), and endowed with the graph norm is a Banach space [6]. In this situation, the nonlinearity in (LABEL:ivp) is assumed to be defined on , for some , and to be locally Lipschitz in a strip along the exact solution. Thus, there exists such that

 ∥f(t,η)−f(t,ξ)∥≤L∥η−ξ∥α,η,ξ∈Xα,0≤t≤T, \hb@xt@.01(2.1)

for .

For , the -step method is constructed from the variation of constants formula in the interval as follows. Taking a stepsize , , and the corresponding time levels , , the solution of (LABEL:ivp) at is given by

 u(tn+k)=ekhAu(tn)+h∫k0e(k−σ)hAf(tn+σh,u(tn+σh))dσ. \hb@xt@.01(2.2)

Given approximations , , the approximation is obtained after replacing in (LABEL:vcf) by the Lagrange interpolation polynomial of degree , through the points and integrating exactly. Writing

 Pn,k−1(tn+σh)=k−1∑j=0(σj)Δjfn, \hb@xt@.01(2.3)

with , , and the standard forward difference operator, the approximation is computed from

 un+k=ϕ0(k,hA)un+hk−1∑j=0ϕj+1(k,hA)Δjfn, \hb@xt@.01(2.4)

where, for , , and , the mappings are given in (LABEL:phis_multistep). The methods defined in (LABEL:multistep) are explicit and, as we already mentioned in the Introduction, they require the evaluation of .

In [2, Theorem 1] it is shown that if the nonlinearity belongs to the space and the starting values satisfy

 ∥u(tj)−uj∥α≤C0hk,0≤j≤k−1,

the method defined in (LABEL:multistep) exhibits full order , i.e., there exists such that

 ∥u(tn)−un∥α≤Chk∥f(k)∥∞,0≤n≤N. \hb@xt@.01(2.5)

The constant depends on , in (LABEL:Lipschitz), and in (LABEL:sectorial), but it is independent of and .

### 2.2 Exponential Runge–Kutta methods

Explicit exponential Runge–Kutta methods are presented in [10] for the time integration of semilinear parabolic problems. For , , and , the approximations to , with , are given by

 Uni=ecihAun+hi−1∑j=1aij(hA)f(tn+cjh,Unj),un+1=ehAun+hs∑i=1bi(hA)f(tn+cih,Uni), \hb@xt@.01(2.6)

with (). The coefficients and are linear combinations of and with

 φk(λ)=∫10e(1−σ)λσk−1(k−1)!dσ,λ∈C,k≥1,t>0. \hb@xt@.01(2.7)

Setting , we see that the implementation of (LABEL:rk) requires the evaluation of and , for and several values of .

As in [2], the nonlinearity in (LABEL:ivp) is assumed to satisfy a local Lipschitz condition like (LABEL:Lipschitz) and the solution and are assumed to be sufficiently smooth so that the composition

 f∗:[0,T]×Xα→X:t→f∗(t)=f(t,u(t))

is a sufficiently often differentiable mapping, too. Under these assumptions, stiff order conditions are derived and exponential Runge–Kutta methods of the form (LABEL:rk) are constructed up to order four in [10].

## 3 The numerical inversion of sectorial Laplace transforms

In this section we review the numerical inversion method for sectorial Laplace transforms presented in [14], which is a further development of [13].

For a locally integrable mapping , bounded by

 ∥f(t)∥≤Ctν−1eγt,for some  γ∈R, ν>0, \hb@xt@.01(3.1)

we denote its Laplace transform

 F(z)=L[f](z)=∫∞0e−tzf(t)dt,Rez>γ. \hb@xt@.01(3.2)

When satisfies (LABEL:LTsectorial), the method in [14] allows to approximate the values of from few evaluations of . This is achieved by means of a suitable quadrature rule to discretize the inversion formula

 f(t)=12πi∫ΓeztF(z)dz, \hb@xt@.01(3.3)

where is a contour in the complex plane, running from to and laying in the analyticity region of . Due to (LABEL:LTsectorial), can be taken so that it begins and ends in the half plane . Following [14], in (LABEL:invLT) we choose as the left branch of a hyperbola parameterized by

 R→Γ : x↦T(x)=μ(1−sin(α+ix))+γ, \hb@xt@.01(3.4)

where is a scale parameter, is the shift in (LABEL:LTsectorial), and . Thus, is the left branch of the hyperbola with center at , foci at , , and with asymptotes forming angles with the real axis, so that remains in the sector of analyticity of , . In Figure LABEL:fig:hyp we show the action of the conformal mapping on the real axis.

After parameterizing (LABEL:invLT), the function is approximated by applying the truncated trapezoidal rule to the resulting integral along the real axis, i.e.,

 f(t)=12πi∫ΓetzF(z)dz≈K∑ℓ=−KwℓetzℓF(zℓ), \hb@xt@.01(3.5)

with quadrature weights and nodes given by

 wℓ=−τ2πiT′(ℓτ) ,zℓ=T(ℓτ),−K≤ℓ≤K, \hb@xt@.01(3.6)

and a suitable step length parameter. We notice that the minus sign in the formula for the weights comes from setting the proper orientation in the parametrization of . In case of symmetry, the sum in (LABEL:numinvLT) can be halved to

 f(t)≈Re(K∑ℓ=0w∗ℓetzℓF(zℓ)), \hb@xt@.01(3.7)

with and , . The good behavior of the quadrature formula (LABEL:numinvLT) is due to the good properties of the trapezoidal rule when the integrand can be analytically extended to a horizontal strip around the real axis [17, 18].

During the last years, different choices of contours and parameterizations have been studied for the numerical inversion of sectorial Laplace transforms. Apart from the approach in [13, 14], which is the one we follow, the choice of a hyperbola has been studied in [4, 5, 15, 16, 21]. The choice of as a parabola has been considered recently in [4, 5, 21] and we refer also to Talbot’s method [19, 20] for another kind of integration contour , with horizontal asymptotes as .

The error analysis for (LABEL:numinvLT) shows that the role played by the round-off errors is very important. Several ways to avoid error propagation are proposed in [14], depending on the available information about the errors in the evaluations of and the elementary mappings involved in (LABEL:numinvLT). Our algorithm for exponential methods uses (LABEL:quadgeneral) to approximate the required operators in (LABEL:Phiexponencial) and (LABEL:Phigeneral), i.e., (LABEL:numinvLT) with Laplace transforms of the form (LABEL:LTPhi_general). Thus, the required Laplace transforms will be in practice evaluated by performing the inversion of few matrices of the form , with being normally a discrete version of the operator in (LABEL:ivp) and the time step length. Taking into account that this is likely to be accomplished by means of some auxiliary routine, we propose to use in this context the least accurate version of the method in [14], where the error propagation is avoided without using any information about the errors in the evaluations of (LABEL:numinvLT).

The following result provides an error bound for the proposed version of the inversion method like , with the number of quadrature nodes.

###### Theorem 3.1

[14] Assume that the Laplace transform satisfies the sectorial condition (LABEL:LTsectorial) and let and be such that

 0<α−d<α+d<π2−δ. \hb@xt@.01(3.8)

For , and , we select the parameters

 τ=a(K)K,μ=2πdΛt0a(K), \hb@xt@.01(3.9)

with .

Then, the error in the approximation (LABEL:numinvLT) to with quadrature weights and nodes in (LABEL:pesos_nodos) is bounded by

 ∥EK(t)∥≤MΠQe2πd/a(K)tν−1(ε+e−2πdK/a(K)1−e−2πdK/a(K)), \hb@xt@.01(3.10)

uniformly for , where and are the constants in (LABEL:LTsectorial), is the precision in the evaluations of the Laplace transform and the elementary operations in (LABEL:numinvLT),

 Π=1π√1+sin(α+d)(1−sin(α+d))2ν−1

and

 Q=max{4L(λt0sin(α−d)),τ+L(λt0sinα)},

with .

In case we have some reliable information about the errors in the computation of the matrices , an -depending selection of and improves the error bound (LABEL:cota) to , i.e., a pure exponentially decaying bound in . In this situation, given , , and fulfilling (LABEL:angulos), the parameters and are given by

 τ=a(θ∗)K,μ=2πdK(1−θ∗)Λt0a(θ∗), \hb@xt@.01(3.11)

where, for , is the mapping

 a(θ)=arccosh(Λ(1−θ)sinα), \hb@xt@.01(3.12)

and

 θ∗=minθ∈(0,1)(εe2πdK(1−θ)/a(θ)+e−2πdKθ/a(θ)). \hb@xt@.01(3.13)

Given , the expression to be minimize in (LABEL:theta_opt) represents the main part of the error bound obtained in [14] for any fixed . By choosing for every the optimal value in (LABEL:theta_opt), an error bound like is attained. We notice that the error bound stated in Theorem LABEL:thm:err is obtained for in (LABEL:para_eps), instead of in (LABEL:theta_opt). We refer to [14] for details. There, the case in (LABEL:LTsectorial) is also studied.

## 4 Evaluation of the vector-valued mappings

In this section we apply some Laplace transformation formulas to obtain a suitable representation of the operators , , and required in (LABEL:multistep) and (LABEL:rk).

Let us denote the Laplace transform of a mapping by , and the inverse Laplace transform by .

### 4.1 Evaluation of the mappings required by the multistep methods

For in (LABEL:phis_multistep) with , it holds

 ϕj(k,λ)=∫k0e(k−σ)λ(σj)dσ=L−1[L[f0(⋅,λ)]×L[fj]](k),

where, for ,

 f0(σ,λ)=eσλandfj(σ)=(σj). \hb@xt@.01(4.1)

For every and , we define

 Φj(z,λ)=L[f0(⋅,λ)](z)×L[fj](z)=1z−λ×L[fj](z). \hb@xt@.01(4.2)

Then, for every and ,

 ϕj(k,λ)=L−1[Φj(⋅,λ)](k). \hb@xt@.01(4.3)

For

 ϕ0(k,λ)=ekλ=L−1(1⋅−λ)(k),

and thus we define

 Φ0(z,λ)=1z−λ. \hb@xt@.01(4.4)

For scalar, the mappings , with are given by

 Φ1(z,λ)=1z(z−λ),Φ2(z,λ)=1z2(z−λ),Φ3(z,λ)=2−z2z3(z−λ),Φ4(z,λ)=3−3z+z23z4(z−λ). \hb@xt@.01(4.5)

In order to evaluate , , we propose to use the formulas in (LABEL:scalarformulae) with instead of and perform the inversion of the Laplace transform to approximate the original mappings at . In this way, the Laplace transforms we need to invert are:

 Φ0(z,hA)=(zI−hA)−1,Φ1(z,hA)=1z(zI−hA)−1,Φ2(z,hA)=1z2(zI−hA)−1,Φ3(z,hA)=2−z2z3(zI−hA)−1,Φ4(z,hA)=3−3z+z23z4(zI−hA)−1. \hb@xt@.01(4.6)

Although the formulas in (LABEL:operator_multistep) are derived just formally, we notice that they can be justified by combining the Cauchy integral formula with the inversion formula for the Laplace transform. More precisely, for suitable contours and in the complex plane, both laying in the resolvent set of , it holds

 ϕj(k,hA)=12πi∫Γ1ϕj(k,λ)(λI−hA)−1dλ=12πi∫Γ1(12πi∫Γ2ekξΦj(ξ,λ)dξ)(λI−hA)−1dλ=12πi∫Γ2ekξ(12πi∫Γ1Φj(ξ,λ)(λI−hA)−1dλ)dξ=12πi∫Γ2ekξΦj(ξ,hA)dξ=L−1[Φj(⋅,hA)](k). \hb@xt@.01(4.7)

Due to (LABEL:sectorial), all the Laplace transforms in (LABEL:operator_multistep) are sectorial, since they satisfy (LABEL:LTsectorial) with and , for and in (LABEL:sectorial). We notice that, for all , the resulting bounds in (LABEL:LTsectorial) are independent of .

Thus, we can compute the operators , , by using the method described in Section LABEL:sec_LT to compute the inverse Laplace transforms of the mappings in (LABEL:LTPhi_j). We notice that the inverse Laplace transforms need to be approximated only at the fix value , which is specially favorable for the application of the inversion method (see the bound in Theorem LABEL:thm:err). Then, we set , , and select the parameters and following Theorem LABEL:thm:err. The selection of and is more heuristic and a good choice is and slighly smaller than . For example, if in (LABEL:LTsectorial), good values are around and . Next, we compute the quadrature weights and nodes in (LABEL:pesos_nodos) and approximate the operators in (LABEL:multistep) by

 ϕj(k,hA)≈K∑ℓ=−KwℓekzℓΦj(zℓ,hA). \hb@xt@.01(4.8)

The sum in (LABEL:aproxphis) can be halved in case of symmetry like in (LABEL:numinvLT_real).

As we already mentioned in the Introduction, the computation of all the required operators , , can be carried out before the time stepping of the exponential method begins. Thus, if we use the method of lines and apply the exponential method to some spatial discretization of (LABEL:ivp), only the matrix-vector products in (LABEL:multistep) need to be computed at every time step.

### 4.2 Evaluation of the mappings required by the Runge–Kutta methods

For in (LABEL:phis_rk), and , we have

 φj(λ)=∫10e(1−σ)λσj−1(j−1)!dσ=L−1[L[g0(⋅,λ)]×L[gj]](1),

where, for and ,

 g0(σ,λ)=eσλandgj(σ)=σj−1(j−1)!. \hb@xt@.01(4.9)

For every and , we define

 Ψj(z,λ)=L[g0(⋅,λ)](z)×L[gj](z)=1zj(z−λ) \hb@xt@.01(4.10)

and . Then, for every and ,

 φj(λ)=L−1[Ψj(⋅,λ)](1). \hb@xt@.01(4.11)

The same argument as in (LABEL:proofformula) justifies the computation of the operators and , , , by performing the inversion of the Laplace transforms

 Ψj(z,βhA)=1zj(zI−βhA)−1,j≥0,β=1, cl, \hb@xt@.01(4.12)

to approximate the original mappings at . If (LABEL:sectorial) holds, the Laplace transforms in (LABEL:operator_rk) are also the sectorial in the sense of (LABEL:LTsectorial) and we can use the inversion method of [14].

As in the case of the methods in (LABEL:multistep), the computation of all the required operators and , , can be carried out before the time stepping of the exponential method begins.

###### Remark 1

In general, we can always evaluate a mapping of the form of (LABEL:Phigeneral) by using the numerical inversion of the Laplace transform, just by noticing that is the inverse Laplace transform at of a mapping like in (LABEL:LTPhi_general),

 Φ(z,hA)=R(z)(zI−hA)−1,

with , a scalar rational function of .

The above Remark implies that our algorithm can be used to implement other kinds of exponential methods, different than those in [2, 10], as long as they require the evaluation of mappings of the form of (LABEL:Phiexponencial) and (LABEL:Phigeneral).

## 5 Evaluation of the scalar mappings

As we already mentioned in the Introduction, we can also apply the inversion of the Laplace transform to evaluate with accuracy the scalar mappings in (LABEL:Phigeneral). In this section we consider with some detail the evaluation of the mappings

 gj(m,λ)=∫m0e(m−σ)λσj−1dσ,j≥1,m∈N, \hb@xt@.01(5.1)

by means of the quadrature formula (LABEL:quadscalar). The result provided by (LABEL:quadscalar) does not depend on the size of , but the formula is only useful in principle for values of inside a sector of the form . However, using that

 emλg1(m,−λ)=g1(m,λ),λ∈C, \hb@xt@.01(5.2)

and

 gj+1(m,λ)=jgj(m,λ)−mjλ,j≥1,m∈N, \hb@xt@.01(5.3)

it is easy to see by induction that, for and ,

 emλgj(m,−λ)=j∑ℓ=1(j−1ℓ−1)(−1)ℓ−1mj−ℓgℓ(m,λ),j≥1. \hb@xt@.01(5.4)

Thus, we can compute

 gj(m,−λ)=e−mλL−1[G∗j(⋅,λ)](m), \hb@xt@.01(5.5)

with

 \hb@xt@.01(5.6)

which provides a stable formula to approximate for inside a proper sector and moderate size.

In Table LABEL:tab1 we show the error obtained in the evaluation of in (LABEL:vfi_1) for different values of in the interval . For , we applied the inversion formula (LABEL:numinvLT) with and

 F(z)=G1(z,λ)=1z(z−λ), \hb@xt@.01(5.7)

which, for these values of , fulfils (LABEL:LTsectorial) with , , and . We assumed that the evaluations of can be carried out in MATLAB up to machine accuracy and thus we set . Then, we computed the quadrature weights and nodes in (LABEL:pesos_nodos) following (LABEL:para_eps)–(LABEL:theta_opt) with . Setting and , we obtained , for , and , for . In Table LABEL:tab1 we can see that is enough to attain almost the machine accuracy of MATLAB in the evaluations of . For positive values of , we used (LABEL:g1_minus) with .

## 6 Numerical illustrations

In this section we test our algorithm by considering the same examples as in [2] and [10].

### 6.1 Example for the multistep exponential methods

Our first example is the problem considered in [2]

 ut(x,t)=uxx(x,t)+(∫10u(s,t)ds)ux(x,t)+g(x,t), \hb@xt@.01(6.1)

for and , subject to homogeneous Dirichlet boundary conditions and with such that the exact solution to (LABEL:ex1_cp) is

The spatial discretization of (LABEL:ex1_cp) is carried out by using standard finite differences with spatial nodes, centered for the approximation of . The nonlocal term is approximated by means of the composite Simpson’s formula.

To integrate in time the semidiscrete problem we use (LABEL:multistep) with and , so that is the matrix

 A=J2tridiag([1,−2,1]).

We approximate the matrices , , required in (LABEL:multistep) by applying the quadrature rule (LABEL:aproxphis). To avoid an extra source of error, the initial values are computed from the exact solution. In a less academic example, these values can be computed by means of a one-step method of sufficiently high order or by performing the fix point iteration proposed in [2].

In Figure LABEL:fig:ex1_cp we show the error versus the stepsize at , measured in a discrete version of the norm , for and in (LABEL:aproxphis). We see that for the full precision is achieved for all the methods implemented; cf. [2, Section 6]. In Figure LABEL:fig:ex1_cp we also show lines of slope 1, 2, 3 and 4, to visualize the order of convergence. In fact, we can observe that the order of convergence for this example is slightly higher than the one predicted in [2], approximately 2.15, 3.15 and 4.15 for and , respectively, instead of sharp order . A further study of this behavior is beyond the scope of this paper.

### 6.2 Examples for the exponential Runge–Kutta methods

For the following two examples we consider the problems and some of the integrators proposed in [10]. Following the notation in [10], in the Butcher tableaus below we use the abbreviations

 φi=φi(hA),and φi,j=φi,j(hA)=φi(cjhA),2≤j≤s. \hb@xt@.01(6.2)

Thus, our second example is

 ut(x,t)=uxx(x,t)+11+u(x,t)2+g(x,t), \hb@xt@.01(6.3)

for and , subject to homogeneous Dirichlet boundary conditions and with such that the exact solution to (LABEL:ex1_ho) is again .

We discretize this problem in space by standard finite differences with grid points. For the time integration of the semidiscrete problem, we implemented (LABEL:rk) with , the second-order method

 01212φ1,20φ1 \hb@xt@.01(6.4)

the third-order method

 01313φ1,22323φ1,3−43φ2,343φ2,3φ1−32φ2032φ2 \hb@xt@.01(6.5)

and the fourth-order one

 01212φ1,21212φ1,3−φ2,3φ2,31φ1,4−2φ2,4φ2,4φ2,41212φ1,5−2a5,2−a5,4a5,2a5,2a5,4φ1−3φ2+4φ300−φ2+4φ34φ2−8φ3 \hb@xt@.01(6.6)

with

 a5,2=12φ2,5−φ3,4+14φ2,4−12φ3,5

and

 a5,4=14φ2,5−a5,2.

For the implementation of (LABEL:rk2) we need to invert four different Laplace transforms of the form of (LABEL:operator_rk), to approximate , , and . The implementation of both (LABEL:rk3) and (LABEL:rk4) requires the inversion of eight Laplace transforms.

In Figure LABEL:fig:ex1_ho we show the error at versus the stepsize, measured in the maximum norm. The expected order of convergence for this example is for the -order method. In order to check our algorithm, we added lines with the corresponding slopes in Figure LABEL:fig:ex1_ho. We can see that also for this kind of methods we attain full precision for in (LABEL:aproxphis).

Finally we consider the second example in [10]

 ut(x,t)=uxx(x,t)+∫10u(s,t)ds+g(x,t), \hb@xt@.01(6.7)

for and , subject also to homogeneous Dirichlet boundary conditions and with such that the exact solution to (LABEL:ex2_ho) is .

We discretize this problem in space as in the previous example (LABEL:ex1_ho), and use the composite Simpson’s rule for the approximation of the nonlocal term. For the time integration, we use (LABEL:rk) with