# On the implementation of exponential methods for semilinear parabolic equations

###### 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

\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

\hb@xt@.01(1.2) |

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

\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

\hb@xt@.01(1.4) |

or

\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

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

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

\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

\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

\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

\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

\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

\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

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

\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

\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

\hb@xt@.01(2.3) |

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

\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

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

\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

\hb@xt@.01(2.6) |

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

\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

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

\hb@xt@.01(3.1) |

we denote its Laplace transform

\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

\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

\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.,

\hb@xt@.01(3.5) |

with quadrature weights and nodes given by

\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

\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

\hb@xt@.01(3.8) |

For , and , we select the parameters

\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

\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),

and

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

\hb@xt@.01(3.11) |

where, for , is the mapping

\hb@xt@.01(3.12) |

and

\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

where, for ,

\hb@xt@.01(4.1) |

For every and , we define

\hb@xt@.01(4.2) |

Then, for every and ,

\hb@xt@.01(4.3) |

For

and thus we define

\hb@xt@.01(4.4) |

For scalar, the mappings , with are given by

\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:

\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

\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

\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

where, for and ,

\hb@xt@.01(4.9) |

For every and , we define

\hb@xt@.01(4.10) |

and . Then, for every and ,

\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

\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),

with , a scalar rational function of .

## 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

\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

\hb@xt@.01(5.2) |

and

\hb@xt@.01(5.3) |

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

\hb@xt@.01(5.4) |

Thus, we can compute

\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

\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 .

-1 | 1.5050e-12 | 1.3323e-15 | 1 | 1.5050e-12 | 3.3307e-15 |

-1e-1 | 1.5227e-12 | 3.2196e-15 | 1e-1 | 1.5227e-12 | 3.5527e-15 |

-1e-2 | 1.4243e-12 | 4.4409e-15 | 1e-2 | 1.4243e-12 | 4.6629e-15 |

-1e-3 | 1.3750e-12 | 1.3323e-15 | 1e-3 | 1.3750e-12 | 1.3323e-15 |

-1e-4 | 1.3738e-12 | 1.7764e-15 | 1e-4 | 1.3738e-12 | 1.7764e-15 |

-1e-5 | 1.3747e-12 | 3.6637e-15 | 1e-5 | 1.3747e-12 | 3.7748e-15 |

-1e-6 | 1.3748e-12 | 3.6637e-15 | 1e-6 | 1.3748e-12 | 3.7748e-15 |

-1e-7 | 1.3695e-12 | 1.9984e-15 | 1e-7 | 1.3695e-12 | 1.9984e-15 |

-1e-8 | 1.3717e-12 | 1.1102e-16 | 1e-8 | 1.3717e-12 | 2.2204e-16 |

-1e-9 | 1.3715e-12 | 1.1102e-16 | 1e-9 | 1.3715e-12 | 0 |

-1e-10 | 1.3711e-12 | 0 | 1e-10 | 1.3711e-12 | 0 |

-1e-11 | 1.3711e-12 | 0 | 1e-11 | 1.3711e-12 | 0 |

-1e-12 | 1.3715e-12 | 1.1102e-16 | 1e-12 | 1.3715e-12 | 0 |

-1e-13 | 1.3712e-12 | 0 | 1e-13 | 1.3712e-12 | 2.2204e-16 |

## 6 Numerical illustrations

### 6.1 Example for the multistep exponential methods

Our first example is the problem considered in [2]

\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

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

\hb@xt@.01(6.2) |

Thus, our second example is

\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

\hb@xt@.01(6.4) |

the third-order method

\hb@xt@.01(6.5) |

and the fourth-order one

\hb@xt@.01(6.6) |

with

and

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]

\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