Exponential quadrature rules without order reduction

Due to the recent development and improvement of Krylov methods [7, 11], exponential quadrature rules have become a valuable tool to integrate linear initial value partial differential problems [9]. This is because of the fact that the linear and stiff part of the problem can be ‘exactly’ integrated in an efficient way through exponential operators. Moreover, when the source term is nontrivial, a variations-of-constants formula and the interpolation of that source term in several nodes leads to the appearance of some $\varphi_{j}$-operators, to which Krylov techniques can also be applied.

However, in the literature [9], these methods have always been applied and analysed either on the initial value problem or on the initial boundary value one with vanishing or periodic boundary conditions. More precisely, when the linear and stiff operator is the infinitesimal generator of a strongly continuous semigroup in a certain Banach space. In such a case, it has been proved [9] that the exponential quadrature rule converges with global order $s$ if $s$ is the number of nodes being used for the interpolation of the source term.

There are no results concerning specifically how to deal with these methods when integrating the most common non-vanishing and time-dependent boundary conditions case. The only related reference is that of Lawson quadrature rules [3, 4, 10] which differ from these methods in the fact that, not only the source term is interpolated, but all the integrand which turns up in the variations-of-constants formula. In such a way, with the latter methods, $\{\varphi_{j}\}$-operators (with $j\geq 1$) do not turn up. Only exponential functions (those corresponding to $j=0$) are present. In [3] a thorough error analysis is given which studies the strong order reduction which turns up with Lawson methods even in the vanishing boundary conditions case unless some even more artificial additional vanishing boundary conditions are satisfied. Nevertheless, in [4], a technique is suggested to avoid that order reduction under homogeneous boundary conditions and to even tackle the time-dependent boundary case without order reduction. Moreover, the analysis there also includes the error coming from the space discretization.

The aim of this paper is to generalize that technique to the most common quadrature rules which are used in the literature and which also use $\{\varphi_{j}\}$-operators. Besides, we also include the error coming from the space discretization not only when avoiding order reduction but also with the classical approach. We will see that, with the technique which is suggested in this paper, we manage to get the order of the classical quadrature interpolatory rule, which implies that, by choosing the $s$ nodes carefully, we can manage to get even order $2s$ in time.

The paper is structured as follows. Section 2 gives some preliminaries on the abstract framework in Banach spaces which is used for the time integration of the problem, on the definition of the exponential quadrature rules and on the general hypotheses which are used for the abstract space discretization. Then, Section 3 describes and makes a thorough error analysis of the classical method of lines, which integrates the problem firstly in space and then in time. Both vanishing and nonvanishing boundary conditions are considered there and several different results are obtained depending on the specific accuracy of the quadrature rule and on whether a parabolic assumption is satisfied. After that, the technique which is suggested in the paper to improve the order of accuracy in time is well described in Section 4. It consists of discretizing firstly in time with suitable boundary conditions and then in space. Therefore, the analysis is firstly performed on the local error of the time semidiscretization and then on the local and global error of the full scheme (27),(32),(33). Finally, Section 5 shows some numerical results which corroborate the theoretical results of the previous sections.

As in [4], we consider the linear abstract initial boundary value problem in the complex Banach space $X$

$\displaystyle\begin{array}[]{rcl}u^{\prime}(t)&=&Au(t)+f(t),\quad 0\leq t\leq T,\\ u(0)&=&u_{0},\\ \partial u(t)&=&g(t),\end{array}$

(4)

where $D(A)$ is a dense subspace of $X$ and the linear operators $A$ and $\partial$ satisfy the following assumptions:

1. (A1)

   The boundary operator $\partial:D(A)\subset X\to Y$ is onto, where $Y$ is another Banach space.

2. (A2)

   Ker($\partial$) is dense in $X$ and $A_{0}:D(A_{0})=Ker(\partial)\subset X\to X$, the restriction of $A$ to Ker($\partial$), is the infinitesimal generator of a $C_{0}$- semigroup $\{e^{tA_{0}}\}_{t\geq 0}$ in $X$, which type is denoted by $\omega$, and we assume that it is negative.

3. (A3)

   If $z\in\mathbb{C}$ satisfies $Re(z)>\omega$ and $v\in Y$, then the steady state problem

   	$\displaystyle Ax$	$\displaystyle=$	$\displaystyle zx,$
   	$\displaystyle\partial x$	$\displaystyle=$	$\displaystyle v,$

   possesses a unique solution denoted by $x=K(z)v$. Moreover, the linear operator $K(z):Y\to D(A)$ satisfies

   $\displaystyle\|K(z)v\|_{X}\leq L\|v\|_{Y},$

   where the constant $L$ holds for any $z$ such that $Re(z)>\omega_{0}>\omega$.

As precisely stated in [4, 12], with these assumptions, problem (4) is well-posed in $BV/L^{\infty}$ sense. Moreover, as $\omega$ is negative, the semigroup $e^{tA_{0}}$ decays exponentially when $t\to 0$ and the operator $A_{0}$ is invertible.

We also remind that, because of hypothesis (A2), $\{\varphi_{j}(tA_{0})\}_{j=0}^{\infty}$ are bounded operators for $t>0$, where $\{\varphi_{j}\}$ are the standard functions which are used in exponential methods [9]:

$\displaystyle\varphi_{j}(tA_{0})=\frac{1}{t^{j}}\int_{0}^{t}e^{(t-\tau)A_{0}}\frac{\tau^{j-1}}{(j-1)!}d\tau.$

(5)

It is well-known that they can be calculated in a recursive way through the formula

$\displaystyle\varphi_{k+1}(z)=\frac{\varphi_{k}(z)-1/k!}{z},\quad z\neq 0,\qquad\varphi_{k+1}(0)=\frac{1}{(k+1)!},\qquad\varphi_{0}(z)=e^{z},$

(6)

and that these functions are bounded on the complex plane when $\mbox{Re}(z)\leq 0$.

We also assume that the solution of (4) satisfies that, for a natural number $p$,

$\displaystyle A^{p+1-j}u^{(j)}\in C([0,T],X),\quad 0\leq j\leq p+1.$

(7)

When $A$ is a differential space operator, this assumption implies that the time derivatives of the solution are regular in space, but without imposing any restriction on the boundary. Because of Theorem 3.1 in [1], this assumption is satisfied when the data $u_{0}$, $f$ and $g$ are regular and satisfy certain natural compatibility conditions at the boundary. Moreover, as from (4),

$\displaystyle u^{(j)}(t)=\sum_{l=0}^{j-1}A^{l}f^{(j-1-l)}(t)+A^{j}u(t),\quad 0\leq j\leq p+1,$

(8)

the following crucial boundary values

$\displaystyle\partial A^{j}u(t)=g^{(j)}(t)-\sum_{l=0}^{j-1}\partial A^{l}f^{(j-1-l)}(t),\quad 0\leq j\leq p+1,$

can be calculated from the given data in problem (4).

We will center on exponential quadrature rules [9] to time integrate (4). When applied to a finite-dimensional linear problem like

$\displaystyle U^{\prime}(t)=BU(t)+F(t),$

(9)

where $B$ is a matrix, these rules correspond to interpolating $F$ in $s$ nodes $\{c_{i}\}_{i=1}^{s}$ in the integral in the equality

$\displaystyle U(t_{n+1})=e^{kB}U(t_{n})+k\int_{0}^{1}e^{k(1-\theta)B}F(t_{n}+\theta k)d\theta,$

(10)

which is satisfied by the solutions of (9) when $t_{n+1}=t_{n}+k$. This yields

$\displaystyle U^{n+1}=e^{kB}U^{n}+k\sum_{i=1}^{s}b_{i}(kB)F(t_{n}+c_{i}k),$

with weights

$\displaystyle b_{i}(kB)=\int_{0}^{1}e^{k(1-\theta)B}l_{i}(\theta)d\theta,$

where $l_{i}$ are the Lagrange interpolation polynomials corresponding to the nodes $\{c_{i}\}_{i=1}^{s}$. We will define the values $\{a_{ij}\}_{i,j=1}^{s}$ in such a way that

$\displaystyle l_{i}(\theta)=a_{i,1}+a_{i,2}\theta+a_{i,3}\frac{\theta^{2}}{2}+\dots+a_{i,s}\frac{\theta^{s-1}}{(s-1)!}.$

(11)

From this,

$$b_{i}(\tau B)=\int_{0}^{1}e^{\tau(1-\theta)B}l_{i}(\theta)d\theta=\frac{1}{\tau}\int_{0}^{\tau}e^{(\tau-\sigma)B}l_{i}(\frac{\sigma}{\tau})d\sigma=\sum_{j=1}^{s}a_{i,j}\varphi_{j}(\tau B),$$

for the functions $\varphi_{j}$ in (5), and the final formula for the integration of (9) is

$\displaystyle U^{n+1}=e^{kB}U^{n}+k\sum_{i,j=1}^{s}a_{i,j}\varphi_{j}(kB)F(t_{n}+c_{i}k).$

(12)

We will consider an abstract spatial discretization which satisfies the same hypotheses as in [4] (Section 4.1) and which includes a big range of techniques. In such a way, for each parameter $h$ in a sequence $\{h_{j}\}_{j=1}^{\infty}$ such that $h_{j}\to 0$, $X_{h}\subset X$ is a finite dimensional space which approximates $X$ when $h_{j}\to 0$ and the elements in $D(A_{0})$ are approximated in a subspace $X_{h,0}$. The norm in $X_{h}$ is denoted by $\|\cdot\|_{h}$. The operator $A$ is then approximated by $A_{h}$, $A_{0}$ by $A_{h,0}$ and the solution of the elliptic problem

$$Aw=F,\qquad\partial w=g,$$

is approximated by $R_{h}w+Q_{h}g$, where $R_{h}w\in X_{h,0}$ is called the elliptic projection, $Q_{h}g\in X_{h}$ discretizes the boundary values and the following is satisfied:

$\displaystyle A_{h,0}R_{h}w+A_{h}Q_{h}g=L_{h}F,$

(13)

for a projection operator $L_{h}:X\to X_{h,0}$. We will also use $P_{h}=L_{h}-L_{h}Q_{h}\partial$ and we remind part of hypothesis (H3) in [4], which states that, for a subspace $Z$ of $X$ with norm $\|\cdot\|_{Z}$, whenever $u\in Z$,

$\displaystyle\|A_{h,0}(R_{h}-P_{h})u\|_{h}\leq\varepsilon_{h}\|u\|_{Z},$

(14)

for $\varepsilon_{h}$ decreasing with $h$ and, therefore, this gives a bound for the error in the space discretization of operator $A$.

Moreover, we will assume that this additional hypothesis is satisfied:

1. (HS)

   $\|A_{h,0}^{-1}A_{h}Q_{h}\|_{h}$ is bounded independently of $h$ for small enough $h$. Considering (13), this in fact corresponds to a discrete maximum principle, which would be simulating the continuous maximum principle which is satisfied because of (A3) when $z=0$.

When considering vanishing boundary conditions in (4) (which has been classically done in the literature with exponential methods [9]), discretizing first in space and then in time leads to the following semidiscrete problem in $X_{h,0}$:

	$\displaystyle U_{h}^{\prime}(t)$	$\displaystyle=$	$\displaystyle A_{h,0}U_{h}(t)+L_{h}f(t),$
	$\displaystyle U_{h}(0)$	$\displaystyle=$	$\displaystyle L_{h}u(0).$

When integrating this problem with an exponential quadrature rule which is based on $s$ nodes (12), the following scheme arises:

$\displaystyle U_{h}^{n+1}=e^{kA_{h,0}}U_{h}^{n}+k\sum_{i,j=1}^{s}a_{i,j}\varphi_{j}(kA_{h,0})L_{h}f(t_{n}+c_{i}k).$

(15)

Denoting by $\rho_{h,n+1}$ to $U_{h}(t_{n+1})-\bar{U}_{h}^{n+1}$ where $\bar{U}_{h}^{n+1}$ is the result of applying (15) from $U_{h}(t_{n})$ instead of $U_{h}^{n}$; and $e_{h,n+1}$ to $U_{h}(t_{n+1})-U_{h}^{n+1}$ the following result follows:

We also have this finer result, which implies global order $s+1$ under more restrictive hypotheses.

On the other hand, when $g\not\equiv 0$ in (4), the semidiscretized problem which arises is

$\displaystyle\left.\begin{array}[]{lcl}U^{\prime}_{h}(t)&=&A_{h,0}U_{h}(t)+A_{h}Q_{h}g(t)+L_{h}Q_{h}(\partial f(t)-g^{\prime}(t))+P_{h}f(t),\\ U_{h}(0)&=&P_{h}u(0).\end{array}\right.$

(20)

In a similar way as before, the local error would be given by

	$\displaystyle k\int_{0}^{1}e^{k(1-\theta)A_{h,0}}\bigg{[}A_{h}Q_{h}[g(t_{n}+k\theta)-I(g(t_{n}+k\theta))]$				(21)
			$\displaystyle\hskip 71.13188pt+L_{h}Q_{h}[\partial f(t_{n}+k\theta)-g^{\prime}(t_{n}+k\theta)-I(\partial f(t_{n}+k\theta)-g^{\prime}(t_{n}+k\theta))$
			$\displaystyle\hskip 71.13188pt+P_{h}[f(t_{n}+k\theta)-I(f(t_{n}+k\theta))]\bigg{]}d\theta.$

Again, when $g\in C^{s+1}([0,T],Y)$ and $f\in C^{s}([0,T],X)$, the error of interpolation will be $O(k^{s})$. However, although $L_{h}Q_{h}$ and $P_{h}$ are bounded [4], $A_{h}Q_{h}$ is not bounded any more. That is why we state the following result which bounds in fact $A_{h,0}^{-1}\rho_{h,n}$ by using (HS) and which proof for the global error is the same as in Theorem 2.

In any case, we want to remark in this section that accuracy has been lost with respect to the vanishing boundary conditions case since order reduction turns up at least for the local error and, in many cases, also for the global error.

In this section, we directly tackle the nonvanishing boundary conditions case by discretizing in a suitable way firstly in time and then in space. We will see that we manage to get at least the same order as with the classical approach when vanishing boundary conditions are present, but even a much higher order some times.

Let us suggest how to apply the exponential quadrature rule (12) directly to (4). When $g=0$, $B$ in (12) is directly substituted by $A_{0}$ and there is no problem because $e^{kA_{0}}$ and $\varphi_{j}(kA_{0})$ have perfect sense over $X$. However, it has no sense to do that when $g\neq 0$ because $A$ is not $A_{0}$ any more. For Lawson methods, for which just exponential functions appear, instead of $e^{\tau A_{0}}\alpha$, it was suggested in [4] to consider $v_{0}(\tau)$ as the solution of

	$\displaystyle v_{0}^{\prime}(\tau)$	$\displaystyle=$	$\displaystyle Av_{0}(\tau),$
	$\displaystyle v_{0}(0)$	$\displaystyle=$	$\displaystyle\alpha,$
	$\displaystyle\partial v_{0}(\tau)$	$\displaystyle=$	$\displaystyle\sum_{l=0}^{p}\frac{\tau^{l}}{l!}\partial A^{l}\alpha,$		(22)

whenever $\alpha\in D(A^{p})$. In such a way, if $\alpha\in D(A^{p+1})$,

$\displaystyle v_{0}(\tau)=\sum_{l=0}^{p}\frac{\tau^{l}}{l!}A^{l}\alpha+\tau^{p+1}\varphi_{p+1}(\tau A_{0})A^{p+1}\alpha,$

(23)

which resembles the formal analytic expansion of the exponential of $\tau A$ applied over $\alpha$. In this manuscript then, whenever $\alpha\in D(A^{p})$, for $j=1,\dots,s$, instead of $\varphi_{j}(\tau A_{0})\alpha$, we suggest to consider the following functions :

$\displaystyle v_{j}(\tau)=\sum_{l=0}^{p-1}\frac{\tau^{l}}{(l+j)!}A^{l}\alpha+\tau^{p}\varphi_{p+j}(\tau A_{0})A^{p}\alpha.$

(24)

This resembles the formal analytic expansion of $\varphi_{j}$ when evaluated at $\tau A$ and applied over $\alpha$. (Notice that, for $j=0$, this would correspond to (23) changing $p$ by $p-1$. As the functions $\varphi_{j}$ are multiplied by $k$ in (12), we need one less term in this expansion.)

Therefore, imitating (12), we suggest to consider as continuous numerical approximation $u_{n+1}$ from the previous $u_{n}$,

	$\displaystyle u_{n+1}=\hbox to0.0pt{$\displaystyle\sum_{l=0}^{p}\frac{k^{l}}{l!}A^{l}u_{n}+k^{p+1}\varphi_{p+1}(kA_{0})A^{p+1}u_{n}$\hss}$				(25)
			$\displaystyle+k\sum_{i,j=1}^{s}a_{i,j}\bigg{[}\sum_{l=0}^{p-1}\frac{k^{l}}{(l+j)!}A^{l}f(t_{n}+c_{i}k)+k^{p}\varphi_{p+j}(kA_{0})A^{p}f(t_{n}+c_{i}k)\bigg{]}.$

It is not practical to discretize directly this formula in space since numerical differentiation is a badly-posed problem. However, if we seek a differential equation which the functions (24) satisfy and apply a numerical integrator to it, there should be no problems. For that, let us first consider the following lemma.

From here, the next result follows:

With Lawson methods [4], starting from a previous approximation $U_{h,n}$ to $P_{h}u(t_{n})$ and discretizing (22) in space with $\alpha=u(t_{n})$ in $\partial v_{0}(\tau)$ led to a term like

	$\displaystyle V_{h,n,0}(k)$	$\displaystyle=$	$\displaystyle e^{kA_{h,0}}U_{h,n}+\sum_{j=1}^{p}k^{j}\varphi_{j}(kA_{h,0})[A_{h}Q_{h}\partial A^{j-1}u(t_{n})-L_{h}Q_{h}\partial A^{j}u(t_{n})]$		(27)
			$\displaystyle+k^{p+1}\varphi_{p+1}(kA_{h,0})A_{h}Q_{h}\partial A^{p}u(t_{n}),$

which is the approximation which corresponds to the first term in (12).

For the rest of the terms in (12), we suggest to discretize (26) in space with $\alpha=f(t_{n}+c_{i}k)$. In such a way, the following system turns up:

	$\displaystyle V_{h,j,n,i}^{\prime}(\tau)+L_{h}Q_{h}\partial\hat{v}_{j,n,i}^{\prime}(\tau)$	$\displaystyle=$	$\displaystyle A_{h,0}V_{h,j,n,i}(\tau)+A_{h}Q_{h}\partial\hat{v}_{j,n,i}(\tau)-\frac{j}{\tau}[V_{h,j,n,i}(\tau)+L_{h}Q_{h}\partial\hat{v}_{j,n,i}(\tau)]$
			$\displaystyle+\frac{1}{(j-1)!\tau}[P_{h}f(t_{n}+c_{i}k)+L_{h}Q_{h}\partial f(t_{n}+c_{i}k)],$
	$\displaystyle V_{h,j,n,i}(0)+L_{h}Q_{h}\partial\hat{v}_{j,n,i}(0)$	$\displaystyle=$	$\displaystyle\frac{1}{j!}L_{h}f(t_{n}+c_{i}k),$

where

$\displaystyle\hat{v}_{j,n,i}(\tau)=\sum_{l=0}^{p-1}\frac{\tau^{l}}{(l+j)!}A^{l}f(t_{n}+c_{i}k).$

This can be rewritten as

	$\displaystyle\hskip 5.69046ptV_{h,j,n,i}^{\prime}(\tau)$	$\displaystyle=$	$\displaystyle(A_{h,0}-\frac{j}{\tau}I)V_{h,j,n,i}(\tau)+A_{h}Q_{h}\partial\hat{v}_{j,n,i}(\tau)+\frac{1}{(j-1)!\tau}P_{h}f(t_{n}+c_{i}k)$
			$\displaystyle+L_{h}Q_{h}\partial[\frac{1}{(j-1)!\tau}f(t_{n}+c_{i}k)-\frac{j}{\tau}\hat{v}_{j,n,i}(\tau)-\hat{v}_{j,n,i}^{\prime}(\tau)],$
	$\displaystyle\hskip 5.69046ptV_{h,j,n,i}(0)$	$\displaystyle=$	$\displaystyle\frac{1}{j!}P_{h}f(t_{n}+c_{i}k).$		(28)

With the same arguments as in Lemma 6, $\varphi_{j}(\tau A_{h,0})P_{h}f(t_{n}+c_{i}k)$ is the solution of

	$\displaystyle W_{h,j,n,i}^{\prime}(\tau)$	$\displaystyle=$	$\displaystyle(A_{h,0}-\frac{j}{\tau}I)W_{h,j}(\tau)+\frac{1}{(j-1)!\tau}P_{h}f(t_{n}+c_{i}k),$
	$\displaystyle W_{h,j,n,i}(0)$	$\displaystyle=$	$\displaystyle\frac{1}{j!}P_{h}f(t_{n}+c_{i}k).$

Therefore, in order to solve (28), we are interested in finding

$\displaystyle Z_{h,j,n,i}(\tau)=V_{h,j,n,i}(\tau)-W_{h,j,n,i}(\tau),$

(29)

which is the solution of

	$\displaystyle Z_{h,j,n,i}^{\prime}(\tau)$	$\displaystyle=$	$\displaystyle(A_{h,0}-\frac{j}{\tau}I)Z_{h,j,n,i}(\tau)+A_{h}Q_{h}\partial\hat{v}_{j,n,i}(\tau)$
			$\displaystyle+L_{h}Q_{h}\partial[\frac{1}{(j-1)!\tau}f(t_{n}+c_{i}k)-\frac{j}{\tau}\hat{v}_{j,n,i}(\tau)-\hat{v}_{j,n,i}^{\prime}(\tau)],$
	$\displaystyle Z_{h,j,n,i}(0)$	$\displaystyle=$	$\displaystyle 0.$		(30)

Now, using the first line of (26) for the boundary with $\alpha=f(t_{n}+c_{i}k)$,

	$\displaystyle\partial[\frac{1}{(j-1)!\tau}f(t_{n}+c_{i}k)-\frac{j}{\tau}\hat{v}_{j,n,i}(\tau)-\hat{v}_{j,n,i}^{\prime}(\tau)]=-\partial Av_{j,n,i}(\tau)$
		$\displaystyle=$	$\displaystyle-\sum_{l=0}^{p-1}\frac{\tau^{l}}{(l+j)!}\partial A^{l+1}f(t_{n}+c_{i}k)-\tau^{p}\partial A_{0}\varphi_{p+j}(\tau A_{0})A^{p}f(t_{n}+c_{i}k),$

the fact that

	$\displaystyle\tau^{p}A_{0}\frac{1}{\tau^{p+j}}\int_{0}^{\tau}e^{(\tau-\sigma)A_{0}}\frac{\sigma^{p+j-1}}{(p+j-1)!}A^{p}f(t_{n}+c_{i}k)d\sigma$
			$\displaystyle=-\frac{1}{\tau^{j}}e^{(\tau-\sigma)A_{0}}\frac{\sigma^{p+j-1}}{(p+j-1)!}|_{\sigma=0}^{\sigma=\tau}A^{p}f(t_{n}+c_{i}k)+\frac{1}{\tau^{j}}\int_{0}^{\tau}e^{(\tau-\sigma)A_{0}}\frac{\sigma^{p+j-2}}{(p+j-2)!}A^{p}f(t_{n}+c_{i}k)d\sigma$
			$\displaystyle=-\frac{\tau^{p-1}}{(p+j-1)!}A^{p}f(t_{n}+c_{i}k)+\tau^{p-1}\varphi_{p+j-1}(\tau A_{0})A^{p}f(t_{n}+c_{i}k),$

and that the boundary of the second term vanishes, it follows that

$\displaystyle\partial[\frac{1}{(j-1)!\tau}f(t_{n}+c_{i}k)-\frac{j}{\tau}\hat{v}_{j,n,i}(\tau)-\hat{v}_{j,n,i}^{\prime}(\tau)]=-\sum_{l=0}^{p-2}\frac{\tau^{l}}{(l+j)!}\partial A^{l+1}f(t_{n}+c_{i}k).$

Using this in (30),

	$\displaystyle Z_{h,j,n,i}(\tau)$				(31)
		$\displaystyle=$	$\displaystyle\int_{0}^{\tau}e^{\int_{\theta}^{\tau}(A_{h,0}-\frac{j}{\sigma}I)d\sigma}\sum_{l=0}^{p-2}\frac{\theta^{l}}{(l+j)!}\bigg{[}A_{h}Q_{h}\partial A^{l}f(t_{n}+c_{i}k)-L_{h}Q_{h}\partial A^{l+1}f(t_{n}+c_{i}k)$
			$\displaystyle\hskip 170.71652pt+\frac{\theta^{p-1}}{(p-1+j)!}A_{h}Q_{h}\partial A^{p-1}f(t_{n}+c_{i}k)\bigg{]}d\theta$
		$\displaystyle=$	$\displaystyle\sum_{l=0}^{p-2}\int_{0}^{\tau}e^{A_{h,0}(\tau-\theta)}\frac{\theta^{j+l}}{\tau^{j}(l+j)!}[A_{h}Q_{h}\partial A^{l}f(t_{n}+c_{i}k)-L_{h}Q_{h}\partial A^{l+1}f(t_{n}+c_{i}k)]d\theta$
			$\displaystyle+\int_{0}^{\tau}e^{A_{h,0}(\tau-\theta)}\frac{\theta^{p-1+j}}{\tau^{j}(p-1+j)!}A_{h}Q_{h}\partial A^{p-1}f(t_{n}+c_{i}k)d\theta$
		$\displaystyle=$	$\displaystyle\sum_{l=0}^{p-2}\tau^{l+1}\varphi_{j+l+1}(\tau A_{h,0})[A_{h}Q_{h}\partial A^{l}f(t_{n}+c_{i}k)-L_{h}Q_{h}\partial A^{l+1}f(t_{n}+c_{i}k)]$
			$\displaystyle+\tau^{p}\varphi_{p+j}(\tau A_{h,0})A_{h}Q_{h}\partial A^{p-1}f(t_{n}+c_{i}k).$

Therefore, using (29),

	$\displaystyle\hskip-28.45274ptV_{h,j,n,i}(k)=\varphi_{j}(kA_{h,0})P_{h}f(t_{n}+c_{i}k)$				(32)
			$\displaystyle+\sum_{l=0}^{p-2}k^{l+1}\varphi_{j+l+1}(kA_{h,0})[A_{h}Q_{h}\partial A^{l}f(t_{n}+c_{i}k)-L_{h}Q_{h}\partial A^{l+1}f(t_{n}+c_{i}k)],$
			$\displaystyle+k^{p}\varphi_{p+j}(kA_{h,0})A_{h}Q_{h}\partial A^{p-1}f(t_{n}+c_{i}k),$

and the overall exponential quadrature rule would be given by

	$\displaystyle U_{h,0}$	$\displaystyle=$	$\displaystyle P_{h}u_{0},$
	$\displaystyle U_{h,n+1}$	$\displaystyle=$	$\displaystyle V_{h,n,0}(k)+k\sum_{i,j=1}^{s}a_{i,j}V_{h,j,n,i}(k),$		(33)

with $V_{h,n,0}(k)$ in (27) and $V_{h,j,n,i}(k)$ in (32).