# Avoiding order reduction with explicit Runge-Kutta exponential methods in
nonlinear initial boundary value problems

B. Cano
IMUVA, Departamento de Matemática Aplicada,
Facultad de Ciencias, Universidad de Valladolid,
Paseo de Belén 7, 47011 Valladolid,
Spain
and
M. J. Moreta
IMUVA, Departamento de Análisis Económico,
Facultad de Ciencias Económicas y Empresariales, Universidad Complutense de
Madrid,
Campus de Somosaguas, Pozuelo de Alarcón, 28223 Madrid
Spain Corresponding author. Email: bcano<EMAIL_ADDRESS>
###### Abstract

In this paper a technique is given to recover the classical order of the
method when explicit exponential Runge-Kutta methods integrate reaction-
diffusion problems. Although methods of high stiff order for problems with
vanishing boundary conditions can be constructed, that may imply increasing
the number of stages and therefore, the computational cost seems bigger than
the technique which is suggested here, which just adds some terms with
information on the boundaries. Moreover, time-dependent boundary conditions
are directly tackled here.

Keywords– Exponential Runge-Kutta methods nonlinear reaction-diffusion
problems avoiding order reduction in time

Mathematics Subject Classification– 65M12 65M20

##  1 Introduction

Order reduction is a well-known phenomenom that affects not only standard
integrators with stages when integrating stiff problems but also exponential-
type methods [5, 11, 16]. The latter have the advantage of integrating the
linear and stiff part of some differential problems in an exact way. This
allows to construct methods which are linearly stable and ‘explicit’ at the
same time, which is not possible with standard integrators.

Different ideas to avoid order reduction with standard integrators and
splitting methods (both exponential and not exponential) have been given in
the literature but, up to our knowledge, order greater than 2 has not been
achieved with them when integrating reaction-diffusion problems with time-
dependent boundary conditions. (Look at [3, 8, 9, 10, 19] for some of the
results with splitting methods)

In this paper, we center on trying to increase accuracy among explicit
exponential Runge-Kutta methods (EERK). Among them, Lawson ones are the
easiest to construct since, given any standard Runge-Kutta (RK) method, there
is directly a Lawson method associated to it through a change of variables
given by an exponential function [18]. However, some more methods were
constructed in the literature considering other functions defined from
exponentials and denoted by $\varphi_{j}$ ($j\geq 1$) [16]. While Lawson
methods show a severe order reduction when integrating initial boundary value
problems [1], with some other EERK methods, the order reduction is not so big.
Furthermore, in [15, 20], stiff order conditions are given when integrating
semilinear parabolic problems with vanishing boundary conditions and then,
some methods are constructed explicitly in order to achieve that order (see
also [21]). However, imposing those conditions usually implies increasing the
number of stages and thus, the computational cost.

For linear problems, a technique has already been suggested to avoid order
reduction with both Lawson and other EERK methods [2, 4] not only in parabolic
problems with vanishing boundary conditions, but also for time-dependent
boundary values. This technique is based on applying the method of lines
integrating firstly in time and then in space. Then, in the first procedure,
the exponential-type operators applied over functions with time-dependent
boundary conditions are substituted by the solution of suitable initial
boundary value problems with appropriate boundary conditions which can be
calculated exactly in terms of data of the problem. The corresponding formulas
after full discretization are also shown in [2, 4] and they are called
exponential quadrature rules because, in fact, they do not depend at all on
the stages of the method due to the fact that the problem is linear. The
exponential RK quadrature rules in [4] are in general more efficient than
Lawson ones because, in the former, all the exponential-type operators are
evaluated at the same argument and, because of that, their calculation through
Krylov subroutines in [22] is quicker [4]. Moreover, in the same paper, for
some particular methods, it is shown that the suggested technique to avoid
order reduction with quadrature rules corresponding to few nodes is more
efficient than considering the classical approach with higher stiff order
because of a higher number of nodes.

For nonlinear problems and Lawson methods, a technique has also been suggested
in [6]. The main difference with the linear case is that, in order to achieve
local order $3$ and higher, numerical differentiation is required to
approximate the necessary boundary values of the intermediate mentioned
problems. Because of that, in order to prove convergence, a CFL condition of
the form (34) is required, although this condition is much less restrictive
than that needed with standard explicit RK methods. For example, when $A$ is a
second-order space derivative operator, $\gamma=2$ would be required in (34)
with the latter. However, with Lawson methods, $\gamma=1$ is just required
with Dirichlet boundary conditions and no condition at all with Neumann ones
because $\gamma=0$ in such a case. Moreover, the additional terms which need
to be calculated in order to avoid order reduction seem to be very cheap to
calculate with Krylov techniques [6].

The aim of this paper is to see how to apply a similar technique with the EERK
methods in [15], for which $\varphi_{j}$ functions ($j\geq 1$) turn up,
instead of just $\varphi_{0}(z)=e^{z}$, as it happens with Lawson methods.
According to [12], the more recent Krylov techniques are quicker with
$\varphi_{j}$-operators when $j$ increases. Moreover, as distinct from Lawson
methods, at least for the formula concerning the calculation of the numerical
solution, the $\varphi_{j}$-operators ($j\geq 0$) are always evaluated at the
same argument. Furthermore, the same applies for the calculation of the stages
when, in (12), $r$ just takes the value $i$. As with exponential quadrature
rules, this makes Krylov subroutines much quicker [22] and, because of that,
these methods are more efficient than Lawson ones. The analysis is much more
complicated than with Lawson methods but the final conclusion is the same.
Under a very mild CFL condition, order reduction can be avoided in a quite
cheap way, since the additional calculations in fact just correspond to
boundary values. That seems much cheaper than constructing methods with a high
number of stages to increase the stiff order.

The paper is structured as follows. Firstly, a preliminaries section
establishes the abstract framework for the problem to integrate and the
required assumptions to deduce the main results in the theorems. Then, the
appropriate time semidiscretization is suggested in Section 3, so that order
reduction is avoided in the local error. In Section 4, the full discretization
is performed considering some hypotheses on the space integrator. This
includes simple finite differences, but also more complicated ones (which
require a mass matrix, as the nine-point formula) and even finite elements.
The formulas to be implemented are given in (29)-(30) and, afterwards, the
result on the local error after full discretization is given, assuming that
the required boundaries in (29)-(30) are calculated exactly. As this is not
always possible in terms of data, Section 4 offers a first simplification of
those boundaries in Tables 1 and 2, which guarantee that the optimal local
error order after full discretization is conserved when searching for local
order $2$, $3$ and $4$ in the timestepsize (respectively $p=1,2,3$ in Theorem
3). Even with those simplifications, just some approximation of those
boundaries can be given when $p=1$ with Neumann/Robin boundary conditions.
Besides, when $p\geq 2$, with both Dirichlet and Neumann/Robin boundary
conditions, numerical differentiation is required. Considering then the error
which comes from this approximation, the results on the final global error
after full discretization is given in Subsection 4.6 for $p=1,2$ and $3$.
Moreover, under the parabolic assumption (35), it is justified that the global
order on the timestepsize can be the same as the local one. Finally, in
Section 5, some numerical experiments corroborate the previous results and
show that the suggested technique to avoid order reduction with EERK methods
of less stiff order can be more efficient than the classical approach with
higher stiff order ones. As the proofs of the theorems are very technical, we
have postponed the proofs to an appendix.

##  2 Preliminaries

Let $X$ and $Y$ be Banach spaces and let $A:D(A)\subset X\to X$ and
$\partial:X\to Y$ be linear operators. Our goal is to avoid order reduction
when integrating in time through EERK methods the nonlinear abstract non
homogeneous initial boundary value problem

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

The abstract setting (4) permits to cover a wide range of nonlinear
evolutionary problems governed by partial differential equations. We use the
following hypotheses, similar to the ones in [3] when avoiding the same kind
of problems with exponential splitting methods.

  1. (A1)

The boundary operator $\partial:D(A)\subset X\to Y$ is onto and $g\in
C^{1}([0,T],Y)$.

  2. (A2)

Denoting $A_{0}:D(A_{0})=Ker(\partial)\subset X\to X$, the restriction of $A$
to Ker($\partial$), at least one of these assumptions is satisfied:

    1. (a)

Ker($\partial$) is dense in $X$ and $A_{0}$ is the infinitesimal generator of
a $C_{0}$\- semigroup $\\{e^{tA_{0}}\\}_{t\geq 0}$ in $X$ of negative type
$\omega$.

    2. (b)

$D(A)$ is dense in $X$ and $A_{0}$ generates a bounded holomorphic semigroup
$\\{e^{tA_{0}}\\}_{t\geq 0}$ in $X$ of negative type $\omega$.

  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\|\leq C\|v\|,$ (5)

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

  4. (A4)

The nonlinear source $f$ belongs to $C^{1}([0,T]\times X,X)$.

  5. (A5)

The solution $u$ of (4) satisfies $u\in C^{1}([0,T],X)$, $u(t)\in D(A)$ for
all $t\in[0,T]$ and $Au\in C([0,T],X)$.

The well-posedness of problem (4) is assured by (A1)-(A4), as explained in
[3]. Besides, as it was also justified there, (A4) is quite restrictive if
$X=L^{p}(\Omega)$ with $\Omega$ a bounded domain in $\mathbb{R}^{d}$. However,
if the supremum norm is chosen, (A4) is satisfied whenever $f$ has the form

$\displaystyle f(t,u)=\phi(u)+h(t),$ (6)

with $\phi\in C^{1}(\mathbb{C},\mathbb{C})$ and $h\in C^{1}([0,T],X)$. For
simplicity, we will assume from now on that $f$ has the form (6), apart from
the following hypotheses:

  1. (A6)

There exists a natural value $\bar{m}(A)\geq 1$ such that, whenever $w\in
D(A^{l})$ and $\phi\in C^{m+\bar{m}(A)l}(\mathbb{C},\mathbb{C})$ for natural
$l$ and $m$, $\phi^{(m)}(w)\in D(A^{l})$.

  2. (A7)

For every natural $l\geq 0$, there exists a norm $\|\cdot\|_{l}$ in
$D(A^{l})\subset X$ such that, for every natural $m\geq 1$ and $u,v\in
D(A^{l})$, whenever $\phi\in C^{m+\bar{m}(A)l}(\mathbb{C},\mathbb{C})$,
$\phi^{(m)}(u)v^{m}\in D(A^{l})$ and

$\|A^{l}\big{[}\phi^{(m)}(u)v^{m}]\big{]}\|\leq C(\phi,u)\|v\|_{l}^{m},$

for some constant $C(\phi,u)$ which depends on
$\phi^{(m)}(u),\dots,\phi^{(m+\bar{m}(A)l)}(u)$ and $\|u\|_{l}$.

  3. (A8)

For every natural $l\geq 1$ and the norm $\|\cdot\|_{l}$ in (A7), for every
$u,v\in D(A^{l})$, whenever $\phi\in
C^{(1+\bar{m}(A)l)}(\mathbb{C},\mathbb{C})$,

$\|A_{0}^{-1}A^{l}\big{[}\phi^{\prime}(u)v\big{]}\|\leq\bar{C}(\phi,u)\|v\|_{l-1},$

for some constant $\bar{C}(\phi,u)$ which depends on
$\phi^{\prime}(u),\dots,\phi^{(1+\bar{m}(A)l)}(u)$ and $\|u\|_{l}$.

Notice that, when $A$ is a space differential operator, (A6) means that $\phi$
and some of its derivatives do not reduce the regularity in space of their
arguments for regular enough $\phi$. As for (A7), it assures that Taylor
expansions with the unbounded operator $A$ can be performed when there is
enough regularity. Finally, (A8) is used to apply a summation-by-parts
argument which explains that the local order is the same as the global one.
(These assumptions will be justified for the example in the numerical
experiments.)

In the remaining of the paper, we always suppose that (A1)-(A8) are satisfied.
However, more regularity will be assumed in certain results.

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 [16] and which are defined by

$\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,\quad
j\geq 1.$ (7)

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

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

For the time integration, we will center on EERK methods which, when applied
to a finite-dimensional nonlinear problem like

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

with $M$ a certain matrix, read like this at each step

$\displaystyle K_{n,i}$ $\displaystyle=$ $\displaystyle
e^{c_{i}kM}U_{n}+k\sum_{j=1}^{i-1}a_{ij}(kM)F(t_{n}+c_{j}k,K_{n,j}),\quad
i=1,\dots,s,$ (10) $\displaystyle U_{n+1}$ $\displaystyle=$ $\displaystyle
e^{kM}U_{n}+k\sum_{i=1}^{s}b_{i}(kM)F(t_{n}+c_{i}k,K_{n,i}).$ (11)

Here, $k>0$ is the time stepsize, $t_{n}=t_{0}+nk$ and, in principle,
$a_{ij}(kM)$ and $b_{i}(kM)$ are linear combinations of
$\\{\varphi_{l}(ckM)\\}_{l=1}^{s}$, for some real value $c$. More precisely,
we will assume that

$\displaystyle a_{ij}(z)$ $\displaystyle=$
$\displaystyle\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\varphi_{l}(c_{r}z),$ (12)
$\displaystyle b_{i}(z)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{s}\mu_{i,l}\varphi_{l}(z),$

for some constants $\lambda_{i,j,l,r}$ and $\mu_{i,l}$. We notice that,
usually, in the formula for $a_{ij}$, $l$ just varies from $1$ to $i-1$ and
$r$ just takes the value $i$, but we consider the more general case because it
does not imply any difficulty in the analysis. We also notice that, usually,
the underlying method satisfies $\sum_{j=1}^{s}a_{ij}(0)=c_{i}$, which implies
that $\sum_{j,l,r=1}^{s}\lambda_{i,j,l,r}/l!=c_{i}$. We will assume this from
now on for the simplicity of some expressions in Section 4.4.

##  3 Suggestion for the time semidiscretization

As already mentioned in the introduction, we will integrate (4) firstly in
time in order to avoid order reduction. However, (4) is not an ODE system like
(9), but an initial boundary value problem. Because of that, for each of the
terms in (10)-(11), we must suggest the solution of a suitable initial
boundary value problem. Simulating what was already done with exponential
functions and Lawson methods [2], whenever $\alpha\in X$ and $\beta\in
D(A^{r+1})$, we consider the generalised solution of this problem, which we
will denote by $\phi_{0,r,\alpha,\beta}(\tau)$,

$\displaystyle v^{\prime}(\tau)$ $\displaystyle=$ $\displaystyle Av(\tau),$
$\displaystyle v(0)$ $\displaystyle=$ $\displaystyle\alpha,$
$\displaystyle\partial v(\tau)$ $\displaystyle=$
$\displaystyle\partial(\sum_{l=0}^{r}\frac{\tau^{l}}{l!}A^{l}\beta).$ (13)

On the other hand, in a similar way to what was done for
$\varphi_{j}$-functions and exponential quadrature rules [4], and under the
same assumptions on $\alpha$ and $\beta$, we also consider the solution of
this other problem, which we will denote by $\phi_{j,r,\alpha,\beta}(\tau)$,

$\displaystyle v^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A-\frac{j}{\tau})v(\tau)+\frac{1}{(j-1)!\tau}\alpha,$
$\displaystyle v(0)$ $\displaystyle=$ $\displaystyle\frac{1}{j!}\alpha,$
$\displaystyle\partial v(\tau)$ $\displaystyle=$
$\displaystyle\partial(\sum_{l=0}^{r}\frac{\tau^{l}}{(l+j)!}A^{l}\beta).$ (14)

As it was already proved in Lemma 3.1 of [2] for $\phi_{0,r,\alpha,\beta}$,
and as it can be proved in a similar way for $\phi_{j,r,\alpha,\beta}$,
considering Lemmas 6 and 7 in [4], for $r\geq 0$,

$\displaystyle\phi_{0,r,\alpha,\beta}(\tau)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{r}\frac{\tau^{l}}{l!}A^{l}\beta+e^{\tau
A_{0}}(\alpha-\beta)+\tau^{r+1}\varphi_{r+1}(\tau A_{0})A^{r+1}\beta,$ (15)
$\displaystyle\phi_{j,r,\alpha,\beta}(\tau)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{r}\frac{\tau^{l}}{(l+j)!}A^{l}\beta+\varphi_{j}(\tau
A_{0})(\alpha-\beta)+\tau^{r+1}\varphi_{j+r+1}(\tau A_{0})A^{r+1}\beta.$ (16)

Then, when the EERK method has non-stiff order $p$, our suggestion in order to
avoid order reduction, when the problem is regular enough, is to consider as
stages, from the previous numerical value $u_{n}$,

$\displaystyle
K_{n,i}=\phi_{0,p-1,u_{n},u(t_{n})}(c_{i}k)+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\phi_{l,p-2,F_{n,j},\bar{F}_{n,j}}(c_{r}k),$
(17)

where

$\displaystyle
F_{n,j}=f(t_{n}+c_{j}k,K_{n,j}),\quad\bar{F}_{n,j}=f(t_{n}+c_{j}k,\bar{K}_{n,j}),$
(18)

with $\bar{K}_{n,j}$ obtained as $K_{n,j}$ but substituting $u_{n}$ by
$u(t_{n})$ and, recursively, $F_{n,l}$ by $\bar{F}_{n,l}$. Then, we suggest as
the numerical solution after timestepsize $k$,

$\displaystyle
u_{n+1}=\phi_{0,p,u_{n},u(t_{n})}(k)+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\phi_{l,p-1,F_{n,i},\bar{F}_{n,i}}(k).$
(19)

We notice that, when $p=1$, instead of (17), it should be used

$\displaystyle
K_{n,i}=\phi_{0,0,u_{n},u(t_{n})}(c_{i}k)+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\varphi_{l}(c_{r}kA_{0})F_{n,j},$
(20)

since $\phi_{l,p-2,F_{n,j},\bar{F}_{n,j}}$ should be the solution of (14),
with vanishing boundary conditions.

###  3.1 Local error

Before stating and proving the result on the local error which assures that
order reduction would be avoided in such a way under enough assumptions of
regularity, we remind that the conditions for non-stiff order $p$ come from
considering the Taylor series of $U(t_{n+1})$ in terms of $U(t_{n})$ when
$U(t)$ is the solution of the ordinary differential system (9), i.e.,

$\displaystyle U(t_{n+1})$ $\displaystyle=$ $\displaystyle
U(t_{n})+k[MU(t_{n})+F_{n}]$ (21)
$\displaystyle+\frac{k^{2}}{2}[M^{2}U(t_{n})+MF_{n}+F_{t,n}+F_{u,n}MU(t_{n})+F_{u,n}F_{n}]$
$\displaystyle+\sum_{j=3}^{p}\frac{k^{j}}{j!}\sum_{m=1}^{n(j)}\gamma_{m,j}ED_{m,j,n}(F,M)(t_{n},U(t_{n}))+O(k^{p+1}),$

where $F_{*,n}$ denotes the corresponding evaluation of $F$ or some of its
derivatives when evaluated at $(t_{n},U(t_{n}))$ and
$\\{ED_{m,j,n}(F,M)\\}_{m=1}^{n(j)}$ are the independent elementary
differentials associated to the power $k^{j}$ ($0\leq j\leq p$) and also
evaluated at $(t_{n},U(t_{n}))$ [13]. The order conditions come from equating
those coefficients $\gamma_{m,j}$ with the ones corresponding to the Taylor
series of the corresponding numerical solution (10)-(11) starting from
$U(t_{n})$. More precisely, the result of developing in Taylor series around
$(t_{n},U(t_{n}))$,

$\displaystyle\bar{U}_{n+1}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{p}\frac{k^{l}}{l!}M^{l}U(t_{n})$ (22)
$\displaystyle+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\sum_{j=0}^{p-1}\frac{k^{j}}{(j+l)!}M^{j}F(t_{n}+c_{i}k,\bar{K}_{n,i})+O(k^{p+1}),$

where, recursively,

$\displaystyle\bar{K}_{n,i}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{p-1}\frac{(c_{i}k)^{l}}{l!}M^{l}U(t_{n})$ (23)
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\sum_{ll=0}^{p-2}\frac{(c_{r}k)^{ll}}{(ll+l)!}M^{ll}F(t_{n}+c_{j}k,\bar{K}_{n,j})+O(k^{p}).$

For the numerical integration of (4) through (17)-(19), we then have the
following result for the local error $\rho_{n}=\bar{u}_{n+1}-u(t_{n+1})$,
where $\bar{u}_{n+1}$ is defined also through (17)-(19), but starting from
$u(t_{n})$ instead of $u_{n}$.

######  Theorem 1.

Under hypotheses (A1)-(A8), if the explicit exponential Runge-Kutta method
(10)-(11) has non-stiff order $\geq p$, when integrating (4) through (17)-(19)
with $u\in C([0,T],D(A^{p+1}))\cap C^{p+1}([0,T],X)$, $\phi\in
C^{1+(p-1)\bar{m}(A)}(\mathbb{C},\mathbb{C})$, $h\in C^{p}([0,T],D(A^{p})),$
it happens that the local error satisfies $\rho_{n}=O(k^{p+1})$.

Moreover, if in formulas (17)-(19), we had substituted $p$ by $\hat{p}$ with
$1\leq\hat{p}\leq p-1$, whenever $u\in C([0,T],D(A^{\hat{p}+2}))\cap
C^{\hat{p}+2}([0,T],X)$, $\phi\in
C^{1+\hat{p}\bar{m}(A)}(\mathbb{C},\mathbb{C})$ and $h\in
C^{\hat{p}+1}([0,T],D(A^{\hat{p}+1}))$, it happens that not only the local
error satisfies $\rho_{n}=O(k^{\hat{p}+1})$ but also
$A_{0}^{-1}\rho_{n}=O(k^{\hat{p}+2})$.

The result of the previous theorem corresponding to $\hat{p}<p$ is interesting
because, for the non-stiff order $p$, some times the corresponding boundaries
in (13) and (14) cannot be calculated exactly in terms of the data of the
problem. We can always resort to numerical differentiation to approximate
them, but the more terms which need to be calculated in such a way, the more
complicated and unstable the algorithm may become. That will be studied in
detail afterwards.

##  4 Suggestion for the full discretization

As distinct from the standard approach when integrating partial differential
problems, we suggest to discretize firstly in time as described in the
previous section and then to integrate in space the resulting auxiliary
problems. For the suggestion of a space discretization, we must first fix a
certain space $X$. As in [3] and [6], we take $X=C(\overline{\Omega})$ for a
certain bounded domain $\Omega\in\mathbb{R}^{d}$ and the maximum norm.

###  4.1 Space discretization

In this subsection we consider a certain grid $\Omega_{h}$ (of $\Omega$) over
which the approximated numerical solution will be defined. Therefore, the
numerical approximation belongs to $\mathbb{C}^{N}$, where $N$ is the number
of nodes in the grid, and the maximum norm $\|\cdot\|_{h}$ in this space will
also be considered. Notice that, with Dirichlet boundary conditions, nodes on
the boundary are not part of the approximation while, when using Neumann or
Robin boundary conditions, the approximation at those nodes must be
considered.

We will make similar hypotheses for the space discretization as those in [3]
and [4]. We state them here for the sake of clarity. We denote the projection
operator

$\displaystyle P_{h}:X\to\mathbb{C}^{N},$

to that which takes a function to its values over the grid $\Omega_{h}$.
Notice that the following bound then holds

$\displaystyle\|P_{h}u\|_{h}\leq\|u\|_{h}.$ (24)

Then, we assume that the discretization of

$Au=F,\qquad\partial u=g,$

is given by

$\displaystyle A_{h,0}U_{h}+C_{h}g=P_{h}F+D_{h}\partial F,$ (25)

where $A_{h,0}$ is the matrix which discretizes $A_{0}$ and
$C_{h},D_{h}:Y\to\mathbb{C}^{N}$ are other operators associated to the
discretization of $A$, which take into account the information on the boundary
of $u$ and $F$.

Notice that, as $f$ is given by (6), it has sense as a function from
$[0,T]\times\mathbb{C}^{N}$ on $\mathbb{C}^{N}$ and, for each $t\in[0,T]$ and
$u\in X$,

$\displaystyle P_{h}f(t,u)=f(t,P_{h}u).$ (26)

Moreover, we consider the following hypotheses:

  1. (H1)

The matrix $A_{h,0}$ satisfies

    1. (a)

$A_{h,0}$ is invertible and $\|A_{h,0}^{-1}\|_{h}\leq C$ for some
$h$-independent constant $C$.

    2. (b)

$\|\varphi_{j}(tA_{h,0})\|_{h}\leq C_{j},\quad j=0,\dots,s,\quad t\in[0,T]$.

  2. (H2)

We define the elliptic projection $R_{h}:D(A)\to\mathbb{C}^{N}$ as the
solution of

$\displaystyle A_{h,0}R_{h}u+C_{h}\partial u=P_{h}Au+D_{h}\partial Au.$ (27)

    1. (a)

There exists a subspace $Z\subset D(A)$ such that, for $u\in Z$,
$A_{0}^{-1}u\in Z$ and $e^{tA_{0}}u\in Z$, for $t\in[0,T]$. Besides, there
exists a natural value $\bar{m}(Z)\geq 1$ such that, whenever $w\in D(A^{l})$
for natural $l$ and $A^{l}w\in Z$, it happens that $A^{l}\phi(w)\in Z$ if
$\phi\in C^{\bar{m}(Z)+l\bar{m}(A)}(\mathbb{C},\mathbb{C})$. Moreover, for
some $\varepsilon_{h}$ and $\eta_{h}$ which are both small with $h$, if $u\in
Z$,

$\left\|A_{h,0}({P_{h}u-R_{h}u})\right\|_{h}\leq\varepsilon_{h}\left\|u\right\|_{Z},\quad\left\|P_{h}u-R_{h}u\right\|_{h}\leq\eta_{h}\left\|u\right\|_{Z}.$
(28)

(Although obviously, because of (H1a), $\eta_{h}$ could be taken as
$C\varepsilon_{h}$, for some discretizations $\eta_{h}$ can decrease more
quickly with $h$ than $\varepsilon_{h}$ and that leads to better error bounds
in the following sections.)

    2. (b)

$\|A_{h,0}^{-1}C_{h}\|_{h}\leq C^{\prime\prime}$ for some constant
$C^{\prime\prime}$ which does not depend on $h$. This resembles the continuous
maximum principle which is satisfied because of (5) when $z=0$.

    3. (c)

$\|D_{h}\|_{h}$ is uniformly bounded on $h$.

  3. (H3)

The nonlinear source $f$ belongs to
$C^{1}([0,T]\times\mathbb{C}^{N},\mathbb{C}^{N})$ and the derivative with
respect to the variable in $\mathbb{C}^{N}$ is uniformly bounded in a
neighbourhood of the solution where the numerical approximation stays.

###  4.2 Full discretization formulas depending on $\partial A^{l}u(t_{n})$
and $\partial A^{l}\bar{F}_{n,i}$

As in (17) and (19), $\alpha=u_{n}$ in the part which corresponds to (13) and
$\alpha=F_{n,i}$ in that corresponding to (14), we suggest to discretize (13)
and (14) in space considering (25) and taking as initial condition for the
space discretized system $U_{h}^{n}$ for (13) with boundary $\partial
u(t_{n})$ and $F_{n,i,h}$ for (14) with boundary $\partial\bar{F}_{n,i}$.

More precisely, for $r=p-1$ or $r=p$, depending on whether we are calculating
the stages or the solution, we consider the space discretized system related
to (13):

$\displaystyle V_{h}^{\prime}(\tau)$ $\displaystyle=$ $\displaystyle
A_{h,0}V_{h}(\tau)+C_{h}\partial[\sum_{l=0}^{r}\frac{\tau^{l}}{l!}A^{l}u(t_{n})]-D_{h}\partial[\sum_{l=1}^{r}\frac{\tau^{l-1}}{(l-1)!}A^{l}u(t_{n})],$
$\displaystyle V_{h}(0)$ $\displaystyle=$ $\displaystyle U_{h}^{n}.$

By the variation-of-constants formula and the definition of $\varphi_{l}$ (7),
it is clear that

$\displaystyle\phi_{0,r,U_{h}^{n},u(t_{n})}^{h}(\tau)$ $\displaystyle\equiv$
$\displaystyle e^{\tau
A_{h,0}}U_{h}^{n}+\sum_{l=0}^{r-1}\tau^{l+1}\varphi_{l+1}(\tau
A_{h,0})[C_{h}\partial A^{l}u(t_{n})-D_{h}\partial A^{l+1}u(t_{n})]$
$\displaystyle+\tau^{r+1}\varphi_{r+1}(\tau A_{h,0})C_{h}\partial
A^{r}u(t_{n})$

solves this problem. (We consider this notation for similarity with the
continuous one in (15).) In a similar way, for the space discretized system
related to (14), taking $r=p-2$ or $r=p-1$ depending on whether we are
calculating the stages or the solution, we consider

$\displaystyle V_{h}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A_{h,0}-\frac{j}{\tau})V_{h}(\tau)+\frac{1}{(j-1)!\tau}F_{n,i,h}+C_{h}\partial[\sum_{l=0}^{r}\frac{\tau^{l}}{(l+j)!}A^{l}\bar{F}_{n,i}]$
$\displaystyle-
D_{h}\partial[\sum_{l=1}^{r}\frac{l\tau^{l-1}}{(l+j)!}A^{l}\bar{F}_{n,i}]-jD_{h}\partial[\sum_{l=0}^{r}\frac{\tau^{l-1}}{(l+j)!}A^{l}\bar{F}_{n,i}]+\frac{1}{(j-1)!\tau}D_{h}\partial\bar{F}_{n,i},$

which can be written like this for $r>-1$,

$\displaystyle V_{h}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A_{h,0}-\frac{j}{\tau})V_{h}(\tau)+\frac{1}{(j-1)!\tau}F_{n,i,h}$
$\displaystyle+\sum_{l=0}^{r-1}\frac{\tau^{l}}{(l+j)!}[C_{h}\partial
A^{l}\bar{F}_{n,i}-D_{h}\partial
A^{l+1}\bar{F}_{n,i}]+\frac{\tau^{r}}{(r+j)!}C_{h}\partial A^{r}\bar{F}_{n,i}$
$\displaystyle V_{h}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{j!}F_{n,i,h}.$

Then, the solution of this, in the same way that it was argued for exponential
quadrature rules in [4], can be written like this for $\tau>0$ and $r>-1$,

$\displaystyle\phi_{j,r,F_{n,i,h},\bar{F}_{n,i}}^{h}(\tau)$
$\displaystyle\equiv$ $\displaystyle\varphi_{j}(\tau
A_{h,0})F_{n,i,h}+\sum_{l=0}^{r-1}\tau^{l+1}\varphi_{j+l+1}(\tau
A_{h,0})[C_{h}\partial A^{l}\bar{F}_{n,i}-D_{h}\partial A^{l+1}\bar{F}_{n,i}]$
$\displaystyle+\tau^{r+1}\varphi_{j+r+1}(\tau A_{h,0})C_{h}\partial
A^{r}\bar{F}_{n,i}.$

Considering all this, the suggestion for the full discretization formulas is,
in principle:

$\displaystyle
K_{n,i,h}=e^{c_{i}kA_{h,0}}U_{h}^{n}+\sum_{l=0}^{p-2}(c_{i}k)^{l+1}\varphi_{l+1}(c_{i}kA_{h,0})[C_{h}\partial
A^{l}u(t_{n})-D_{h}\partial A^{l+1}u(t_{n})]$ (29) $\displaystyle\hskip
85.35826pt+(c_{i}k)^{p}\varphi_{p}(c_{i}kA_{h,0})C_{h}\partial
A^{p-1}u(t_{n})$
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}[\varphi_{l}(c_{r}kA_{h,0})F_{n,j,h}+\sum_{ll=0}^{p-3}(c_{r}k)^{ll+1}\varphi_{l+ll+1}(c_{r}kA_{h,0})[C_{h}\partial
A^{ll}\bar{F}_{n,j}-D_{h}\partial A^{ll+1}\bar{F}_{n,j}]$ $\displaystyle\hskip
85.35826pt+(c_{r}k)^{p-1}\varphi_{l+p-1}(c_{r}kA_{h,0})C_{h}\partial
A^{p-2}\bar{F}_{n,j}],$ $\displaystyle
U_{h}^{n+1}=e^{kA_{h,0}}U_{h}^{n}+\sum_{l=0}^{p-1}k^{l+1}\varphi_{l+1}(kA_{h,0})[C_{h}\partial
A^{l}u(t_{n})-D_{h}\partial A^{l+1}u(t_{n})]$ (30) $\displaystyle\hskip
85.35826pt+k^{p+1}\varphi_{p+1}(kA_{h,0})C_{h}\partial A^{p}u(t_{n})$
$\displaystyle+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\bigg{[}\varphi_{l}(kA_{h,0})F_{n,i,h}+\sum_{ll=0}^{p-2}k^{ll+1}\varphi_{l+ll+1}(kA_{h,0})[C_{h}\partial
A^{ll}\bar{F}_{n,i}-D_{h}\partial A^{ll+1}\bar{F}_{n,i}]$ $\displaystyle\hskip
85.35826pt+k^{p}\varphi_{l+p}(kA_{h,0})C_{h}\partial
A^{p-1}\bar{F}_{n,i}\bigg{]},$

where $F_{n,i,h}=f(t_{n}+c_{i}k,K_{n,i,h})$ and where we notice that, for
$p=1$, due to (20), (29) should read

$\displaystyle K_{n,i,h}$ $\displaystyle=$ $\displaystyle
e^{c_{i}kA_{h,0}}U_{h}^{n}+c_{i}k\varphi_{1}(c_{i}kA_{h,0})C_{h}\partial
u(t_{n})$ (31)
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\varphi_{l}(c_{r}kA_{h,0})F_{n,j,h}.$

###  4.3 Local error of the full discretization

We define the local error of the full discretization by
$\rho_{n,h}=\bar{U}_{h}^{n+1}-P_{h}u(t_{n+1})$, where $\bar{U}_{h}^{n+1}$ is
defined through (29)-(30) but substituting $U_{h}^{n}$ by $P_{h}u(t_{n})$. The
following theorem states how this error behaves.

######  Theorem 2.

Under the first set of hypotheses of Theorem 1, (H1)-(H3) and assuming also
that $\phi\in C^{\bar{m}(Z)+p\bar{m}(A)}(\mathbb{C},\mathbb{C})$,

$\displaystyle A^{l}u\in C([0,T],Z),\quad l=0,\dots,p+1,\quad A^{l}h\in
C([0,T],Z),\quad l=0,\dots,p,$ (32)

it happens that $\rho_{n,h}=O(k^{p+1}+k\varepsilon_{h})$.

Moreover, under the second set of hypotheses of Theorem 1,(H1)-(H3) and
assuming also that $\phi\in
C^{\bar{m}(Z)+\hat{p}\bar{m}(A)}(\mathbb{C},\mathbb{C})$ and (32), but with
$\hat{p}$ instead of $p$, it happens that, not only
$\rho_{n,h}=O(k^{\hat{p}+1}+k\varepsilon_{h})$ but also
$A_{h,0}^{-1}\rho_{n,h}=O(k^{\hat{p}+2}+k\eta_{h}+k^{2}\varepsilon_{h})$.

###  4.4 Simplification of the required boundaries without losing order

As it also happened with Lawson methods and nonlinear problems [6], the exact
calculation of the required boundaries is not usually possible in terms of the
data of the problem and just some numerical approximation of those values can
be considered. Because of that, before considering any numerical approximation
of them, we firstly suggest some simplifications of the previously suggested
boundary values which do not alter the local error after full discretization.
More precisely, for $p=1,2,3$ in (29)-(30), we suggest some simplifications
which do not alter the respective order $p+1$ for the local error after full
discretization, which has just been proved in the previous subsection. More
precisely, we suggest the simplifications for $\partial A^{l}\bar{F}_{n,i}$ in
Tables 1 and 2, which come from developing into Taylor series $\bar{K}_{n,i}$
in (23) around $t_{n}$. (Notice that the simplification is different depending
on whether the term is in the stages (29) or in the approximation (30)).

In the following, we will state more precisely under which assumptions the
order of the full discretization local error after that simplification (which
we will denote by $\rho_{n,h}^{simp}$) will be conserved.

Table 1: Simplifications for the boundaries in (29) $p$ | Simplification in $K_{n,i,h}$
---|---
2 | $\partial\bar{F}_{n,i}\approx\partial f(t_{n},u(t_{n}))$
3 | $\begin{array}[]{rcl}\partial\bar{F}_{n,i}&\approx&\partial f(t_{n}+c_{i}k,u(t_{n})+c_{i}k\dot{u}(t_{n}))\\\ \partial A\bar{F}_{n,i}&\approx&\partial Af(t_{n},u(t_{n}))\end{array}$
Table 2: Simplifications for the boundaries in (30) $p$ | Simplification in $U_{h}^{n+1}$
---|---
1 | $\partial\bar{F}_{n,i}\approx\partial f(t_{n},u(t_{n}))$
2 | $\begin{array}[]{rcl}\partial\bar{F}_{n,i}&\approx&\partial f(t_{n}+c_{i}k,u(t_{n})+c_{i}k\dot{u}(t_{n}))\\\ \partial A\bar{F}_{n,i}&\approx&\partial Af(t_{n},u(t_{n}))\end{array}$
3 | $\begin{array}[]{rcl}\partial\bar{F}_{n,i}&\approx&\partial f\bigg{(}t_{n}+c_{i}k,u(t_{n})+kc_{i}Au(t_{n})+\frac{(c_{i}k)^{2}}{2}A^{2}u(t_{n})\\\ &&\hskip 28.45274pt+k\sum_{j,l,r}\lambda_{i,j,l,r}\big{[}\frac{1}{l!}f\big{(}t_{n}+c_{j}k,u(t_{n})+c_{j}k\dot{u}(t_{n}))+\frac{c_{r}k}{(l+1)!}Af(t_{n},u(t_{n}))\big{]}\bigg{)}\\\ \partial A\bar{F}_{n,i}&\approx&\partial Af(t_{n}+c_{i}k,u(t_{n})+c_{i}k\dot{u}(t_{n}))\\\ \partial A^{2}\bar{F}_{n,i}&\approx&\partial A^{2}f(t_{n},u(t_{n}))\end{array}$

######  Theorem 3.

Whenever $p=1,2,3$, under hypotheses (A1)-(A8) and (H1)-(H3), if the EERK
method (10)-(11) has non-stiff order $\geq p$, when integrating (4) through
(29)-(30) with $u\in C([0,T],D(A^{p+1}))\cap C^{p+1}([0,T],X)$, $\phi\in
C^{\bar{m}(Z)+p\bar{m}(A)}(\mathbb{C},\mathbb{C})$, $h\in
C^{p}([0,T],D(A^{p}))$ and (32), it happens that the local error of the full
discretization with the simplified boundaries for $U_{h}^{n+1}$ in Tables 1
and 2 satisfies $\rho_{n,h}^{simp}=O(k^{p+1}+k\varepsilon_{h})$.

Furthermore, if the method has non-stiff order $\geq p+1$, $u\in
C([0,T],D(A^{p+2}))\cap C^{p+2}([0,T],X)$ and $h\in C^{p+1}([0,T],D(A^{p+1}))$
and, just for the case $p\geq 2$, the following bound holds

$\displaystyle\|A_{0}^{-1}[\phi^{\prime}(u(t))A_{0}w]\|\leq C\|w\|,\mbox{ for
every }t\in[0,T],\,w\in D(A_{0}),$ (33)

it happens that
$A_{h,0}^{-1}\rho_{n,h}^{simp}=O(k^{p+2}+k\eta_{h}+k^{2}\varepsilon_{h})$.

###  4.5 Calculation of the required boundaries

In this subsection, we will state the error which is committed when
approximating not only the corresponding simplified boundaries in (29)-(30)
described in Tables 1 and 2, but also the error in the approximation of
$\partial A^{l}u(t_{n})$ for the different values of $l$. The latter and part
of the former was already studied in detail in [6] for Dirichlet and
Robin/Neumann boundary conditions when considering Lawson methods. There are
some differences in the arguments where the terms of the form $\partial
A^{l}f$ must be evaluated now, but the error which is committed when
approximating these boundaries in terms of the data of the problem (4),
through the numerical approximation of the problem itself, or using numerical
differentiation, is very similar. We state it here without entering into the
details but reminding that $e_{n,h}=U_{h}^{n}-P_{h}u(t_{n})$ (and therefore a
bound for this error turns up when substituting $u(t_{n})$ at some node in the
boundary by the numerical approximation itself), $\mu_{k,1}$ is the norm of
the error of some numerical differentiation formula for approximating the
first derivative in time, $\nu_{h}$ the norm of the error of some numerical
differentiation formula for approximating a certain $\gamma$th-derivative in
space and $\mu_{k,2}$ the norm of the error of some numerical differentiation
formula for approximating the second derivative in time.

For simplicity, we will separate the cases $p=1,2,3$ (which correspond
respectively to Tables 3,4 and 5), distinguishing between the error which is
committed for every term in the boundaries for the stages and for the
numerical solution, and taking into account the term to be approximated in
formulas (29)-(30) through their simplifications in Tables 1 and 2. We also
notice that, in 1-dimensional problems, when $A=D_{xx}$, as justified in [6]
for $p=2,3$, $\gamma=1$ with Dirichlet boundary conditions and $\gamma=0$ for
Robin/Neumann ones.

Table 3: Errors which are committed at each step when approximating the corresponding simplifications of the boundary terms in Tables 1 and 2 to achieve local order $2$ $p=1$ |  | Dirichlet | Robin/Neumann
---|---|---|---
$K_{n,i,h}$ | $\partial u(t_{n})$ | - | -
$U_{h}^{n+1}$ | $\partial u(t_{n})$ | - | -
$\partial Au(t_{n})/\partial\bar{F}_{n,i}$ | - | $\displaystyle{O(\|e_{n,h}\|_{h})}$
Table 4: Errors which are committed at each step when approximating the corresponding simplifications of the boundary terms in Tables 1 and 2 to achieve local order $3$ $p=2$ |  | Dirichlet | Robin/Neumann
---|---|---|---
$K_{n,i,h}$ | $\partial u(t_{n})$ | - | -
$\partial Au(t_{n})/\partial\bar{F}_{n,i}$ | - | $O(\|e_{n,h}\|_{h})$
$U_{h}^{n+1}$ | $\partial u(t_{n})$ | - | -
$\partial Au(t_{n})$ | - | $O(\|e_{n,h}\|_{h})$
$\partial\bar{F}_{n,i}$ | - | $O(\|e_{n,h}\|_{h}+k\mu_{k,1})$
$\partial A^{2}u(t_{n})/\partial A\bar{F}_{n,i}$ | $O(\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})$ | $\displaystyle{O(\mu_{k,1}+\frac{\|e_{n,h}\|_{h}}{k}+\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})}$
Table 5: Errors which are committed at each step when approximating the corresponding simplifications of the boundary terms in Tables 1 and 2 to achieve local order $4$ $p=3$ |  | Dirichlet | Robin/Neumann
---|---|---|---
$K_{n,i,h}$ | $\partial u(t_{n})$ | - | -
$\partial Au(t_{n})$ | - | $O(\|e_{n,h}\|_{h})$
$\partial\bar{F}_{n,i}$ | - | $O(\|e_{n,h}\|_{h}+k\mu_{k,1})$
$\partial A^{2}u(t_{n})/\partial A\bar{F}_{n,i}$ | $O(\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})$ | $O(\mu_{k,1}+\frac{\|e_{n,h}\|_{h}}{k}+\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})$
$U_{h}^{n+1}$ | $\partial u(t_{n})$ | - | -
$\partial Au(t_{n})$ | - | $O(\|e_{n,h}\|_{h})$
$\partial\bar{F}_{n,i}$ | $O(k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}\|e_{n,h}\|_{h})$ |  $O(\|e_{n,h}\|_{h}+k\mu_{k,1}+k^{2}\mu_{k,2}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}\|e_{n,h}\|_{h}$)
$\partial A^{2}u(t_{n})$ | $O(\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})$ | $O(\mu_{k,1}+\frac{\|e_{n,h}\|_{h}}{k}+\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})$
$\partial A\bar{F}_{n,i}$ | $O(\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}}+\frac{k\mu_{k,1}}{h^{\gamma}})$ | $O(\mu_{k,1}+k\mu_{k,2}+\frac{\|e_{n,h}\|_{h}}{k}+\nu_{h}+\frac{\|e_{n,h}\|_{h}}{h^{\gamma}})$
$\partial A^{3}u(t_{n})/\partial A^{2}\bar{F}_{n,i}$ | $O(\nu_{h}+\frac{\|e_{n,h}\|_{h}}{kh^{\gamma}}+\frac{\mu_{k,1}}{h^{\gamma}})$ | $O(\mu_{k,1}+\mu_{k,2}+\frac{\|e_{n,h}\|_{h}}{k^{2}}+\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}})$

###  4.6 Global error of the full discretization

In this subsection we show how the global error $e_{n,h}$ behaves with our
suggestion for the implementation of EERK methods, which is based in formulas
(29)-(30), after approximating the boundaries of those formulas as suggested
in the previous subsections. For that, the classical argument which relates
local error to global error can be performed, but the error coming from the
approximation of the simplified boundaries must also be considered. Besides,
under certain assumptions related to more regularity and parabolicity, a
summation-by-parts argument can be applied and therefore, the order can
increase.

######  Theorem 4.

Let us assume the first set of hypotheses of Theorem 3 and also, just for
$p=2,3,$ that, for a certain constant $C$,

$\displaystyle\frac{k}{h^{\gamma}}\leq C,$ (34)

where $\gamma$ is the order of the space derivative which must be approximated
through numerical differentiation to calculate the necessary boundaries of the
suggested method. Then, it happens that, under Dirichlet boundary conditions,

  * $\bullet$

For $p=1$, $\|e_{n,h}\|_{h}=O(k+\varepsilon_{h})$,

  * $\bullet$

For $p=2$, $\|e_{n,h}\|_{h}=O(k^{2}+k\nu_{h}+\varepsilon_{h})$,

  * $\bullet$

For $p=3$, $\|e_{n,h}\|_{h}=O(k^{3}+k\nu_{h}+k\mu_{k,1}+\varepsilon_{h})$,

and, under R/N boundary conditions,

  * $\bullet$

For $p=1$, $\|e_{n,h}\|_{h}=O(k+\varepsilon_{h})$,

  * $\bullet$

For $p=2$, $\|e_{n,h}\|_{h}=O(k^{2}+k\mu_{k,1}+k\nu_{h}+\varepsilon_{h})$,

  * $\bullet$

For $p=3$,
$\|e_{n,h}\|_{h}=O(k^{3}+k\mu_{k,1}+k^{2}\mu_{k,2}+k\nu_{h}+\varepsilon_{h})$.

Assuming also the second set of hypotheses of the same theorem, that the
following condition holds for a constant $C$ which is independent of $k$ and
$h$

$\displaystyle\|kA_{h,0}\sum_{r=1}^{n-1}e^{rkA_{h,0}}\|_{h}\leq C,\quad 0\leq
nk\leq T,$ (35)

and that, for $t\in[0,T]$, $\dot{u}(t)\in D(A^{p+1})$ and

$\displaystyle A^{l}\dot{u}(t)\in Z,\quad l=0,1,\dots,p+1,$ (36)

it happens that, under Dirichlet boundary conditions,

  * $\bullet$

For $p=1$, $\|e_{n,h}\|_{h}=O(k^{2}+k\varepsilon_{h}+\eta_{h})$,

  * $\bullet$

For $p=2$, $\|e_{n,h}\|_{h}=O(k^{3}+k\nu_{h}+k\varepsilon_{h}+\eta_{h})$,

  * $\bullet$

For $p=3$,
$\|e_{n,h}\|_{h}=O(k^{4}+k\nu_{h}+k\mu_{k,1}+k\varepsilon_{h}+\eta_{h})$,

and, under R/N boundary conditions,

  * $\bullet$

For $p=1$, $\|e_{n,h}\|_{h}=O(k^{2}+k\varepsilon_{h}+\eta_{h})$,

  * $\bullet$

For $p=2$,
$\|e_{n,h}\|_{h}=O(k^{3}+k\mu_{k,1}+k\nu_{h}+k\varepsilon_{h}+\eta_{h})$,

  * $\bullet$

For $p=3$,
$\|e_{n,h}\|_{h}=O(k^{4}+k\mu_{k,1}+k^{2}\mu_{k,2}+k\nu_{h}+k\varepsilon_{h}+\eta_{h})$.

Notice that condition (35) was also used in [3] and was proved in [14] for
analytic semigroups, covering the case in which the linear operator in (4)
corresponds to a parabolic problem.

On the other hand, we remark that, for $p=2,3$, a CFL condition (34) which, in
any case, is much less restrictive than that needed with explicit Runge-Kutta
methods is required.

##  5 Numerical experiments

In this section, we have made some numerical experiments in order to
corroborate the results of the previous sections.

###  5.1 One-dimensional problem

Table 6: Local and global error when integrating problem (37) with nonvanishing boundary conditions (38) with method (43) without avoding order reduction. k | 1/20 | 1/40 | 1/80 | 1/160
---|---|---|---|---
Local error | 2.1370e-2 | 1.0531e-2 | 5.2115e-3 | 2.5816e-3
Order |  | 1.02 | 1.01 | 1.01
Global error | 2.2693e-2 | 1.1280e-2 | 5.6078e-3 | 2.7848e-3
Order |  | 1.01 | 1.01 | 1.01
Table 7: Local and global error when integrating problem (37) with nonvanishing boundary conditions (39) with method (43) without avoding order reduction. k | 1/20 | 1/40 | 1/80 | 1/160
---|---|---|---|---
Local error | 1.1214e-3 | 3.9698e-4 | 1.4089e-4 | 4.9963-5
Order |  | 1.50 | 1.49 | 1.50
Global error | 2.0719e-2 | 1.0369e-2 | 5.1725e-3 | 2.5729e-3
Order |  | 0.99 | 1.00 | 1.01
Table 8: Local and global error when integrating problem (37) with nonvanishing boundary conditions (38) with method (43) using $p=1$ in (29)-(30) and the simplification for $\partial\bar{F}_{n,i}$ in $U_{h}^{n+1}$ in Table 2 for that value of $p$. k | 1/20 | 1/40 | 1/80 | 1/160
---|---|---|---|---
Local error | 1.2438e-3 | 3.1031e-4 | 7.7348e-5 | 1.9256e-5
Order |  | 2.00 | 2.00 | 2.01
Global error | 7.5946e-4 | 1.9141e-4 | 4.8515e-5 | 1.2292e-5
Order |  | 1.99 | 1.98 | 1.98
Table 9: Local and global error when integrating problem (37) with nonvanishing boundary conditions (39) with method (43) using $p=1$ in (29)-(30) and the simplification for $\partial\bar{F}_{n,i}$ in $U_{h}^{n+1}$ in Table 2 for that value of $p$. k | 1/20 | 1/40 | 1/80 | 1/160
---|---|---|---|---
Local error | 1.2438e-3 | 3.1031e-4 | 7.7348e-5 | 1.9256e-5
Order |  | 2.00 | 2.00 | 2.01
Global error | 8.1236e-4 | 2.0456e-4 | 5.1368e-5 | 1.2857e-5
Order |  | 1.99 | 1.99 | 2.00
Table 10: Local and global error when integrating problem (37) with nonvanishing boundary conditions (38) with method (43) using $p=2$ in (29)-(30) and the simplification for the boundaries in Tables 1 and 2 for that value of $p$.  k | 1/20 | 1/40 | 1/80 | 1/160
---|---|---|---|---
Local error | 4.9174e-5 | 6.5311e-6 | 8.3433e-7 | 1.0548e-7
Order |  | 2.91 | 2.97 | 2.98
Global error | 4.2309e-5 | 1.0515e-5 | 2.6288e-6 | 6.5912e-7
Order |  | 2.01 | 2.00 | 2.00
Table 11: Local and global error when integrating problem (37) with nonvanishing boundary conditions (39) with method (43) using $p=2$ in (29)-(30) and the simplification for the boundaries in Tables 1 and 2 for that value of $p$.  k | 1/20 | 1/40 | 1/80 | 1/160
---|---|---|---|---
Local error | 4.9884e-5 | 6.5394e-6 | 8.3434e-7 | 1.0548e-7
Order |  | 2.93 | 2.97 | 2.98
Global error | 2.0216e-4 | 4.9199e-5 | 1.2062e-5 | 2.8962e-6
Order |  | 2.04 | 2.03 | 2.06
Figure 1: Error against CPU time when integrating problem (37) with
nonvanishing boundary conditions (38), using (43) without avoiding order
reduction (green circles), the suggested technique (31)-(30) with $p=1$ and
the simplifications in Table 2 (blue asterisks), the suggested technique
corresponding to $p=2$ and the simplifications in Tables 1 and 2 (45) (magenta
squares) and the method with stiff order $2$ (49) (black diamonds). Figure 2:
Error against CPU time when integrating problem (37) with nonvanishing
boundary conditions (39), using (43) without avoiding order reduction (green
circles), the suggested technique (31)-(30) with $p=1$ and the simplifications
in Table 2 (blue asterisks), the suggested technique corresponding to $p=2$
and the simplifications in Tables 1 and 2 (45) (magenta squares) and the
method with stiff order $2$ (49) (black diamonds).

In the first place, we have considered the one-dimensional problem

$\displaystyle u_{t}=u_{xx}+u^{2}+h(x,t),\quad x\in[0,1],t\in[0,1],$ (37)

with function $h$ and initial and boundary conditions so that the exact
solution is $u(x,t)=\cos(x+t)$. On the one hand, we have considered Dirichlet
boundary conditions at both sides

$\displaystyle u(0,t)=g_{0}(t),\quad u(1,t)=g_{1}(t),$ (38)

and, on the other hand, Dirichlet at $x=0$ and Neumann at $x=1$

$\displaystyle u(0,t)=g_{0}(t),\quad u_{x}(1,t)=g_{1}(t).$ (39)

In such a way, (A1)-(A8) are clearly satisfied with $X=C[0,1]$, $\|\cdot\|$
the maximum norm (see the example in page 231 in [3]), $\phi(u)=u^{2}$,
$h(t)=-\sin(x+t)+\cos(x+t)-\cos^{2}(x+t)$, $A$ the second-order differential
operator in space, $\bar{m}(A)=2$ in (A6) and

$\|v\|_{l}=\max(\|v\|_{\infty},\|v^{\prime}\|_{\infty},\dots,\|v^{(2l)}\|_{\infty}),\quad\mbox{
for }v\in D(A^{l})=C^{2l}([0,1]),$

in (A7) and (A8).

Then, for the space discretization, we have taken the standard second-order
difference scheme, for which $D_{h}\equiv 0$ in (25) and, in the case of two
Dirichlet boundary conditions,

$A_{h,0}=\mbox{tridiag}(1,-2,1)/h^{2},\quad
C_{h}[g_{0}(t),g_{1}(t)]^{T}=[g_{0}(t),0,\dots,0,g_{1}(t)]^{T}/h^{2},$

while, in the Dirichlet and Neumann one, $A_{h,0}$ has one more row and column
since the value at the final node is also numerically approximated. More
precisely, by considering a centered difference in that last node to
approximate the Neumann condition,

$A_{h,0}=\frac{1}{h^{2}}\left(\begin{array}[]{cccccc}-2&1&0&\dots&\dots&0\\\
1&-2&1&0&\dots&\vdots\\\ 0&1&-2&1&\ddots&\vdots\\\
\vdots&0&\ddots&\ddots&\ddots&0\\\ \vdots&&\ddots&1&-2&1\\\
0&\dots&\dots&0&2&-2\end{array}\right),\quad
C_{h}\left(\begin{array}[]{c}g_{0}(t)\\\
g_{1}(t)\end{array}\right)=\left(\begin{array}[]{c}\frac{1}{h^{2}}g_{0}(t)\\\
0\\\ \vdots\\\ \frac{2}{h}g_{1}(t)\end{array}\right).$

It can be proved (see [3] and [23]) that (H1)-(H3) are satisfied for these
space discretizations when $Z=C^{4}([0,1])$, $\varepsilon_{h}=O(h^{2})$ for
the completely Dirichlet case and $\varepsilon_{h}=O(h)$ for the Dirichlet and
Neumann one, $\eta_{h}=O(h^{2})$ and $\bar{m}(Z)=4$.

As for the time integration, we have considered the EERK method

$\displaystyle\begin{array}[]{c|cc}0&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,2}&\\\ \hline\cr&0&\varphi_{1}\end{array},$
(43)

which was proved to lead to order $\geq 1$ when applied to problems with
vanishing boundary conditions [15]. (We notice that
$\varphi_{1,2}(z)=\varphi_{1}(z/2)$ and that, with the notation in (11),
$a_{21}(z)=\varphi_{1}(z/2)/2$, $b_{1}(z)=0$ and $b_{2}(z)=\varphi_{1}(z)$).
When this method is applied directly to

$\displaystyle\dot{U}_{h}(t)$ $\displaystyle=$ $\displaystyle
A_{h,0}U_{h}(t)+C_{h}g(t)+U_{h}.^{2}+P_{h}h(t),$ $\displaystyle U_{h}(0)$
$\displaystyle=$ $\displaystyle P_{h}u_{0},$ (44)

where $U_{h}.^{2}$ denotes the vector of the square of components of $U_{h}$.
We are in fact taking $M=A_{h,0}$ and $F(t,U)=C_{h}g(t)+U_{h}.^{2}+P_{h}h(t)$
in (10)-(11), and thus local and global time order around 1 can be observed
with Dirichlet boundary conditions (38), as Table 6 shows when considering
$h=1/1000$, so that the error in space is negligible. As for Dirichlet and
Neumann boundary conditions case (39), a bit higher local order can be
observed in Table 7, although the global order continues to be around 1.

As this method has nonstiff order 2, the results in this paper can be applied
for both $p=1$ and $p=2$. More precisely, when implementing the corresponding
formulas in (29) and (30) with $p=1$ and the simplification for
$\partial\bar{F}_{n,i}$ in $U_{h}^{n+1}$ in Table 2 for that value of $p$,
Theorem 3 assures local order 2 in the timestepsize when the error in space is
small enough and Theorem 4 global order 2. When implementing the same formulas
with $p=2$ using Tables 1 and 2 for the simplification of some boundaries,
Theorem 3 assures local order 3 and Theorem 4 global order 2, in principle
under a CFL condition of the form (34) with $\gamma=1$ due to Dirichlet
boundary conditions. We remind that, as stated in Section 4.5, this CFL
condition comes from the numerical approximation of $\partial A^{2}u(t_{n})$
and $\partial Af(t_{n},u(t_{n}))$ through numerical differentiation in space.
However, it happens that, for this particular method, this numerical
differentiation is not required. This is due to the fact that, for $p=2$,
(29)-(30) read like this after the simplification on the boundaries in Tables
1 and 2:

$\displaystyle K_{n,1,h}$ $\displaystyle=$ $\displaystyle U_{h}^{n},$
$\displaystyle K_{n,2,h}$ $\displaystyle=$ $\displaystyle
e^{\frac{k}{2}A_{h,0}}U_{h}^{n}+\frac{k}{2}\varphi_{1}(\frac{k}{2}A_{h,0})C_{h}\partial
u(t_{n})+\frac{k^{2}}{4}\varphi_{2}(\frac{k}{2}A_{h,0})C_{h}\partial
Au(t_{n}),$
$\displaystyle+\frac{k}{2}[\varphi_{1}(\frac{k}{2}A_{h,0})F_{n,1,h}+\frac{k}{2}\varphi_{2}(\frac{k}{2}A_{h,0})C_{h}\partial
f(t_{n},u(t_{n}))],$ $\displaystyle U_{h}^{n+1}$ $\displaystyle=$
$\displaystyle e^{kA_{h,0}}U_{h}^{n}+k\varphi_{1}(kA_{h,0})C_{h}\partial
u(t_{n})+k^{2}\varphi_{2}(kA_{h,0})C_{h}\partial
Au(t_{n})+k^{3}\varphi_{3}(kA_{h,0})C_{h}\partial A^{2}u(t_{n})$ (45)
$\displaystyle+k[\varphi_{1}(kA_{h,0})F_{n,2,h}+k\varphi_{2}(kA_{h,0})C_{h}\partial
f(t_{n}+\frac{k}{2},u(t_{n})+\frac{k}{2}\dot{u}(t_{n}))$ $\displaystyle\hskip
28.45274pt+k^{2}\varphi_{3}(kA_{h,0})C_{h}\partial Af(t_{n},u(t_{n}))].$

Then, because of (4),

$A^{2}u(t_{n})+Af(t_{n},u(t_{n}))=A\dot{u}(t_{n})=\ddot{u}(t_{n})-f_{t}(t_{n},u(t_{n}))-f_{u}(t_{n},u(t_{n}))\dot{u}(t_{n}),$

which means that the whole term which multiplies $k^{3}\varphi_{3}(kA_{h,0})$
in (45) can be exactly calculated in terms of the data of the problem, when
considering Dirichlet boundary conditions on both sides. Therefore, for this
particular problem and time integrator, numerical differentiation is not
required and condition (34) is not necessary. On the other hand, when
considering the Neumann boundary condition on one side, no space numerical
differentiation is either required, but time numerical integration is needed
indeed. For that, we have used a second-order Taylor expansion for the first
time step and then, a 2-BDF formula for the next ones.

Tables 8, 9, 10 and 11 show that the expected numerical order is in fact
achieved. When implementing the different techniques which have been described
here, putting together all the terms which multiply each $\varphi_{j}$ and
using Krylov subroutines [22], we can see that procedure (45) is the most
efficient for both types of boundary conditions. (See Figures 1 and 2, where
it is clear that, for a same amount of CPU time, the technique with the
smallest error is the one corresponding to the suggested technique with $p=2$.
A thorough explanation for that is given in [7]). Although it is not an aim of
this paper, we have also compared in the same figures with another EERK method
in [15],

$\displaystyle\begin{array}[]{c|cc}0&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,2}&\\\
\hline\cr&\varphi_{1}-2\varphi_{2}&2\varphi_{2}\end{array},$ (49)

which happens to have stiff order $2$ for problem (37). (Notice that this
method and (43) have the same underlying RK method.) When implementing this
method by integrating directly (44) with the formulas in (11), order $2$ is
observed, but is less efficient than the procedure in (45).

###  5.2 Bidimensional problem

Table 12: Local and global error when integrating Dirichlet problem (50) with nonvanishing boundary conditions with method (56) without avoding order reduction. k | 1/8 | 1/16 | 1/32 | 1/64 | 1/128 | 1/256
---|---|---|---|---|---|---
Local error | 3.6983e-5 | 4.9435e-6 | 6.6113e-7 | 8.6583e-08 | 1.1152e-8 | 1.4159e-9
Order |  | 2.90 | 2.90 | 2.93 | 2.96 | 2.98
Global error | 3.4712e-5 | 4.4791e-6 | 5.9470e-7 | 7.7735e-8 | 1.0010e-8 | 1.2770e-9
Order |  | 2.95 | 2.91 | 2.94 | 2.96 | 2.97
Table 13: Local and global error when integrating Dirichlet problem (50) with nonvanishing boundary conditions with method (56) using $p=3$ in (29)-(30) and the simplification for the boundaries in Tables 1 and 2 for that value of $p$. k | 1/8 | 1/16 | 1/32 | 1/64 | 1/128 | 1/256
---|---|---|---|---|---|---
Local error | 8.2450e-5 | 4.9697e-6 | 3.0088e-7 | 1.8147e-8 | 1.0834e-9 | 6.3605e-11
Order |  | 4.05 | 4.05 | 4.05 | 4.07 | 4.09
Global error | 1.0842e-4 | 6.9957e-6 | 4.3923e-7 | 2.7503e-8 | 1.7108e-9 | 1.0553e-10
Order |  | 3.95 | 3.99 | 4.00 | 4.01 | 4.02
Figure 3: Error against CPU time when integrating problem (50), using (56)
with the suggested technique (29)-(30) with $p=3$ and the simplifications in
Tables 1 and 2 (magenta squares), and the method in [21] with stiff order $4$
when applied directly to (57) (black diamonds).

In a second place, we have considered the bidimensional problem

$\displaystyle u_{t}=u_{xx}+u_{yy}+u^{2}+h(x,y,t),\quad x,y\in[0,1],\quad
t\in[0,1],$ (50)

with function $h$ and initial and Dirichlet boundary conditions such that the
exact solution is $u(x,y,t)=\cos(t+x+y)$ and again (A1)-(A8) are satisfied. We
notice that

$h(x,y,t)=-\sin(t+x+y)+2\cos(t+x+y)-\cos^{2}(t+x+y).$

As space discretization, we have considered the $9$-point formula [23], which
determines certain matrices $A_{h,0}$, $C_{h}$ and $D_{h}$ in (25). It can be
observed that (H1)-(H3) are again satisfied (see [2, 23]) for
$Z=C^{6}([0,1]\times[0,1])$ with $\varepsilon_{h}=O(h^{4})$,
$\eta_{h}=O(h^{4})$ and $\bar{m}(Z)=6$. As for the time integration, we
consider Krogstad’s method [17]

$\displaystyle\begin{array}[]{c|cccc}0&&&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,2}&&&\\\
\frac{1}{2}&\frac{1}{2}\varphi_{1,3}-\varphi_{2,3}&\varphi_{2,3}&&\\\
1&\varphi_{1,4}-2\varphi_{2,4}&0&2\varphi_{2,4}&\\\
\hline\cr&\varphi_{1}-3\varphi_{2}+4\varphi_{3}&2\varphi_{2}-4\varphi_{3}&2\varphi_{2}-4\varphi_{3}&4\varphi_{3}-\varphi_{2}\end{array}.$
(56)

This method was proved to lead to order $3$ when applied to stiff problems
with vanishing boundary conditions in [15], although its nonstiff order is
well known to be 4 [17]. We have implemented the method when integrating
directly

$\displaystyle\dot{U}_{h}(t)$ $\displaystyle=$ $\displaystyle
A_{h,0}U_{h}(t)+C_{h}g(t)+U_{h}.^{2}+P_{h}h(t)+D_{h}[(g(t)^{2}+\partial
h(t)-\dot{g}(t)],$ $\displaystyle U_{h}(0)$ $\displaystyle=$ $\displaystyle
P_{h}u(0),$ (57)

through formulas (10)-(11) and Krylov methods [22] to calculate the
exponential operators (We notice that, as the matrix $A_{h,0}$ is not so
sparse in this case because it is just the product of the inverse of a sparse
mass matrix times another sparse matrix, a modification of the subroutines in
[22] has been used following the lines of Section 7.1 in [12], so as to take
profit of the highly sparse structure of the underlying matrices.) Local and
global order $3$ then turn up, as Table 9 shows for the space grid $h=1/160$.
We have then applied the suggested formulas (29)-(30) with $p=3$ with this
method, considering the simplifications in Tables 1 and 2. We notice that, for
this particular method, no numerical differentiation is required for the
stages since the terms in $\partial A^{2}u(t_{n})$ and $\partial
A\bar{F}_{n,i}$ sum up to a term in $\partial A\dot{u}$, which can be
calculated in terms of data because
$A\dot{u}=\ddot{u}-\dot{h}-\phi^{\prime}(u)\dot{u}$. However, for the
approximation of $\partial A^{2}u(t_{n})/\partial\bar{F}_{n,i}$ in
$U_{h}^{n+1}$, $u_{x}$ and $u_{y}$ are required at the boundary, for which we
have used a 4th-order BDF formula in the sides of the square in which we could
not calculate it directly in terms of $g$. As for the approximation of
$\partial A^{3}u(t_{n})/\partial A\bar{F}_{n,i}/\partial A^{2}\bar{F}_{n,i}$,
we required not only $u_{x}$, $u_{y}$, but also $\dot{u}_{x}$, $\dot{u}_{y}$
at the boundary. To calculate them, we have firstlyapproximated $\dot{u}$ in
some points near the boundary, but at the interior of the domain. We have done
that by using Taylor expansions of order $3$ near the initial condition, i.e.,

$u_{t}(x,y,lk)\approx
u_{t}(x,y,0)+lku_{tt}(x,y,0)+\frac{(lk)^{2}}{2}u_{ttt}(x,y,0),\quad l=0,1,2,$

where the terms above can be calculated in terms of $u_{0}$ by using equation
(50). Then, after the third step, we have used the 3rd-order BDF formula

$u_{t}(x,y,lk)\approx\frac{1}{k}[\frac{11}{6}u(x,y,lk)-3u(x,y,(l-1)k)+\frac{3}{2}u(x,y,(l-2)k)-\frac{1}{3}u(x,y,(l-3)k)].$

Once those approximations of $\dot{u}$ were calculated, we have approximated
$\dot{u}_{x}$, $\dot{u}_{y}$ with a 4th-order BDF formula in space. In such a
way, $\mu_{k,1}=O(k^{3})$ and $\nu_{h}=O(h^{4})$ in Theorem 4, so that the
error in space is easy to be negligible and the fourth order in time can be
seen. In particular, condition (34) does not seem very restrictive since local
and global order 4 turn up, as it was justified by Theorems 3 and 4, and Table
10 corroborates.

Again, although it is not an aim of this paper, we have made a comparison in
CPU time of this method with another one which shows order 4 without avoiding
order reduction. We firstly tested with method (5.19) in [15] because it was
suggested there as a $4$th-order method for stiff problems with vanishing
boundary conditions, but it happened to show just order $3$ with non-vanishing
boundary values. Because of that, we decided to consider the $5$th-order
method for stiff problems and vanishing boundary conditions which is
constructed in [21]. This method is $4$th-order accurate when applied directly
to the space discretized system (57). The comparison in CPU time is shown in
Figure 3, where it is shown that the latter is about $1/2$ more expensive than
the suggested modification of Krogstad’s method when avoiding order reduction.

## Acknowledgements

This research has been supported by Ministerio de Ciencia e Innovación and
Regional Development European Funds through project PGC2018-101443-B-I00 and
by Junta de Castilla y León and Feder through project VA169P20.

## References

  * [1] I. Alonso–Mallo, B. Cano and N. Reguera, Analysis of order reduction when integrating linear initial boundary value problems with Lawson methods, Applied Numerical Mathematics 118, pp. 64-74 (2017).
  * [2] I. Alonso–Mallo, B. Cano and N. Reguera, Avoiding order reduction when integrating linear initial boundary value problems with Lawson methods, IMA Journal of Numerical Analysis 37 (2017), pp. 2091–2119.
  * [3] I. Alonso–Mallo, B. Cano and N. Reguera, Avoiding order reduction when integrating reaction-diffusion boundary value problems with exponential splitting methods, J. Comput. Appl. Math. 357 (2019), 228–250.
  * [4] B. Cano and M. J. Moreta, Exponential quadrature rules without order reduction for integrating linear initial boundary value problems, SIAM J. Num. Anal. 56-3 (2018), 1187–1209.
  * [5] B. Cano and N. Reguera, CMMSE: Analysis of order reduction when Lawson methods integrate nonlinear initial boundary value problems, Mathematical Methods in the Applied Sciences (2022), https://doi.org/10.1002/mma.8451
  * [6] B. Cano and N. Reguera, How to avoid order reduction when Lawson methods integrate nonlinear initial boundary value problems, BIT Numerical Mathematics 62 (2022), 431–463.
  * [7] B. Cano and N. Reguera, Why improving the accuracy of exponential integrators can decrease their computational cost?, Mathematics 9(9) (2021), 1008.
  * [8] J. Connors, J.W. Banks, J.A. Hittinger, and C.S. Woodward, Quantification of errors for operator-split advection-diffusion calculations, Comput. Methods Appl. Mech. Engrg. 272 (2014), 181–197.
  * [9] L. Einkemmer and A. Ostermann, Overcoming order reduction in diffusion-reaction splitting. Part 1: Dirichlet boundary conditions, SIAM J. Sci. Comput. 37 (3) (2015), A1577–A1592.
  * [10] L. Einkemmer and A. Ostermann, Overcoming order reduction in diffusion-reaction splitting. Part 2: Oblique boundary conditions, SIAM J. Sci. Comput. 38 (2016) A3471-A3757.
  * [11] E. Faou, A. Ostermann and K. Schratz, Analysis of exponential splitting methods for inhomogenoous parabolic equations, IMA J. Numer. Anal. 35 (2015), 161–178.
  * [12] T. Göckler and V. Grimm, Convergence analysis of an extended Krylov subspace method for the approximation of operator functions in exponential integratiors, SIAM J. Numer. Anal. 53 (4) (2013) 2189–2213.
  * [13] E. Hairer, S. Nörsett and G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems, Second Revised Edition, Springer, 2000.
  * [14] M. Hochbruck, A. Ostermann, Exponential Runge-Kutta methods for parabolic problems, Appl. Numer. Math. 53 (2-4) (2005) 323–339.
  * [15] M. Hochbruck and A. Ostermann, Explicit exponential Runge-Kutta methods for semilinear parabolic problems, SIAM J. Num. Anal. 43 (2005), 1069–1090.
  * [16] M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numerica (2010), 209-286.
  * [17] S. Krogstad, Generalized integrating factor methods for stiff PDEs, J. Comp. Phys. 203 (2005), 72–88.
  * [18] J. D. Lawson, Generalized Runge-Kutta processes for stable systems with large Lipschitz constants, SIAM J. Numer. Anal. 4 (1967), 372–380.
  * [19] R. J. LeVeque and J. Oliger, Numerical methods based on additive splittings for hyperbolic partial differential equations, Math. Comp. 40 (1983), no. 162, 469–497.
  * [20] V.T. Luan, A. Ostermann, Stiff order conditions for exponential Runge-Kutta methods of order five, Modeling, Simulation and Optimization of Complex Processes-HPSC (2012), 133–143.
  * [21] V.T. Luan, A. Ostermann, Explicit exponential Runge-Kutta methods of high order for parabolic problems, J. Comput. Appl. Math. 262 (2014), 361–372.
  * [22] J. Niesen and W. M. Wright, Algorithm 919: a Krylov subspace algorithm for evaluating the $\varphi$-functions appearing in exponential integrators, ACM Trans. Math. Software 38, no. 3, Art. 22 (2012).
  * [23] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wadsworth & Brooks, United States of America (1989).

##  6 Appendix

###  6.1 Proof of Theorem 1

Notice that, for the first set of hypotheses, through its definition and using
(15)-(16),

$\displaystyle\bar{u}_{n+1}$ $\displaystyle=$
$\displaystyle\phi_{0,p,u(t_{n}),u(t_{n})}(k)+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(k)$
(58) $\displaystyle=$
$\displaystyle\sum_{l=0}^{p}\frac{k^{l}}{l!}A^{l}u(t_{n})+k^{p+1}\varphi_{p+1}(kA_{0})A^{p+1}u(t_{n})$
$\displaystyle+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}[\sum_{j=0}^{p-1}\frac{k^{j}}{(j+l)!}A^{j}\bar{F}_{n,i}+k^{p}\varphi_{p+l}(kA_{0})A^{p}\bar{F}_{n,i}],$

where $\bar{F}_{n,i}=f(t_{n}+c_{i}k,\bar{K}_{n,i})$, with

$\displaystyle\bar{K}_{n,i}$ $\displaystyle=$
$\displaystyle\phi_{0,p-1,u(t_{n}),u(t_{n})}(c_{i}k)+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\phi_{l,p-2,\bar{F}_{n,j},\bar{F}_{n,j}}(c_{r}k)$
(59) $\displaystyle=$
$\displaystyle\sum_{l=0}^{p-1}\frac{(c_{i}k)^{l}}{l!}A^{l}u(t_{n})+(c_{i}k)^{p}\varphi_{p}(c_{i}kA_{0})A^{p}u(t_{n})$
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}[\sum_{ll=0}^{p-2}\frac{(c_{r}k)^{ll}}{(ll+l)!}A^{ll}\bar{F}_{n,j}+(c_{r}k)^{p-1}\varphi_{p-1+l}(c_{r}kA_{0})A^{p-1}\bar{F}_{n,j}].$

In the following, we will inductively prove that $\bar{K}_{n,i}\in D(A^{p})$
and $A^{p}\bar{K}_{n,i}$ is uniformly bounded in $k$ for small enough $k$, so
that (58) and (59) have perfect sense because of (A6) and the assumed
regularity on $\phi$ and $h$. The result comes from Lemma 3.1 in [2] because,
according to it,

$\displaystyle
A\phi_{0,p-1,u(t_{n}),u(t_{n})}(\tau)=\sum_{l=0}^{p-2}\frac{\tau^{l}}{l!}A^{l+1}u(t_{n})+\tau^{p-1}\varphi_{p-1}(\tau
A_{0})A^{p}u(t_{n}),$ (60)

and, arguing in the same way,

$\displaystyle A^{2}\phi_{0,p-1,u(t_{n}),u(t_{n})}(\tau)$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{p-3}\frac{\tau^{l}}{l!}A^{l+2}u(t_{n})+\tau^{p-2}\varphi_{p-2}(\tau
A_{0})A^{p}u(t_{n}),$ $\displaystyle\vdots$ $\displaystyle=$
$\displaystyle\vdots$ $\displaystyle
A^{p}\phi_{0,p-1,u(t_{n}),u(t_{n})}(\tau)$ $\displaystyle=$
$\displaystyle\varphi_{0}(\tau A_{0})A^{p}u(t_{n}).$ (61)

From this, $\bar{K}_{n,1}\in D(A^{p})$, which implies that $\bar{F}_{n,1}\in
D(A^{p})$. As for $\bar{K}_{n,2}$, in a similar way we notice that

$\displaystyle A\phi_{j,p-2,\bar{F}_{n,1},\bar{F}_{n,1}}(\tau)$
$\displaystyle=$
$\displaystyle\sum_{l=0}^{p-3}\frac{\tau^{l}}{(l+j)!}A^{l+1}\bar{F}_{n,1}+\tau^{p-2}\varphi_{j+p-2}(\tau
A_{0})A^{p-1}\bar{F}_{n,1},$ $\displaystyle\vdots$ $\displaystyle=$
$\displaystyle\vdots$ $\displaystyle
A^{p-1}\phi_{j,p-2,\bar{F}_{n,1},\bar{F}_{n,1}}(\tau)$ $\displaystyle=$
$\displaystyle\varphi_{j}(\tau A_{0})A^{p-1}\bar{F}_{n,1},$ $\displaystyle
A^{p}\phi_{j,p-2,\bar{F}_{n,1},\bar{F}_{n,1}}(\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}[\varphi_{j+1}(\tau
A_{0})-\frac{1}{(j+1)!}I]A^{p-1}\bar{F}_{n,1},\quad\mbox{ for }\tau\neq 0,$
(62)

and therefore,

$kA^{p}\phi_{j,p-2,\bar{F}_{n,1},\bar{F}_{n,1}}(c_{r}k)=\bigg{\\{}\begin{array}[]{ll}\frac{k}{j!}A^{p}\bar{F}_{n,1},&\mbox{
if }c_{r}=0,\\\
\frac{1}{c_{r}}[\varphi_{j+1}(c_{r}kA_{0})-\frac{1}{(j+1)!}I]A^{p-1}\bar{F}_{n,1},&\mbox{
if }c_{r}\neq 0,\end{array}$

from what $\bar{K}_{n,2}\in D(A^{p})$ and $A^{p}\bar{K}_{n,2}$ is uniformly
bounded on $k$ for small enough $k$, which implies the same fact for
$\bar{F}_{n,2}$. Proceeding inductively with the rest of stages, it is
justified that (58) and (59) are true.

Then, developing in Taylor series (which can be done because of (A7) and the
assumed regularity), the terms in $k^{j}$ ($0\leq j\leq p-1$) in (59) have the
same form as those in (23) but substituting $M$ by the operator $A$,
$U(t_{n})$ by $u(t_{n})$ and $F_{*,n}$ by $f_{*}(t_{n},u(t_{n}))$. Moreover,
the terms in $k^{j}(0\leq j\leq p)$ in (58) have the same form as those in
(22) with the same type of substitution. We notice that the residue terms in
the Taylor expansions around $u(t_{n})$ of $A^{j}\phi(\bar{K}_{n,i})$
($j=0,\dots,p-1$) in (58) and (59) correspond to

$\frac{1}{(p-j)!}\int_{0}^{1}A^{j}\big{[}\phi^{(p-j)}\big{(}(1-\lambda)u(t_{n})+\lambda\bar{K}_{n,i}\big{)}(\bar{K}_{n,i}-u(t_{n}))^{p-j}\big{]}d\lambda.$

Using that $(1-\lambda)u(t_{n})+\lambda\bar{K}_{n,i}\in D(A^{j})$ for
$\lambda\in[0,1]$ and that $\|\bar{K}_{n,i}-u(t_{n})\|_{j}=O(k)$, as justified
through (60),(61) and (62), from (A7) it follows that these residues are
$O(k^{p-j})$ whenever $\phi\in C^{p-j+j\bar{m}(A)}$. We notice that this
regularity of $\phi$ for $j=0,\dots,p-1$ is satisfied under the assumed
regularity. On the other hand, the factor $k^{j+1}$ which is multiplying
$A^{j}\bar{F}_{n,i}$ in (58) and (59) leads to the fact that all the residues
are $O(k^{p+1})$.

At the same time, the same happens for the Taylor series of $u(t_{n+1})$,
which terms in $k^{j}$ ($0\leq j\leq p$) also coincide with those in (21)
except for the previous substitution of terms. As those terms coincided for an
ordinary differential system, they also coincide now till the terms in which
we are interested.

On the other hand, if in formulas (17)-(19), $p$ is substituted by
$\hat{p}<p$, the same argument leads to the fact that
$\rho_{n}=O(k^{\hat{p}+1})$. Moreover, when applying $A_{0}^{-1}$, the first
term in $k^{\hat{p}+1}$ in (58) would be

$\displaystyle
A_{0}^{-1}k^{\hat{p}+1}\varphi_{\hat{p}+1}(kA_{0})A^{\hat{p}+1}u(t_{n})$
$\displaystyle=$ $\displaystyle
A_{0}^{-1}k^{\hat{p}+1}[kA_{0}\varphi_{\hat{p}+2}(kA_{0})+\frac{1}{(\hat{p}+1)!}]A^{\hat{p}+1}u(t_{n})$
$\displaystyle=$ $\displaystyle
k^{\hat{p}+2}\varphi_{\hat{p}+2}(kA_{0})A^{\hat{p}+1}u(t_{n})+\frac{k^{\hat{p}+1}}{(\hat{p}+1)!}A_{0}^{-1}A^{\hat{p}+1}u(t_{n}),$

where the formula (8) has been used in the first equality. We notice that the
first term is $O(k^{\hat{p}+2})$ and the second one coincides with one
corresponding to the power $k^{\hat{p}+1}$ in $A_{0}^{-1}u(t_{n+1})$. In a
similar way, the other terms in $k^{\hat{p}+1}$ in $A_{0}^{-1}$ times (58)
would be

$\displaystyle\mu_{i,l}A_{0}^{-1}k^{\hat{p}+1}\varphi_{\hat{p}+l}(kA_{0})A^{\hat{p}}\bar{F}_{n,i}^{\hat{p}}=\mu_{i,l}\bigg{[}k^{\hat{p}+2}\varphi_{\hat{p}+l+1}(kA_{0})A^{\hat{p}}\bar{F}_{n,i}^{\hat{p}}+\frac{k^{\hat{p}+1}}{(\hat{p}+l)!}A_{0}^{-1}A^{\hat{p}}\bar{F}_{n,i}^{\hat{p}}\bigg{]},$

where $\bar{F}_{n,i}^{\hat{p}}=f(t_{n}+c_{i}k,\bar{K}_{n,i}^{\hat{p}})$ with
$\bar{K}_{n,i}^{\hat{p}}$ that in (59) with $p$ substituted by $\hat{p}$. We
notice again that the first term is $O(k^{\hat{p}+2})$ and the second
coincides with one of the others corresponding to the power $k^{\hat{p}+1}$ in
$A_{0}^{-1}u(t_{n+1})$ considering that the difference between $\bar{K}_{n,i}$
and $\bar{K}_{n,i}^{\hat{p}}$ is at least $O(k)$ (due to $\hat{p}\geq 1$),
hypothesis (A4) and taking (22) into account and the fact that the terms of
that type there correspond to the same terms in $U(t_{n+1})$.

###  6.2 Proof of Theorem 2

We notice that $\rho_{n,h}$ can be decomposed as

$\rho_{n,h}=(\bar{U}_{h}^{n+1}-P_{h}\bar{u}_{n+1})+P_{h}\rho_{n},$

where $\rho_{n}$ is the local error associated to the time semidiscretization
and can be bounded as Theorem 1 states. Using (24) and (H3d), it suffices to
prove then that $\bar{U}_{h}^{n+1}-P_{h}\bar{u}_{n+1}=O(k\varepsilon_{h})$.
For that, we firstly take into account that

$\displaystyle\bar{U}_{h}^{n+1}-P_{h}\bar{u}_{n+1}$ $\displaystyle=$
$\displaystyle\phi_{0,p,P_{h}u(t_{n}),u(t_{n})}^{h}(k)-P_{h}\phi_{0,p,u(t_{n}),u(t_{n})}(k)$
(63)
$\displaystyle+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}[\phi_{l,p-1,\bar{F}_{n,i,h},\bar{F}_{n,i}}^{h}(k)-P_{h}\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(k)].$

By definition of $\phi_{0,r,P_{h}u(t_{n}),u(t_{n})}^{h}$ and
$\phi_{0,r,u(t_{n}),u(t_{n})}$ and using (27), it follows that

$\displaystyle\frac{d}{d\tau}[\phi_{0,p,P_{h}u(t_{n}),u(t_{n})}^{h}(\tau)-P_{h}\phi_{0,p,u(t_{n}),u(t_{n})}(\tau)]$
$\displaystyle=$ $\displaystyle
A_{h,0}[\phi_{0,p,P_{h}u(t_{n}),u(t_{n})}^{h}(\tau)-P_{h}\phi_{0,p,u(t_{n}),u(t_{n})}(\tau)]$
$\displaystyle+A_{h,0}(P_{h}-R_{h})\phi_{0,p,u(t_{n}),u(t_{n})}(\tau)$
$\displaystyle\phi_{0,p,P_{h}u(t_{n}),u(t_{n})}^{h}(0)-P_{h}\phi_{0,p,u(t_{n}),u(t_{n})}(0)$
$\displaystyle=$ $\displaystyle 0,$

which implies that

$\displaystyle\phi_{0,p,P_{h}u(t_{n}),u(t_{n})}^{h}(k)-P_{h}\phi_{0,p,u(t_{n}),u(t_{n})}(k)$
(64) $\displaystyle=$
$\displaystyle\int_{0}^{k}e^{(k-s)A_{h,0}}A_{h,0}(P_{h}-R_{h})\phi_{0,p,u(t_{n}),u(t_{n})}(s)ds=O(k\varepsilon_{h}),$

by using (28) and that $\phi_{0,p,u(t_{n}),u(t_{n})}\in C([0,T],Z)$ because of
(15), (32) and (H2a). In a similar way,

$\displaystyle\frac{d}{d\tau}[\phi_{l,p-1,\bar{F}_{n,i,h},\bar{F}_{n,i}}^{h}(\tau)-P_{h}\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(\tau)]=(A_{h,0}-\frac{l}{\tau}I_{h})[\phi_{l,p-1,\bar{F}_{n,i,h},\bar{F}_{n,i}}^{h}(\tau)-P_{h}\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(\tau)]$
$\displaystyle\hskip
85.35826pt+\frac{1}{(l-1)!\tau}[\bar{F}_{n,i,h}-P_{h}\bar{F}_{n,i}]+A_{h,0}(P_{h}-R_{h})\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(\tau)$
$\displaystyle\phi_{l,p-1,\bar{F}_{n,i,h},\bar{F}_{n,i}}^{h}(0)-P_{h}\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(0)=\frac{1}{l!}[\bar{F}_{n,i,h}-P_{h}\bar{F}_{n,i}],$

where $I_{h}$ is the identity matrix with the same size as $A_{h,0}$. Now, we
take into account that the solution of this system is
$V_{h}(\tau)+W_{h}(\tau)$, with $V_{h}(\tau)$ the solution of

$\displaystyle V_{h}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A_{h,0}-\frac{l}{\tau}I_{h})V_{h}(\tau)+A_{h,0}(P_{h}-R_{h})\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(\tau)$
$\displaystyle V_{h}(0)$ $\displaystyle=$ $\displaystyle 0,$

and $W_{h}$ the solution of

$\displaystyle W_{h}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A_{h,0}-\frac{l}{\tau}I_{h})W_{h}(\tau)+\frac{1}{(l-1)!\tau}[\bar{F}_{n,i,h}-P_{h}\bar{F}_{n,i}]$
$\displaystyle W_{h}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\bar{F}_{n,i,h}-P_{h}\bar{F}_{n,i}].$

From this,

$\displaystyle V_{h}(k)$ $\displaystyle=$
$\displaystyle\int_{0}^{k}e^{\int_{s}^{k}(A_{h,0}-\frac{l}{\tau}I)d\tau}A_{h,0}(P_{h}-R_{h})\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(s)ds$
(65) $\displaystyle=$
$\displaystyle\int_{0}^{k}e^{(k-s)A_{h,0}}(\frac{s}{k})^{l}A_{h,0}(P_{h}-R_{h})\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}(s)ds$
$\displaystyle=$ $\displaystyle O(k\varepsilon_{h}),$

where again (28) has been used noting that, inductively,
$\phi_{l,p-1,\bar{F}_{n,i},\bar{F}_{n,i}}\in C([0,T],Z)$ because of (16),
(32), the regularity of $\phi$, (61), (62) and (H2a). On the other hand, we
notice that, using Lemma 6 in [4],

$\displaystyle
W_{h}(k)=\varphi_{l}(kA_{h,0})[\bar{F}_{n,i,h}-P_{h}\bar{F}_{n,i}].$ (66)

As $\varphi_{l}(kA_{h,0})$ is bounded because of (H1b) and, because of (26),
$P_{h}\bar{F}_{n,i}=f(t_{n}+c_{i}k,P_{h}\bar{K}_{n,i})$, using (H3), what we
need is to bound

$\displaystyle\bar{K}_{n,i,h}-P_{h}\bar{K}_{n,i}$ $\displaystyle=$
$\displaystyle\phi_{0,p-1,P_{h}u(t_{n}),u(t_{n})}^{h}(c_{i}k)-P_{h}\phi_{0,p-1,u(t_{n}),u(t_{n})}(c_{i}k)$
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}[\phi_{l,p-2,\bar{F}_{n,j,h},\bar{F}_{n,j}}^{h}(c_{r}k)-P_{h}\phi_{l,p-2,\bar{F}_{n,j},\bar{F}_{n,j}}(c_{r}k)].$

Using the same arguments as above, it can be proved, recursively in
$i=1,\dots,s$, that $\bar{K}_{n,i,h}-P_{h}\bar{K}_{n,i}=O(k\varepsilon_{h})$,
from what, considering that in (66) and inserting it in (63) together with
(65) and (64), the first result follows.

As for the result concerning the second set of hypotheses, it suffices to
decompose again

$A_{h,0}^{-1}\rho_{n,h}=A_{h,0}^{-1}(\bar{U}_{n+1,h}-P_{h}\bar{u}_{n+1})+A_{h,0}^{-1}P_{h}\rho_{n}.$

Then, in a similar way as above, but using now the second formula of (28) when
possible, it can be proved that
$A_{h,0}^{-1}(\bar{U}_{n+1,h}-P_{h}\bar{u}_{n+1})=O(k\eta_{h}+k^{2}\varepsilon_{h})$.
On the other hand, as $\partial A_{0}^{-1}\rho_{n+1}=0$, by using (27) and the
fact that, at the next step, the boundary of the numerical solution will be
considered as the exact, $A_{h,0}R_{h}(A_{0}^{-1}\rho_{n})=P_{h}\rho_{n}$,
from what

$A_{h,0}^{-1}P_{h}\rho_{n}=P_{h}(A_{0}^{-1}\rho_{n})+(R_{h}-P_{h})A_{0}^{-1}\rho_{n}.$

Now, the first term on the right-hand side is $O(k^{\hat{p}+2})$ by using the
second part of Theorem 1 and (24). As for the second term, it is
$O(k\eta_{h})$ because, for some $t^{*}\in(t_{n},t_{n+1})$,
$(u(t_{n+1})-u(t_{n}))/k=\dot{u}(t^{*})=Au(t^{*})+f(t^{*},u(t^{*}))\in Z$ and
$(\bar{u}_{n+1}-u(t_{n}))/k\in Z$, due to (32) and (H2a). Therefore,
$A_{0}^{-1}\rho_{n}/k\in Z$ and (28) can be applied.

###  6.3 Proof of Theorem 3

We firstly consider $p=1$ and notice that the same decomposition as in the
proof of Theorem 2 can be performed, i.e,

$\displaystyle\rho_{n,h}^{simp}=(\bar{U}_{h}^{n+1,simp}-P_{h}\bar{u}_{n+1}^{simp})+P_{h}\rho_{n}^{simp},$
(67)

where $\bar{U}_{h}^{n+1,simp}$ is that corresponding to (30) with $U_{h}^{n}$
substituted by $P_{h}u(t_{n})$, $p=1$ and the simplified value
$\partial\bar{F}_{n,i}=\partial f(t_{n},u(t_{n}))$. On the other hand,

$\bar{u}_{n+1}^{simp}=\phi_{0,1,u(t_{n}),u(t_{n})}(k)+k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\phi_{l,0,\bar{F}_{n,i},\bar{F}_{n,i}^{u,1,simp}}(k),$

where $\bar{F}_{n,i}^{u,1,simp}=f(t_{n},u(t_{n}))$, and then
$\rho_{n}^{simp}=\bar{u}_{n+1}^{simp}-u(t_{n+1})$. (We notice that, for $p=1$,
there is no difference in the required boundaries for the stages and therefore
$\bar{F}_{n,i}$ is the same which would correspond to $\bar{u}_{n+1}$.)

Then, in (67), the first parenthesis can be bounded in the same way than in
the proof of Theorem 2 and, in order to bound the second term, it suffices to
consider that

$\rho_{n}^{simp}=(\bar{u}_{n+1}^{simp}-\bar{u}_{n+1})+\rho_{n},$

that $\rho_{n}$ can be bounded as in Theorem 1 with $p=1$ and that, using
(16),

$\displaystyle\bar{u}_{n+1}-\bar{u}_{n+1}^{simp}$ $\displaystyle=$
$\displaystyle
k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\big{[}\frac{1}{l!}[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]$
(68) $\displaystyle\hskip
56.9055pt-\varphi_{l}(kA_{0})[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]$
$\displaystyle\hskip
56.9055pt+k\varphi_{l+1}(kA_{0})A[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]\big{]}.$

Then, by (59), $\|\bar{K}_{n,i}-u(t_{n})\|=O(k)$ and therefore, because of the
assumed regularity of $\phi$, $h$ and (A7), the term in small brackets in (68)
is $O(k)$ and $A$ times that bracket is at least bounded. From this, it is
clear that $\|\bar{u}_{n+1}-\bar{u}_{n+1}^{simp}\|=O(k^{2})$, so that the
first part of the theorem follows.

As for the second part of the theorem when $p=1$, it suffices to follow the
same arguments as in the proof of Theorem 2 again, taking into account that
$A_{0}^{-1}\rho_{n}^{simp}=O(k^{3})$ because, using (8),

$\displaystyle A_{0}^{-1}(\bar{u}_{n+1}-\bar{u}_{n+1}^{simp})$
$\displaystyle=$ $\displaystyle
k\sum_{i=1}^{s}\sum_{l=1}^{s}\mu_{i,l}\bigg{[}-k\varphi_{l+1}(kA_{0})[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]$
$\displaystyle\hskip
28.45274pt+k^{2}\varphi_{l+2}(kA_{0})A[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]$
$\displaystyle\hskip
28.45274pt+\frac{k}{(l+1)!}A_{0}^{-1}A[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]\bigg{]}.$

It suffices to notice that $\bar{K}_{n,i}-u(t_{n})\in D(A_{0})$ and because of
(A8) and (A7),

$\displaystyle\|A_{0}^{-1}A[f(t_{n}+c_{i}k,\bar{K}_{n,i})-f(t_{n},u(t_{n}))]\|$
$\displaystyle=\|\int_{0}^{1}A_{0}^{-1}A[\phi^{\prime}(\lambda\bar{K}_{n,i}+(1-\lambda)u(t_{n}))(\bar{K}_{n,i}-u(t_{n}))]d\lambda+c_{i}kA_{0}^{-1}A\dot{h}(t^{*})\|=O(k),$

from what the result follows.

As for $p=2$, by using again (67), the key point is to bound appropriately
$\rho_{n}^{simp}$ under the first new set of hypotheses and
$A_{0}^{-1}\rho_{n}^{simp}$ under the second set. However, now the stages are
calculated differently than in (23) since the approximation
$\partial\bar{F}_{n,i}\approx\partial\bar{F}_{n,i}^{K,2,simp}$ at each stage
with $\bar{F}_{n,i}^{K,2,simp}=f(t_{n},u(t_{n}))$ must be used. More
precisely,

$\displaystyle\bar{K}_{n,i}^{simp}$ $\displaystyle=$
$\displaystyle\phi_{0,1,u(t_{n}),u(t_{n})}(c_{i}k)$ (69)
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\phi_{l,0,\bar{F}_{n,j}^{simp},\bar{F}_{n,j}^{K,2,simp}}(c_{r}k),\quad
i=1,\dots,s,$

with $\bar{F}_{n,j}^{simp}=f(t_{n}+c_{j}k,\bar{K}_{n,j}^{simp})$. On the other
hand, as for the definition of $\bar{u}_{n+1}^{simp}$, the approximation for
the boundaries in Table 2 must be considered, denoting

$\displaystyle\bar{F}_{n,i}^{u,2,simp}$ $\displaystyle=$ $\displaystyle
f\big{(}t_{n}+c_{i}k,u(t_{n})+k[c_{i}Au(t_{n})+\sum_{j,l,r}\frac{\lambda_{i,j,l,r}}{l!}f(t_{n},u(t_{n}))]\big{)},$
$\displaystyle\bar{F}_{n,i}^{u,2,A,simp}$ $\displaystyle=$ $\displaystyle
f(t_{n},u(t_{n})),$

it follows that

$\displaystyle\bar{u}_{n+1}-\bar{u}_{n+1}^{simp}=k\sum_{i,l}\mu_{i,l}\Delta_{i,l}(k),$
(70)

where

$\displaystyle\Delta_{i,l}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A-\frac{l}{\tau}I)\Delta_{i,l}(\tau)+\frac{1}{(l-1)!\tau}[\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}],$
$\displaystyle\Delta_{i,l}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}],$
$\displaystyle\partial\Delta_{i,l}(\tau)$ $\displaystyle=$
$\displaystyle\partial\bigg{[}\frac{1}{l!}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,simp}]+\frac{\tau}{(l+1)!}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,A,simp}]\bigg{]}.$
(71)

From here, denoting by $W_{i,l}(\tau)$ to the term in big brackets in (71),

$\displaystyle(\Delta_{i,l}-W_{i,l})^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A-\frac{l}{\tau}I)(\Delta_{i,l}-W_{i,l})(\tau)+\frac{1}{(l-1)!\tau}[\bar{F}_{n,i}^{u,2,simp}-\bar{F}_{n,i}^{simp}]$
$\displaystyle+\frac{1}{l!}A[\bar{F}_{n,i}^{u,A,2,simp}-\bar{F}_{n,i}^{u,2,simp}]+\frac{\tau}{(l+1)!}A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,simp}],$
$\displaystyle(\Delta_{i,l}-W_{i,l})(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\bar{F}_{n,i}^{u,2,simp}-\bar{F}_{n,i}^{simp}],$
$\displaystyle\partial(\Delta_{i,l}-W_{i,l})(\tau)$ $\displaystyle=$
$\displaystyle 0,$

so that

$\displaystyle(\Delta_{i,l}-W_{i,l})(k)$ $\displaystyle=$
$\displaystyle\varphi_{l}(kA_{0})[\bar{F}_{n,i}^{u,2,simp}-\bar{F}_{n,i}^{simp}]+k\varphi_{l+1}(kA_{0})A[\bar{F}_{n,i}^{u,2,A,simp}-\bar{F}_{n,i}^{u,2,simp}]$
(72)
$\displaystyle+k^{2}\varphi_{l+2}(kA_{0})A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,A,simp}].$

Now, from (69), (15), (16), the assumed regularity and (A7),
$\|\bar{F}_{n,i}^{u,2,simp}-\bar{F}_{n,i}^{simp}\|=O(k^{2})$. On the other
hand, it is also clear that
$\|A[\bar{F}_{n,i}^{u,A,2,simp}-\bar{F}_{n,i}^{u,2,simp}]\|=O(k)$. Therefore,
from (72), $\|(\Delta_{i,l}-W_{i,l})(k)\|=O(k^{2})$. Using also that, from
(59), $\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,s}=O(k^{2})$ and through (A7),
$\|A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,A,simp}]\|=O(k)$, it follows that
$W_{i,l}(k)$ is also $O(k^{2})$. Then,

$\|\Delta_{i,l}(k)\|=\|(\Delta_{i,l}-W_{i,l})(k)+W_{i,l}(k)\|=O(k^{2}),$

and from (70) and the first part of Theorem 1 for $p=2$,
$\rho_{n}^{simp}=O(k^{3})$, which proves the first part of the present
theorem.

In order to prove the second part of the theorem for $p=2$, let us first
notice that, using (16),

$\displaystyle\bar{K}_{n,i}-\bar{K}_{n,i}^{simp}$ $\displaystyle=$
$\displaystyle
k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\big{[}\frac{1}{l!}[\bar{F}_{n,j}-\bar{F}_{n,j}^{K,2,simp}]-\varphi_{l}(c_{r}kA_{0})[\bar{F}_{n,j}^{simp}-\bar{F}_{n,j}^{K,2,simp}]$
$\displaystyle\hskip
85.35826pt+c_{r}k\varphi_{l+1}(c_{r}kA_{0})A[\bar{F}_{n,j}-\bar{F}_{n,j}^{K,2,simp}]\big{]}.$

Inductively on the stages, it is therefore proved that
$\|\bar{K}_{n,i}-\bar{K}_{n,i}^{simp}\|=O(k^{2})$, from what
$\|\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}\|=O(k^{2})$. But, furthermore, using
(8),

$\displaystyle A_{0}^{-1}(\bar{K}_{n,i}-\bar{K}_{n,i}^{simp})$
$\displaystyle=$ $\displaystyle
k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\bigg{[}\frac{1}{l!}A_{0}^{-1}[\bar{F}_{n,j}-\bar{F}_{n,j}^{simp}]-k\varphi_{l+1}(c_{r}kA_{0})[\bar{F}_{n,j}^{simp}-\bar{F}_{n,j}^{K,2,simp}]$
$\displaystyle\hskip
14.22636pt+(c_{r}k)^{2}\varphi_{l+2}(c_{r}kA_{0})A[\bar{F}_{n,j}-\bar{F}_{n,j}^{K,2,simp}]+\frac{c_{r}k}{(l+1)!}A_{0}^{-1}A[\bar{F}_{n,j}-\bar{F}_{n,j}^{K,2,simp}]\bigg{]},$

which implies, using (A7) and (A8), that
$\|A_{0}^{-1}(\bar{K}_{n,i}-\bar{K}_{n,i}^{simp})\|=O(k^{3})$. Now, we notice
that, from (72) and the definition of $W_{i,l}$,

$\displaystyle\Delta_{i,l}(k)$ $\displaystyle=$
$\displaystyle(\varphi_{l}(kA_{0})-\frac{1}{l!}I)[\bar{F}_{n,i}^{u,2,simp}-\bar{F}_{n,i}^{simp}]+\frac{1}{l!}[\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}]$
$\displaystyle+k\varphi_{l+1}(kA_{0})A[\bar{F}_{n,i}^{u,2,A,simp}-\bar{F}_{n,i}^{u,2,simp}]+\frac{k}{(l+1)!}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,simp}]$
$\displaystyle+k^{2}\varphi_{l+2}(kA_{0})A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,A,simp}],$

from what, applying $A_{0}^{-1}$ and then (8),

$\displaystyle A_{0}^{-1}\Delta_{i,l}(k)$ $\displaystyle=$ $\displaystyle
k\varphi_{l+1}(kA_{0})[\bar{F}_{n,i}^{u,2,simp}-\bar{F}_{n,i}^{simp}]+\frac{1}{l!}A_{0}^{-1}[\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}]$
$\displaystyle+k^{2}\varphi_{l+2}(kA_{0})A[\bar{F}_{n,i}^{u,2,A,simp}-\bar{F}_{n,i}^{u,2,simp}]+\frac{k}{(l+1)!}A_{0}^{-1}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,simp}]$
$\displaystyle+k^{2}\varphi_{l+2}(kA_{0})A_{0}^{-1}A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,A,simp}].$

Now, the first and third terms are $O(k^{3})$ because of previously justified
bounds. The second term is $O(k^{3})$ using (33), (A7) and the fact that
$\|A_{0}^{-1}(\bar{K}_{n,i}-\bar{K}_{n,i}^{simp})\|=O(k^{3})$. The fourth term
is $O(k^{3})$ because, using (A8),
$\|A_{0}^{-1}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,simp}]\|=O(k^{2})$. The fifth
term is $O(k^{3})$ because $\|A(\bar{K}_{n,i}-u(t_{n}))\|=O(k)$ as it is
deduced from the proof of Theorem 1 and then, applying (A8),
$\|A_{0}^{-1}A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,2,A,simp}]\|=O(k)$. Summing
up, $\|A_{0}^{-1}\Delta_{i,l}(k)\|=O(k^{3})$, from what
$A_{0}^{-1}(\bar{u}_{n+1}-\bar{u}_{n+1}^{simp})$ and
$A_{0}^{-1}\rho_{n}^{simp}$ are $O(k^{4})$.

As for $p=3$, the proof follows the same lines of the case $p=2$ with the
difference that now,

$\displaystyle\bar{K}_{n,i}^{simp}=\phi_{0,2,u(t_{n}),u(t_{n})}(c_{i}k)+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\tilde{\phi}_{l,j}(c_{r}k),$
(73)

where

$\displaystyle\tilde{\phi}_{l,j}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A-\frac{l}{\tau})\tilde{\phi}_{l,j}(\tau)+\frac{1}{(l-1)!\tau}\bar{F}_{n,j}^{simp},$
$\displaystyle\tilde{\phi}_{l,j}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}\bar{F}_{n,j}^{simp},$
$\displaystyle\partial\tilde{\phi}_{l,j}(\tau)$ $\displaystyle=$
$\displaystyle\partial[\frac{1}{l!}\bar{F}_{n,j}^{K,3,simp}+\frac{\tau}{(l+1)!}A\bar{F}_{n,j}^{K,3,A,simp}],$

with $\bar{F}_{n,j}^{K,3,simp}$ and $\bar{F}_{n,j}^{K,3,A,simp}$ the
corresponding functions taken from Table 1:

$\displaystyle\bar{F}_{n,j}^{K,3,simp}$ $\displaystyle=$ $\displaystyle
f\big{(}t_{n}+c_{j}k,u(t_{n})+k[c_{j}Au(t_{n})+\sum_{jj=1}^{j-1}\sum_{l,r}\frac{\lambda_{j,jj,l,r}}{l!}f(t_{n},u(t_{n}))]\big{)},$
(74) $\displaystyle\bar{F}_{n,j}^{K,3,A,simp}$ $\displaystyle=$ $\displaystyle
f(t_{n},u(t_{n})).$

From this, considering the notation

$\displaystyle
T_{l,j}(\tau)=\frac{1}{l!}\bar{F}_{n,j}^{K,3,simp}+\frac{\tau}{(l+1)!}A\bar{F}_{n,j}^{K,3,A,simp},$

the following is satisfied

$\displaystyle(\tilde{\phi}_{l,j}-T_{l,j})^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A-\frac{l}{\tau})(\tilde{\phi}_{l,j}-T_{l,j})(\tau)+\frac{1}{(l-1)!\tau}[\tilde{F}_{n,j}^{simp}-\bar{F}_{n,j}^{K,3,simp}]$
$\displaystyle+\frac{1}{l!}A[\bar{F}_{n,j}^{K,3,simp}-\bar{F}_{n,j}^{K,3,A,simp}]+\frac{\tau}{(l+1)!}A^{2}\bar{F}_{n,j}^{K,3,A,simp},$
$\displaystyle(\tilde{\phi}_{l,j}-T_{l,j})(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\tilde{F}_{n,j}^{simp}-\tilde{F}_{n,j}^{K,3,simp}],$
$\displaystyle\partial(\tilde{\phi}_{l,j}-T_{l,j})(\tau)$ $\displaystyle=$
$\displaystyle 0.$

Solving this,

$\displaystyle(\tilde{\phi}_{l,j}-T_{l,j})(c_{r}k)$ $\displaystyle=$
$\displaystyle\varphi_{l}(c_{r}kA_{0})[\bar{F}_{n,j}^{simp}-\bar{F}_{n,j}^{K,3,simp}]+c_{r}k\varphi_{l+1}(c_{r}kA_{0})A[\bar{F}_{n,j}^{K,3,simp}-\bar{F}_{n,j}^{K,3,A,simp}]$
$\displaystyle+(c_{r}k)^{2}\varphi_{l+2}(\tau
A_{0})A^{2}\tilde{F}_{n,j}^{K,3,A,simp}$ $\displaystyle=$ $\displaystyle
O(k^{2}),$

where, for the last equality, an induction argument over the stages has been
used, (73), (74), as well as the assumed regularity and (A7). From this,

$\displaystyle\bar{K}_{n,i}^{simp}$ $\displaystyle=$ $\displaystyle
u(t_{n})+c_{i}kAu(t_{n})+\frac{(c_{i}k)^{2}}{2}A^{2}u(t_{n})$ (75)
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\big{[}\frac{1}{l!}\bar{F}_{n,j}^{K,3,simp}+\frac{c_{r}k}{(l+1)!}A\bar{F}_{n,j}^{K,3,A,simp}\big{]}+O(k^{3}).$

On the other hand,

$\displaystyle\bar{u}_{n+1}-\bar{u}_{n+1}^{simp}=k\sum_{i,l=1}^{s}\mu_{i,l}\tilde{\Delta}_{i,l}(k),$
(76)

where

$\displaystyle\tilde{\Delta}_{i,l}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle(A-\frac{l}{\tau})\Delta_{i,l}(\tau)+\frac{1}{(l-1)!\tau}[\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}],$
$\displaystyle\tilde{\Delta}_{i,l}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\bar{F}_{n,i}-\bar{F}_{n,i}^{simp}],$
$\displaystyle\partial\tilde{\Delta}_{i,l}(\tau)$ $\displaystyle=$
$\displaystyle\partial\bigg{[}\frac{1}{l!}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,simp}]+\frac{\tau}{(l+1)!}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A,simp}]+\frac{\tau^{2}}{(l+2)!}A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A^{2},simp}]\bigg{]},$

where $\bar{F}_{n,i}^{u,3,simp}$, $\bar{F}_{n,i}^{u,3,A,simp}$ and
$\bar{F}_{n,i}^{u,3,A^{2},simp}$ are again the corresponding expressions in
Table 2. Considering then the notation

$\displaystyle\tilde{W}_{i,l}(\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,simp}]+\frac{\tau}{(l+1)!}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A,simp}]$
(77)
$\displaystyle+\frac{\tau^{2}}{(l+2)!}A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A^{2},simp}],$

it follows that

$\displaystyle(\tilde{\Delta}_{i,l}-\tilde{W}_{i,l})^{\prime}(\tau)$
$\displaystyle=$
$\displaystyle(A-\frac{l}{\tau})(\tilde{\Delta}_{i,l}-\tilde{W}_{i,l})(\tau)+\frac{1}{l!}A[\bar{F}_{n,i}^{u,3,A,simp}-\bar{F}_{n,i}^{u,3,simp}]$
$\displaystyle+\frac{\tau}{(l+1)!}A^{2}[\bar{F}_{n,i}^{u,3,A^{2},simp}-\bar{F}_{n,i}^{u,3,A,simp}]+\frac{\tau^{2}}{(l+2)!}A^{3}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A^{2},simp}]$
$\displaystyle+\frac{1}{(l-1)!\tau}[\bar{F}_{n,i}^{u,3,simp}-\bar{F}_{n,i}^{simp}],$
$\displaystyle(\tilde{\Delta}_{i,l}-\tilde{W}_{i,l})(0)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}[\bar{F}_{n,i}^{u,3,simp}-\bar{F}_{n,i}^{simp}],$
$\displaystyle\partial(\tilde{\Delta}_{i,l}-\tilde{W}_{i,l})(\tau)$
$\displaystyle=$ $\displaystyle 0.$

From this,

$\displaystyle(\tilde{\Delta}_{i,l}-\tilde{W}_{i,l})(k)$ $\displaystyle=$
$\displaystyle\varphi_{l}(kA_{0})[\bar{F}_{n,i}^{u,3,simp}-\bar{F}_{n,i}^{simp}]+k\varphi_{l+1}(kA_{0})A[\bar{F}_{n,i}^{u,3,A,simp}-\bar{F}_{n,i}^{u,3,simp}]$
(78)
$\displaystyle+k^{2}\varphi_{l+2}(kA_{0})A^{2}[\bar{F}_{n,i}^{u,3,A^{2},simp}-\bar{F}_{n,i}^{u,3,A,simp}]$
$\displaystyle+k^{3}\varphi_{l+3}(kA_{0})A^{3}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A^{2},simp}].$

Considering (75), (A7) and the assumed regularity,
$\|\bar{F}_{n,i}^{u,3,simp}-\bar{F}_{n,i}^{simp}\|=O(k^{3})$. In a similar
way, the rest of the terms which turn up here, as well as those in (77) for
$\tau=k$ are $O(k^{3})$, so that $\tilde{\Delta}_{i,l}(k)=O(k^{3})$.
Therefore, using (76), $\bar{u}_{n+1}-\bar{u}_{n+1}^{simp}=O(k^{4})$, from
what, using Theorem 1 with $p=3$, $\rho_{n}^{simp}=O(k^{4})$ and, finally,
through the same proof of Theorem 2 with $p=3$, the first result of the
present theorem follows.

As for the second part of the theorem, we firstly notice that

$\displaystyle
A_{0}^{-1}(\bar{K}_{n,i}-\bar{K}_{n,i}^{simp})=k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}A_{0}^{-1}\Delta_{j,l}(c_{r}k),$
(79)

where $\Delta_{j,l}$ is the function in the proof of Theorem 3. From there,
$A_{0}^{-1}\Delta_{j,l}(c_{r}k)=O(k^{3})$, which implies that
$A_{0}^{-1}(\bar{K}_{n,i}-\bar{K}_{n,i}^{simp})=O(k^{4})$. Then,

$\displaystyle
A_{0}^{-1}(\bar{u}_{n+1}-\bar{u}_{n+1}^{simp})=k\sum_{i,l}\mu_{i,l}A_{0}^{-1}\tilde{\Delta}_{i,l}(k),$
(80)

and from (78) and (77),

$\displaystyle A_{0}^{-1}\tilde{\Delta}_{i,l}(k)$ $\displaystyle=$
$\displaystyle\frac{1}{l!}A_{0}^{-1}(\bar{F}_{n,i}-\bar{F}_{n,i}^{simp})+k\varphi_{l+1}(kA_{0})[\bar{F}_{n,i}^{u,3,simp}-\bar{F}_{n,i}^{simp}]$
$\displaystyle+k^{2}\varphi_{l+2}(kA_{0})A[\bar{F}_{n,i}^{u,3,A,simp}-\bar{F}_{n,i}^{u,3,simp}]+\frac{k}{(l+1)!}A_{0}^{-1}A[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,simp}]$
$\displaystyle+k^{3}\varphi_{l+3}(kA_{0})A^{2}[\bar{F}_{n,i}^{u,3,A^{2},simp}-\bar{F}_{n,i}^{u,3,A,simp}]+\frac{k^{2}}{(l+2)!}A_{0}^{-1}A^{2}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A,simp}]$
$\displaystyle+k^{3}\varphi_{l+3}(kA_{0})A_{0}^{-1}A^{3}[\bar{F}_{n,i}-\bar{F}_{n,i}^{u,3,A^{2},simp}].$

Now, the first term in the right hand side is $O(k^{4})$ using (33) and (79).
The second, third and fifth terms are also $O(k^{4})$ using the assumed
regularity and (A7). Finally, the fourth, sixth and seventh terms are
$O(k^{4})$ by using (A8) and the fact that
$\|\bar{K}_{n,i}-\bar{K}_{n,i}^{u,3,simp}\|_{h}=O(k^{3})$,
$\|A(\bar{K}_{n,i}-\bar{K}_{n,i}^{u,3,A,simp})\|_{h}=O(k^{2})$ and
$\|A^{2}(\bar{K}_{n,i}-\bar{K}_{n,i}^{u,3,A^{2},simp})\|_{h}=O(k)$,
considering the proof of Theorem 1. Therefore, from (80),
$A_{0}^{-1}(\bar{u}_{n+1}-\bar{u}_{n+1}^{simp})=O(k^{5})$, from what
$A_{0}^{-1}\rho_{n}^{simp}=O(k^{5})$, and the result follows considering the
same proof as that of the second part of Theorem 2 with $p=3$.

###  6.4 Proof of Theorem 4

In the following $\rho_{n}$ and $\rho_{n,h}$ in fact correspond to
$\rho_{n}^{sim}$ and $\rho_{n,h}^{sim}$. By definition,

$\displaystyle e_{n+1,h}=(U_{h}^{n+1}-\bar{U}_{h}^{n+1})+\rho_{n,h}.$ (81)

Then, for $p=1$, using Table 3,

$\displaystyle U_{h}^{n+1}-\bar{U}_{h}^{n+1}$ $\displaystyle=$ $\displaystyle
e^{kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-k\varphi_{1}(kA_{h,0})D_{h}O(e_{n,h})+k^{2}\varphi_{2}(kA_{h,0})C_{h}O(e_{n,h})$
$\displaystyle+k\sum_{i,l=1}^{s}\mu_{i,l}\big{[}\varphi_{l}(kA_{h,0})[F_{n,i,h}-\bar{F}_{n,i,h}]+k\varphi_{l+1}(kA_{h,0})C_{h}O(e_{n,h})\big{]},$

where the term $C_{h}O(e_{n,h})$ is just non-zero for Robin/Neumann boundary
conditions. On the other hand, using (31),

$\displaystyle K_{n,i,h}-\bar{K}_{n,i,h}$ $\displaystyle=$ $\displaystyle
e^{c_{i}kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\varphi_{l}(c_{r}kA_{h,0})[F_{n,j,h}-\bar{F}_{n,j,h}],$

and so, inductively, because of (H3), $F_{n,i,h}-\bar{F}_{n,i,h}=O(e_{n,h})$.
Inserting this in (81), using also (8), (H2b) and (H2c),

$e_{n+1,h}=e^{kA_{h,0}}e_{n,h}+O(ke_{n,h})+\rho_{n+1,h}.$

Then, applying Gronwall lemma, under the first set of hypotheses of Theorem 2,
using the bound for $\rho_{n,h}$ and also Theorem 3, the result for $p=1$
follows. Under the second set of hypotheses, the result for $p=1$ can be
obtained using a summation-by-parts argument, which can be applied in a
similar way to what was done for exponential splitting methods in [3] taking
(35), the regularity (36) and the bound for $A_{h,0}^{-1}\rho_{n,h}$ into
account.

For $p=2$, we will separate the cases corresponding to Dirichlet and to R/N
boundary conditions.

With Dirichlet boundary conditions, considering Table 4,

$\displaystyle
U_{h}^{n+1}-\bar{U}_{h}^{n+1}=e^{kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-k^{2}\varphi_{2}(kA_{h,0})D_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle+k^{3}\varphi_{3}(kA_{h,0})C_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle+k\sum_{i,l=1}^{s}\mu_{i,l}\bigg{[}\varphi_{l}(kA_{h,0})[F_{n,i,h}-\bar{F}_{n,i,h}]-k\varphi_{l+1}(kA_{h,0})D_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle\hskip
56.9055pt+k^{2}\varphi_{l+2}(kA_{h,0})C_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})\bigg{]},$

where $F_{n,i,h}-\bar{F}_{n,i,h}=O(e_{n,h})$ with the same argument of the
previous theorem. Then, using (H2b) and (H2c),

$e_{n+1,h}=e^{kA_{h,0}}e_{n,h}+O(k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h}+ke_{n,h})+\rho_{n+1,h},$

from what, using also (34) and the first set of hypotheses of Theorem 3,

$e_{n+1,h}=e^{kA_{h,0}}e_{n,h}+O(k^{2}\nu_{h}+ke_{n,h}+k^{3}+k\varepsilon_{h}),$

and the result follows by a discrete Gronwall lemma. Under the second set of
hypotheses, the summation-by-parts argument can be applied as in the previous
theorem using the respective results of Theorems 2 and 3.

On the other hand, with R/N boundary conditions and $p=2$, using again Table
4,

$\displaystyle U_{h}^{n+1}-\bar{U}_{h}^{n+1}$ $\displaystyle=$ $\displaystyle
e^{kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-k\varphi_{1}(kA_{h,0})D_{h}O(e_{n,h})$
$\displaystyle+k^{2}\varphi_{2}(kA_{h,0})[C_{h}O(e_{n,h})-D_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})]$
$\displaystyle+k^{3}\varphi_{3}(kA_{h,0})C_{h}O(\mu_{k,1}+\frac{e_{n,h}}{h^{\gamma}}+\nu_{h}+\frac{e_{n,h}}{k})$
$\displaystyle+k\sum_{i,l=1}^{s}\mu_{i,l}\bigg{[}\varphi_{l}(kA_{h,0})[F_{n,i,h}-\bar{F}_{n,i,h}]$
$\displaystyle\hskip
56.9055pt+k\varphi_{l+1}(kA_{h,0})[C_{h}O(e_{n,h}+k\mu_{k,1})-D_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})]$
$\displaystyle\hskip
56.9055pt+k^{2}\varphi_{l+2}(kA_{h,0})C_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})\bigg{]},$

where

$\displaystyle K_{n,i,h}-\bar{K}_{n,i,h}$ $\displaystyle=$ $\displaystyle
e^{c_{i}kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))+(c_{i}k)^{2}\varphi_{2}(c_{i}kA_{h,0})C_{h}O(e_{n,h})$
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\big{[}\varphi_{l}(c_{r}kA_{h,0})[F_{n,j,h}-\bar{F}_{n,j,h}]+c_{r}k\varphi_{l+1}(c_{r}kA_{h,0})C_{h}O(e_{n,h})\big{]},$

and so again, inductively, because of (H2b) and (H3),
$F_{n,i,h}-\bar{F}_{n,i,h}=O(e_{n,h})$. Then,

$e_{n+1,h}=e^{kA_{h,0}}e_{n,h}+O(ke_{n,h}+k^{2}\mu_{k,1}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h})+\rho_{n+1,h},$

and arguing as before, the result follows.

For $p=3$, the theorem follows as before but taking now into account that,
with Dirichlet boundary conditions, using Table 5,

$\displaystyle U_{h}^{n+1}-\bar{U}_{h}^{n+1}$ $\displaystyle=$ $\displaystyle
e^{kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-k^{2}\varphi_{2}(kA_{h,0})D_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle+k^{3}\varphi_{3}(kA_{h,0})[C_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})-D_{h}O(\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}}+\frac{\mu_{k,1}}{h^{\gamma}})]$
$\displaystyle+k^{4}\varphi_{4}(kA_{h,0})C_{h}O(\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}}+\frac{\mu_{k,1}}{h^{\gamma}})$
$\displaystyle+k\sum_{i,l=1}^{s}\mu_{i,l}\bigg{[}\varphi_{l}(kA_{h,0})[F_{n,i,h}-\bar{F}_{n,i,h}]$
$\displaystyle\hskip
56.9055pt+k\varphi_{l+1}(kA_{h,0})[C_{h}O(k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h})-D_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}}+\frac{k\mu_{k,1}}{h^{\gamma}})]$
$\displaystyle\hskip
56.9055pt+k^{2}\varphi_{l+2}(kA_{h,0})[C_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}}+\frac{k\mu_{k,1}}{h^{\gamma}})-D_{h}O(\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}}+\frac{\mu_{k,1}}{h^{\gamma}})]$
$\displaystyle\hskip
56.9055pt+k^{3}\varphi_{l+3}(kA_{h,0})C_{h}O(\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}}+\frac{\mu_{k,1}}{h^{\gamma}})\bigg{]},$

where

$\displaystyle K_{n,i,h}-\bar{K}_{n,i,h}$ $\displaystyle=$ $\displaystyle
e^{c_{i}kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-(c_{i}k)^{2}\varphi_{2}(kA_{h,0})D_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle+(c_{i}k)^{3}\varphi_{3}(c_{i}kA_{h,0})C_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\big{[}\varphi_{l}(c_{r}kA_{h,0})[F_{n,j,h}-\bar{F}_{n,j,h}]-c_{r}k\varphi_{l+1}(c_{r}kA_{h,0})D_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle\hskip
85.35826pt+(c_{r}k)^{2}\varphi_{l+2}(c_{r}kA_{h,0})C_{h}O(\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})\big{]},$

and therefore
$F_{n,i,h}-\bar{F}_{n,i,h}=O(e_{n,h}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h})$.
Using this,

$e_{n+1,h}=e^{kA_{h,0}}e_{n,h}+O(ke_{n,h}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h}+k^{3}\frac{\mu_{k,1}}{h^{\gamma}})+\rho_{n+1,h},$

from what arguing as in previous theorems, but using now the results in
Theorem 3, the result follows under the first and second set of hypotheses.

With R/N boundary conditions, using again Table 5,

$\displaystyle U_{h}^{n+1}-\bar{U}_{h}^{n+1}$ $\displaystyle=$ $\displaystyle
e^{kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-k\varphi_{1}(kA_{h,0})D_{h}O(e_{n,h})$
$\displaystyle+k^{2}\varphi_{2}(kA_{h,0})[C_{h}O(e_{n,h})-D_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})]$
$\displaystyle+k^{3}\varphi_{3}(kA_{h,0})[C_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})-D_{h}O(\mu_{k,1}+\mu_{k,2}+\frac{e_{n,h}}{k^{2}}+\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}})]$
$\displaystyle+k^{4}\varphi_{4}(kA_{h,0})C_{h}O(\mu_{k,1}+\mu_{k,2}+\frac{e_{n,h}}{k^{2}}+\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}})$
$\displaystyle+k\sum_{i,l=1}^{s}\mu_{i,l}\big{[}\varphi_{l}(kA_{h,0})[F_{n,i,h}-\bar{F}_{n,i,h}]$
$\displaystyle\hskip
28.45274pt+k\varphi_{l+1}(kA_{h,0})[C_{h}O(e_{n,h}+k\mu_{k,1}+k^{2}\mu_{k,2}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h})-D_{h}O(\mu_{k,1}+k\mu_{k,2}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})]$
$\displaystyle\hskip
28.45274pt+k^{2}\varphi_{l+2}(kA_{h,0})[C_{h}O(\mu_{k,1}+k\mu_{k,2}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})-D_{h}O(\mu_{k,1}+\mu_{k,2}+\frac{e_{n,h}}{k^{2}}+\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}})]$
$\displaystyle\hskip
28.45274pt+k^{3}\varphi_{l+3}(kA_{h,0})C_{h}O(\mu_{k,1}+\mu_{k,2}+\frac{e_{n,h}}{k^{2}}+\nu_{h}+\frac{e_{n,h}}{kh^{\gamma}})\big{]},$

where

$\displaystyle K_{n,i,h}-\bar{K}_{n,i,h}$ $\displaystyle=$ $\displaystyle
e^{c_{i}kA_{h,0}}(U_{h}^{n}-P_{h}u(t_{n}))-c_{i}k\varphi_{1}(c_{i}kA_{h,0})D_{h}O(e_{n,h})$
$\displaystyle+(c_{i}k)^{2}\varphi_{2}(c_{i}kA_{h,0})[C_{h}O(e_{n,h})-D_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})]$
$\displaystyle+(c_{i}k)^{3}\varphi_{3}(c_{i}kA_{h,0})[C_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})$
$\displaystyle+k\sum_{j=1}^{i-1}\sum_{l,r=1}^{s}\lambda_{i,j,l,r}\bigg{[}\varphi_{l}(c_{r}kA_{h,0})[F_{n,j,h}-\bar{F}_{n,j,h}]$
$\displaystyle\hskip
28.45274pt+c_{r}k\varphi_{l+1}(c_{r}kA_{h,0})[C_{h}O(e_{n,h}+k\mu_{k,1})-D_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})]$
$\displaystyle\hskip
28.45274pt+(c_{r}k)^{2}\varphi_{l+2}(c_{r}kA_{h,0})C_{h}O(\mu_{k,1}+\frac{e_{n,h}}{k}+\nu_{h}+\frac{e_{n,h}}{h^{\gamma}})\bigg{]},$

and therefore
$F_{n,i,h}-\bar{F}_{n,i,h}=O(e_{n,h}+k^{2}\mu_{k,1}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h})$.
Using this,

$e_{n+1,h}=e^{kA_{h,0}}e_{n,h}+O(ke_{n,h}+k^{2}\mu_{k,1}+k^{2}\nu_{h}+\frac{k^{2}}{h^{\gamma}}e_{n,h}+k^{3}\mu_{k,2})+\rho_{n+1,h},$

from what the result follows as before.