# A fourth-order exponential time differencing scheme with dimensional
splitting for non-linear reaction-diffusion systems

E.O. Asante-Asamani<EMAIL_ADDRESS>Department of Mathematics, Clarkson
University
Potsdam NY 13676, USA  A. Kleefeld<EMAIL_ADDRESS>Jülich
Supercomputing Centre, Forschungszentrum Jülich GmbH,
52425 Jülich, Germany  Faculty of Medical Engineering and Technomathematics,
University of Applied Sciences Aachen,
52428 Jülich, Germany  B.A. Wade<EMAIL_ADDRESS>Department of
Mathematics, University of Louisiana at Lafayette,
Lafayette LA 70504, USA

###### Abstract

A fourth-order exponential time differencing (ETD) Runge-Kutta scheme with
dimensional splitting is developed to solve multidimensional non-linear
systems of reaction-diffusion equations (RDE). By approximating the matrix
exponential in the scheme with the A-acceptable Padé (2,2) rational function,
the resulting scheme (ETDRK4P22-IF) is verified empirically to be fourth-order
accurate for several RDE. The scheme is shown to be more efficient than
competing fourth-order ETD and IMEX schemes, achieving up to 20X speed-up in
CPU time. Inclusion of up to three pre-smoothing steps of a lower order
L-stable scheme facilitates efficient damping of spurious oscillations arising
from problems with non-smooth initial/boundary conditions.

###### keywords:

Exponential time differencing , semilinear parabolic problems , reaction-
diffusion equations , fourth-order time stepping.

##  1 Introduction

Time-dependent reaction-diffusion equations (RDE) are mathematical models that
describe the spatio-temporal dynamics of many physical processes in diverse
fields such as cell and developmental biology [1, 2, 3], physics [4],
geological engineering [5], and finance [6]. Mathematically, RDE are described
by the system of partial differential equations (PDE) of the form

$\begin{split}\frac{\partial u}{\partial t}&=D\Delta u+f(u,t),\quad\text{ in
}\Omega_{T}\\\ u(x,t)&=g(x,t),\text{ on }\Gamma_{T}.\end{split}$ (1)

Here, $\Omega$ is a bounded open subset of $\mathbb{R}^{d},d\in\mathbb{N}$,
$\Omega_{T}:=\Omega\times(0,T]$ with parabolic boundary
${\Gamma_{T}:=\overline{\Omega_{T}}\setminus\Omega_{T}}$. The solution
$u:\overline{\Omega_{T}}\rightarrow\mathbb{R}^{s}$, describes the
concentration of $s$ species with diffusion coefficients
$D=\mathrm{diag}(d_{1},d_{2},\ldots,d_{s}),\ d_{i}>0\ \forall i$. The function
$f:\Omega_{T}\rightarrow\mathbb{R}^{s}$ describes the interaction between all
the species and is assumed to be sufficiently smooth.

In practice, the spatial derivatives in (1) are first discretized on a
partition of $\Omega$ with $m^{d}$ points to obtain a system of ordinary
differential equations (ODE)

$\displaystyle\frac{dU}{dt}+AU$ $\displaystyle=F(U(t),t),\quad
U\in\mathbb{R}^{s\times m^{d}}$ (2) $\displaystyle U(x,0)$
$\displaystyle=U_{0}(x)$

where $U(t)=(u_{1}(t),u_{2}(t),\ldots,u_{s}(t),\ldots,u(t)_{s\times
m^{d}})^{T}$, $A$ is an $(s\times m^{d})\times(s\times m^{d})$ matrix
approximation of $-D\Delta$ and $F\in\mathbb{R}^{s\times m^{d}}$ is typically
a non-linear map approximating $f$ on the numerical grid for all $s$ species.
The major computational challenges in solving (2) include numerical stiffness
due to the spatial diffusion, coupling of very fast and very slow reaction
kinetics, simulating problems in high spatial dimensions ($d\geq 2$) with non-
linear source terms, spurious oscillations arising from non-smooth initial
conditions or mismatched initial and boundary conditions as well as
maintaining positive solutions especially for biological applications [7].

Among the host of time stepping methods available to solve stiff ODE systems
[7, 8, 9, 10, 11] are the class of exponential time differencing (ETD) Runge-
Kutta schemes. These schemes approximate the exact solution of the system (2),

$U(t)=\mathrm{e}^{-tA}U(0)+\int_{0}^{t}\mathrm{e}^{-(t-s)A}F(U(s),s)\;\mathrm{d}s\,,\quad
t\in[0,T]$ (3)

by approximating the integral and matrix exponential. In 2002, Cox and
Matthews [12] introduced a class of ETD Runge-Kutta (ETDRK) schemes which
approximate the integral in (3) by interpolating the non-linear function $F$
with polynomials. The major attraction of these schemes is in the accurate
treatment of the linear (diffusion) term via matrix exponential leading to
stable evolution of stiff problems.

The computational time for ETDRK schemes has been improved by employing
rational approximations of the matrix exponential. Several second-order and
fourth-order ETDRK schemes have been developed that approximate the matrix
exponential using Padé rational functions (ETDRK-Padé schemes) [13, 14, 15,
16] as well as a rational function with real and distinct poles (RDP) leading
to the ETDRDP scheme in [17]. To further improve the computational efficiency
of second-order ETDRK schemes when solving multidimensional RDE, we introduced
a class of dimensionally split second-order ETDRK schemes to take advantage of
fast solvers for tridiagonal linear systems. Utilizing Padé (1,1) and RDP
rational functions we developed efficient second-order ETDRK schemes with
dimensional splitting in [18] and [19], respectively. The dimensional
splitting scheme developed with RDP rational functions is referred to as
ETDRDP-IF. In spite of this progress, the savings in CPU time obtained from
dimensional splitting has, until now, been limited to second-order ETDRK
schemes.

In this work, we introduce a new class of fourth-order ETDRK schemes with
dimensional splitting. By approximating the matrix exponential with the
A-acceptable Padé (2,2) rational function, we create an efficient fourth-order
scheme we call ETDRK4P22-IF. Empirical convergence analysis validates the
fourth-order accuracy of the scheme across various prototypical scalar and
systems of RDE with Dirichlet and Neumann boundary conditions posed in
$\mathbb{R}^{2}$. The new scheme is shown to significantly improve the run
time of the existing fourth order ETDRK4P22 scheme (without dimensional
splitting) developed by Yousuf [16, 15], in some cases achieving 20X speed-up
in CPU time with serial implementation of the algorithm. The one-step nature
of the scheme is also shown to provide a substantial computational advantage
over competing fourth-order multistep IMEX schemes. By presmoothing
ETDRK4P22-IF with a few steps of a lower order L-stable ETDRK scheme, we show
that ETDRK4P22-IF can accurately solve RDE with non-smooth initial or boundary
conditions.

The remainder of the manuscript is organized as follows. In the Section 2 we
describe the finite difference scheme used to obtain a fourth-order
approximation the Laplacian. This is followed by details of the dimensional
splitting of the fourth-order ETDRK scheme along with a rational approximation
of the matrix exponential using Padé(2,2) functions to yield the ETDRK4P22-IF
scheme. Serial and parallel algorithms for implementing the ETDRK4P22-IF are
also presented in this section. The convergence and efficiency of ETDRK4P22-IF
is investigated in Section 3 for problems with and without exact solution.
Finally, we wrap up with concluding remarks in Section 4.

##  2 Materials and Methods

###  2.1 Discretization in space

To discretize the Laplacian $\Delta=\partial_{xx}+\partial_{yy}$ on the domain
$[a,b]^{2}$, we partition the spatial domain in each direction into $m+2$
points with a uniform mesh of size $h=\frac{b-a}{m+1}$ with $m\geq 3.$ If we
set $x_{j}=a+jh$ with $j=0,1,2,\ldots,m+1$, then the second-order partial
derivative of a function $w(x,t)$ with respect to $x$ at $x_{j}$ can be
discretized using the fourth-order central difference scheme ([20], p. 339)

$\partial_{xx}w|_{x_{j}}=\frac{1}{12h^{2}}(-W_{j-2}+16W_{j-1}-30W_{j}+16W_{j+1}-W_{j+2})+\mathcal{O}(h^{4})\,$
(4)

$\quad j=2,3,\ldots,m-1$, with $W_{j}\approx w(x_{j},t)$. The approximation
used for the nodes at $j=0,1,m,m+1$ depends on the boundary condition. Since
the solution at $x_{0}$ and $x_{m+1}$ are known for homogeneous Dirichlet
boundary conditions, the derivatives at $x_{1}$ and $x_{m}$ are estimated with
the formulas

$\begin{split}\partial_{xx}w|_{x_{1}}&\approx\frac{1}{12h^{2}}(11W_{0}-20W_{1}+6W_{2}+4W_{3}-W_{4})\,,\\\
\partial_{xx}w|_{x_{m}}&\approx\frac{1}{12h^{2}}(-W_{m-3}+4W_{m-2}+6W_{m-1}-20W_{m}+11W_{m+1})\,\end{split}$
(5)

obtained by extrapolating the values, $W_{-1},W_{m+2}$ using a fourth degree
Lagrange interpolating polynomial centered at $x_{0}$ and $x_{m+1}$.

For homogeneous Neumann boundary conditions, we use the approximation

$\begin{split}\partial_{xx}w|_{x_{1}}&\approx\frac{1}{12h^{2}}(16W_{0}-31W_{1}+16W_{2}-W_{3})\
,\\\
\partial_{xx}w|_{x_{m}}&\approx\frac{1}{12h^{2}}(-W_{m-2}+16W_{m-1}-30W_{m}+16W_{m+1}-W_{m+2})\,\end{split}$
(6)

which is derived by setting $W_{-1}=W_{1}$ and $W_{m}=W_{m+2}$. The
derivatives at $x_{0}$ and $x_{m+1}$ are set to

$\begin{split}\partial_{xx}w|_{x_{0}}&\approx\frac{1}{12h^{2}}(-30W_{0}+32W_{1}-2W_{2})\,,\\\
\partial_{xx}w|_{x_{m+1}}&\approx\frac{1}{12h^{2}}(-2W_{m-1}+32W_{m}-30W_{m+1})\,.\end{split}$
(7)

Here the Neumann boundary condition is discretized using a fourth-order
central difference approximation to the first-order partial derivative (see
[20] p. 339) with the additional assumption that $W_{-1}=W_{1}$ and
$W_{m+2}=W_{m}$ . This discretization is shown empirically to be fourth-order
accurate in our numerical experiments. A similar scheme has been used in [21]
to achieve fourth-order accuracy in irregular domains.

The matrix approximating the second-order partial derivatives in the PDE in
the presence of homogeneous Dirichlet boundary conditions has a bandwidth of
six whereas the matrix for the case of Neumann boundary conditions has a
bandwidth of five. To formulate the two dimensional Laplacian, we take
advantage of the Kronecker product. Let $A=A_{1}+A_{2}$ be the matrix
approximation of the two dimensional (2D) Laplacian with $A_{1}$ and $A_{2}$
as the matrix approximations of $\partial_{xx}$ and $\partial_{yy}$
respectively in 2D. Let $B_{p}$ be the matrix representations of the second-
order partial derivative in 1D with entries defined by equations (4)–(7) and
$I_{p}$ an $p$-dimensional identity matrix. Then $A_{1}=B_{p}\otimes I_{p}$
and $A_{2}=I_{p}\otimes B_{p}$ (see [7] for details). For Dirichlet boundary
conditions and $s$ species, $p=s\cdotp m$, whereas $p=s\cdotp(m+2)$ for
Neumann boundary. We proved in [19, 18] that $A_{1}$ and $A_{2}$ commute. This
is an important property which facilitates the formulation of the dimensional
splitting scheme presented in the next section.

###  2.2 Discretization in time

####  2.2.1 Dimensional splitting of fourth-order ETD Runge-Kutta scheme

In this section, we introduce a dimensional splitting version of the class of
fourth-order ETD Runge-Kutta schemes developed by Cox-Matthews [12] and
Kassam-Trefethen [22] for solving the semi-linear system of ordinary
differential equations in (2). Our motivation is to speed up the simulation
time when solving problems in high spatial dimensions by taking advantage of a
nicer matrix structure resulting from splitting dimensions.

The scheme is

$\begin{split}U_{n+1}&=\mathrm{e}^{-kA}U_{n}+\frac{1}{k^{2}}(-A)^{-3}[-4I+kA+\mathrm{e}^{-kA}(4I+3kA+k^{2}A^{2})]F(U_{n},t_{n})\\\
&+\frac{2}{k^{2}}(-A)^{-3}[2I-kA-\mathrm{e}^{-kA}(2+kA)]\left(F(a_{n},t_{n}+\frac{k}{2})+F(b_{n},t_{n}+\frac{k}{2})\right)\\\
&+\frac{1}{k^{2}}(-A)^{-3}[-4I+3kA-k^{2}A^{2}+\mathrm{e}^{-kA}(4I+kA)]F(c_{n},t_{n}+k)\,,\end{split}$
(8)

where

$\displaystyle a_{n}$
$\displaystyle=\mathrm{e}^{-\frac{k}{2}A}U_{n}-A^{-1}(\mathrm{e}^{-\frac{k}{2}A}-I)F(U_{n},t_{n})$
$\displaystyle b_{n}$
$\displaystyle=\mathrm{e}^{-\frac{k}{2}A}U_{n}-A^{-1}(\mathrm{e}^{-\frac{k}{2}A}-I)F(a_{n},t_{n}+\frac{k}{2})$
$\displaystyle c_{n}$
$\displaystyle=\mathrm{e}^{-\frac{k}{2}A}a_{n}-A^{-1}(\mathrm{e}^{-\frac{k}{2}A}-I)\left[2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n})\right]\,.$

For convenience, we will utilize the following more compact form of the
scheme,

$\begin{split}a_{n}&=\mathrm{e}^{-\frac{k}{2}A}U_{n}+\tilde{P}(kA)F(U_{n},t_{n})\\\
b_{n}&=\mathrm{e}^{-\frac{k}{2}A}U_{n}+\tilde{P}(kA)F(a_{n},t_{n}+\frac{k}{2})\\\
c_{n}&=\mathrm{e}^{-\frac{k}{2}A}a_{n}\tilde{P}(kA)[2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n})]\\\
U_{n+1}&=\mathrm{e}^{-kA}U_{n}+P_{1}(kA)F(U_{n},t_{n})+2P_{2}(kA)\left(F(a_{n},t_{n}+\frac{k}{2})+F(b_{n},t_{n}+\frac{k}{2})\right)\\\
&+P_{3}(kA)F(c_{n},t_{n}+k)\,,\end{split}$ (9)

with

$\displaystyle P_{1}(kA)$
$\displaystyle=\frac{1}{k^{2}}(-A)^{-3}[-4I+kA+\mathrm{e}^{-kA}(4I+3kA+k^{2}A^{2})]$
$\displaystyle P_{2}(kA)$
$\displaystyle=\frac{1}{k^{2}}(-A)^{-3}[2I-kA-\mathrm{e}^{-kA}(2I+kA)]$ (10)
$\displaystyle P_{3}(kA)$
$\displaystyle=\frac{1}{k^{2}}(-A)^{-3}[-4I+3kA-k^{2}A^{2}+\mathrm{e}^{-kA}(4I+kA)]$
$\displaystyle\tilde{P}(kA)$
$\displaystyle=-A^{-1}(\mathrm{e}^{-\frac{k}{2}A}-I).$

Consider the decomposition $A=A_{1}+A_{2}$, where $A_{1}$ and $A_{2}$ commute,
we can rewrite (2) in terms of $A_{1}$ and $A_{2}$ as follows. Consider a new
variable $V(t)=\mathrm{e}^{Bt}U(t)$ where $B$ is also a matrix, then
$V^{{}^{\prime}}(t)=\mathrm{e}^{Bt}U^{{}^{\prime}}(t)+B\mathrm{e}^{Bt}U(t)$.
Replacing $U^{{}^{\prime}}(t)$ with the expression $-AU(t)+F(U(t),t)$ from (2)
we have
$V^{{}^{\prime}}(t)=\mathrm{e}^{Bt}F(U(t),t)-\mathrm{e}^{Bt}AU(t)+B\mathrm{e}^{Bt}U(t)$.
Now, suppose $B$ commutes with $A$ then
$\mathrm{e}^{Bt}AU(t)=A\mathrm{e}^{Bt}U(t)$ and
$V^{{}^{\prime}}(t)=\mathrm{e}^{Bt}F(U(t),t)+(B-A)\mathrm{e}^{Bt}U(t)$.
Setting $B=A_{1}$ and using the decomposition of $A=A_{1}+A_{2}$ we obtain the
following semi-linear problem in the transformed variable $V$

$\displaystyle V^{{}^{\prime}}(t)+A_{2}V$ $\displaystyle=\tilde{F}(V(t),t)$
(11) $\displaystyle V(0)$ $\displaystyle=U(0)$

with $\tilde{F}(V,t)=\mathrm{e}^{A_{1}t}F(\mathrm{e}^{-A_{1}t}V,t)$. The
solution variable $U(t)$ can then be obtained from the transformed variable
$V(t)$ via the relation $U(t)=\mathrm{e}^{-A_{1}t}V(t)$. This reformulated
semi-linear problem is the basis for deriving the dimensional splitting
version of the class of ETDRK4 schemes. We accomplish this by applying the
existing ETDRK4 scheme in (8) to our reformulated problem (11) and the formula
$U(t)=\mathrm{e}^{-A_{1}t}V(t)$ to obtain a numerical solution for $U(t)$. We
used a similar technique in [19, 18] to develop the class of second-order
dimensionally split ETD Runge-Kutta schemes.

Applying the ETDRK4 scheme to (11) we obtain

$\displaystyle a_{n}^{*}$
$\displaystyle=\mathrm{e}^{-\frac{k}{2}A_{2}}V_{n}+\tilde{P}(kA_{2})\tilde{F}(V_{n},t_{n})$
$\displaystyle b_{n}^{*}$
$\displaystyle=\mathrm{e}^{-\frac{k}{2}A_{2}}V_{n}+\tilde{P}(kA_{2})\tilde{F}(a_{n}^{*},t_{n}+\frac{k}{2})$
$\displaystyle c_{n}^{*}$
$\displaystyle=\mathrm{e}^{-\frac{k}{2}A_{2}}a_{n}^{*}+\tilde{P}(kA_{2})\left[2\tilde{F}(b_{n}^{*},t_{n}+\frac{k}{2})-\tilde{F}(V_{n},t_{n})\right]\,$
$\displaystyle V_{n+1}$
$\displaystyle=\mathrm{e}^{-kA_{2}}V_{n}+P_{1}(kA_{2})\tilde{F}(V_{n},t_{n})+2P_{2}(kA_{2})\left(\tilde{F}(a_{n}^{*},t_{n}+\frac{k}{2})+\tilde{F}(b_{n}^{*},t_{n}+\frac{k}{2})\right)$
(12) $\displaystyle+P_{3}(kA_{2})\tilde{F}(c_{n}^{*},t_{n}+k).$

Now, using the fact that $V_{n}=\mathrm{e}^{A_{1}t_{n}}U_{n}$ and
$\tilde{F}(V_{n},t_{n})=\mathrm{e}^{A_{1}t_{n}}F(\mathrm{e}^{-A_{1}t_{n}}V_{n},t_{n})$
and setting ${a_{n}=\mathrm{e}^{-A_{1}t_{n+\frac{1}{2}}}a_{n}^{*}}$,
$b_{n}=\mathrm{e}^{-A_{1}t_{n+\frac{1}{2}}}b_{n}^{*}$,
$c_{n}=\mathrm{e}^{-A_{1}t_{n+1}}c_{n}^{*}$, we can rewrite (2.2.1) in terms
of $U_{n}$,

$\begin{split}a_{n}&=\mathrm{e}^{-\frac{k}{2}A_{2}}\mathrm{e}^{-\frac{k}{2}A_{1}}U_{n}+\tilde{P}(kA_{2})\mathrm{e}^{-\frac{k}{2}A_{1}}F(U_{n},t_{n})\\\
b_{n}&=\mathrm{e}^{-\frac{k}{2}A_{2}}\mathrm{e}^{-\frac{k}{2}A_{1}}U_{n}+\tilde{P}(kA_{2})F(a_{n},t_{n}+\frac{k}{2})\\\
c_{n}&=\mathrm{e}^{-\frac{k}{2}A_{2}}\mathrm{e}^{-\frac{k}{2}A_{1}}a_{n}+\tilde{P}(kA_{2})\left[2\mathrm{e}^{-\frac{k}{2}A_{1}}F(b_{n},t_{n}+\frac{k}{2})-\mathrm{e}^{-kA_{1}}F(U_{n},t_{n})\right]\,\\\
U_{n+1}&=\mathrm{e}^{-kA_{1}}\mathrm{e}^{-kA_{2}}U_{n}+P_{1}(kA_{2})\mathrm{e}^{-kA_{1}}F(U_{n},t_{n})+2P_{2}(kA_{2})\mathrm{e}^{-\frac{k}{2}A_{1}}G(a_{n},b_{n},t_{n}+\frac{k}{2})\\\
&+P_{3}(kA_{2})F(c_{n},t_{n}+k)\end{split}$ (13)

with

$G(a_{n},b_{n},t_{n}+\frac{k}{2})=F(a_{n},t_{n}+\frac{k}{2})+F(b_{n},t_{n}+\frac{k}{2}).$

Equation (13) is the dimensional splitting formulation of the fourth-order ETD
Runge-Kutta scheme first proposed by Cox and Matthews ([12], p. 436) for
diagonal systems and extended to non-diagonal systems by Kassam and Trefethen
([22], p. 1218).

####  2.2.2 Padé approximation of the matrix exponential

Whereas the scheme in (13) is fully discretized and can be applied to solve
the semi-linear problem in (2), it involves high powers of matrix inverses
which can be expensive to compute, at best, or fail to exist, at worst,
particularly in the cases of a coefficient matrix obtained from a PDE with
Neumann boundary conditions. The scheme also requires an efficient
implementation of the matrix exponential for reasonable computational time. To
resolve this, rational functions are commonly used to approximate the matrix
exponential and simplify the scheme in (13). Padé rational functions have been
used successfully to develop the class of ETD-Padé schemes leading to
efficient second-order [13, 14] and fourth-order ETD schemes [16, 15].
Rational functions with real and distinct poles have also been used to develop
the second-order schemes [19, 17]. In the discussion that follows we will
employ the Padé (2,2) rational function to approximate the matrix exponential
in (13) and simplify the scheme. This same rational function was used to
develop the fourth-order ETD-Padé scheme in [16] without dimensional
splitting, referred to here as ETDRK4P22. We denote the Padé (2,2)
approximation of $\mathrm{e}^{-kA}$ and $\mathrm{e}^{-\frac{k}{2}A}$ by
$R_{2,2}(kA)$ and $\tilde{R}_{2,2}(kA)$, respectively, thus obtaining the
representations

$\displaystyle\mathrm{e}^{-kA}$ $\displaystyle\approx
R_{2,2}(kA)=(12I-6kA+k^{2}A^{2})(12I+6kA+k^{2}A^{2})^{-1}\,,$ (14)
$\displaystyle\mathrm{e}^{-\frac{k}{2}A}$
$\displaystyle\approx\tilde{R}_{2,2}(kA)=(48I-12kA+k^{2}A^{2})(48I+12kA+k^{2}A^{2})^{-1}\,.$
(15)

Replacing the exponential functions in (13) by the corresponding rational
approximations above, the functions $P_{1},P_{2},P_{3},\tilde{P}$ can be
expressed as the rational functions

$\displaystyle P_{1}(kA)$
$\displaystyle=\frac{1}{k^{2}}(-A)^{-3}[-4I+kA+R_{2,2}(kA_{1})(4I+3kA+k^{2}A^{2})]$
$\displaystyle=k(2I-kA_{2})(12I+6kA_{2}+k^{2}A^{2})^{-1}$ $\displaystyle
P_{2}(kA)$ $\displaystyle=\frac{1}{k^{2}}(-A)^{-3}[2I-kA-
R_{2,2}(kA_{1})(2I+kA)]$ $\displaystyle=2k(12I+6kA_{2}+k^{2}A_{2}^{2})^{-1}$
$\displaystyle P_{3}(kA)$
$\displaystyle=\frac{1}{k^{2}}(-A)^{-3}[-4I+3kA-k^{2}A^{2}+R_{2,2}(kA_{1})(4I+kA)]$
$\displaystyle=k(2I+kA_{2})(12I+6kA_{2}+k^{2}A_{2}^{2})^{-1}$
$\displaystyle\tilde{P}(kA_{2})$
$\displaystyle=-A_{2}^{-1}(\tilde{R}_{2,2}(kA_{2})-I)=24k(48I+12kA+k^{2}A^{2})^{-1}$

leading to the simplified scheme,

$\displaystyle a_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA_{2})\tilde{R}_{2,2}(kA_{1})U_{n}+\tilde{P}(kA_{2})\tilde{R}_{2,2}(kA_{1})F(U_{n},t_{n})$
(16) $\displaystyle b_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA_{2})\tilde{R}_{2,2}(kA_{1})U_{n}+\tilde{P}(kA_{2})F(a_{n},t_{n}+\frac{k}{2})$
(17) $\displaystyle c_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA_{2})\tilde{R}_{2,2}(kA_{1})a_{n}$
$\displaystyle+\tilde{P}(kA_{2})[2\tilde{R}_{2,2}(kA_{1})F(b_{n},t_{n}+\frac{k}{2})-R_{2,2}(kA_{1})F(U_{n},t_{n})]$
(18) $\displaystyle U_{n+1}$
$\displaystyle=R_{2,2}(kA_{1})R_{2,2}(kA_{2})U_{n}+P_{1}(kA_{2})R_{2,2}(kA_{1})F(U_{n},t_{n})$
(19)
$\displaystyle+2P_{2}(kA_{2})\tilde{R}_{2,2}(kA_{1})G(a_{n},b_{n},t_{n}+\frac{k}{2})+P_{3}(kA_{2})F(c_{n},t_{n}+k)$

We refer to the simplified scheme in (16)–(19) as the ETDRK4P22-IF scheme.

The use of Padé (2,2) rational functions, which are well known to be
A-acceptable, suggests that our scheme will have limited damping properties
for non-smooth problems. For such problems, a few pre-smoothing steps
involving a lower order scheme with an L-acceptable rational approximation to
the exponential was recommended in [16]. We investigate presmoothing the
ETDRK4P22-IF scheme with the third order ETDRK scheme proposed in [16], which
approximates the matrix exponentials in (9)–(2.2.1) with Padé (0,3) rational
functions. Note that this smoothing scheme is not implemented with dimensional
splitting.

####  2.2.3 Partial fraction decomposition of rational functions

In implementing the ETDRK4P22-IF scheme, we utilize the partial fraction
decomposition of the rational functions. This facilitates easy parallelization
on multiple processors. The respective decomposition, whose complex poles and
weights are similar to ones used in [14], are as follows:

$R_{2,2}(z)=1+2\mathrm{Re}\left(\frac{w_{11}}{z-c_{1}}\right),\
c_{1}=-3+\mathrm{i}\sqrt{3},\
w_{11}=-\left(6+\mathrm{i}\frac{18}{\sqrt{3}}\right)\,.$ (20)
$\tilde{R}_{2,2}=1+4\mathrm{Re}\left(\frac{w_{11}}{z-c_{2}}\right),\
c_{2}=-6+2\mathrm{i}\sqrt{3},\
w_{11}=-\left(6+\mathrm{i}\frac{18}{\sqrt{3}}\right)\,.$ (21)
$P_{1}(z)=2k\mathrm{Re}\left(\frac{w_{21}}{z-c_{1}}\right),\
c_{1}=-3+\sqrt{3}\mathrm{i},\
w_{21}=-\left(\frac{1}{2}+\mathrm{i}\frac{5\sqrt{3}}{6}\right)\,.$ (22)
$P_{2}(z)=4k\mathrm{Re}\left(\frac{w_{31}}{z-c_{1}}\right),\
c_{1}=-3+\sqrt{3}\mathrm{i},\ w_{31}=-\mathrm{i}\frac{\sqrt{3}}{6}.$ (23)
$P_{3}(z)=2k\mathrm{Re}\left(\frac{w_{41}}{z-c_{1}}\right),\
c_{1}=-3+\mathrm{i}\sqrt{3},\
w_{41}=\frac{1}{2}+\mathrm{i}\frac{\sqrt{3}}{6}.$ (24)
$\tilde{P}(z)=48k\mathrm{Re}\left(\frac{w_{51}}{z-c_{2}}\right),\
c_{2}=-6+2\mathrm{i}\sqrt{3},\ w_{51}=-\mathrm{i}\frac{\sqrt{3}}{12}.$ (25)

Partial fraction decomposition of rational approximations to the exponential
was also used to develop the ETD Crank Nicolson scheme [13] and ETDRDP scheme
[17], where in the later all poles and weights are real.

####  2.2.4 Implementation of ETDRK4P22 scheme

Before laying out the implementation steps for the new ETDRK4P22-IF scheme, it
is instructive to first describe the implementation of the unsplit ETDRK4P22
scheme as some expressions will be carried over to shorten the discussion.

The ETDRK4P22 scheme expressed in terms of the respective rational
approximations of the matrix exponential is

$\displaystyle a_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA)U_{n}+\tilde{P}(kA)F(U_{n},t_{n})$
$\displaystyle b_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA)U_{n}+\tilde{P}(kA)F(a_{n},t_{n}+\frac{k}{2})$
$\displaystyle c_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA)a_{n}+\tilde{P}(kA)[2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n})]\,$
$\displaystyle U_{n+1}$
$\displaystyle=R_{2,2}(kA)U_{n}+P_{1}(kA)F(U_{n},t_{n})+2P_{2}(kA)G(a_{n},b_{n},t_{n}+\frac{k}{2})$
(26) $\displaystyle+P_{3}(kA)F(c_{n},t_{n}+k)\,,$

To implement this scheme, we first need to compute $a_{n}$, $b_{n}$ and
$c_{n}$. Substituting the decomposition in (21) and (25) into (2.2.4), $a_{n}$
becomes

$\displaystyle a_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA)U_{n}+\tilde{P}(kA)F(U_{n},t_{n})$
$\displaystyle=\left(I+4\mathrm{Re}\left(w_{11}(kA-
c_{2}I)^{-1}\right)\right)U_{n}+48k\mathrm{Re}\left(w_{51}(kA-
c_{2}I)^{-1}\right)F(U_{n},t_{n})$ $\displaystyle=U_{n}+2\mathrm{Re}\left[(kA-
c_{2}I)^{-1}(2w_{11}U_{n})\right]+2k\mathrm{Re}\left[(kA-
c_{2}I)^{-1}(24w_{51}F(U_{n},t_{n}))\right]$
$\displaystyle=U_{n}+2\mathrm{Re}(a_{n1})\,,$ (27)

where $(kA-c_{2}I)a_{n1}=2w_{11}U_{n}+k24w_{51}F(U_{n},t_{n})$. For $b_{n}$ we
get

$\displaystyle b_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA)U_{n}+\tilde{P}(kA)F(a_{n},t_{n}+\frac{k}{2})$
$\displaystyle=U_{n}+\mathrm{Re}\left[(kA-
c_{2}I)^{-1}(4w_{11}U_{n})\right]+k\mathrm{Re}\left[(kA-
c_{2}I)^{-1}(48w_{51}F(a_{n},t_{n}+\frac{k}{2}))\right]$
$\displaystyle=U_{n}+2\mathrm{Re}(b_{n1})\,,$ (28)

where $(kA-c_{2})b_{n1}=2w_{11}U_{n}+24w_{51}kF(a_{n},t_{n}+\frac{k}{2})$.
Finally for $c_{n}$ we have,

$\displaystyle c_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA)a_{n}+\tilde{P}(kA)[2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n})]$
$\displaystyle=\left(I+4\mathrm{Re}\left(w_{11}(kA-
c_{2}I)^{-1}\right)\right)a_{n}+48k\mathrm{Re}\left(w_{51}(kA-
c_{2}I)^{-1}\right)\times$
$\displaystyle[2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n})]$
$\displaystyle=a_{n}+2\mathrm{Re}\left[(kA-
c_{2}I)^{-1}\left(2w_{11}a_{n}+24w_{51}k[2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n})]\right)\right]$
$\displaystyle=a_{n}+2\mathrm{Re}(c_{n1})\,,$ (29)

where $(kA-
c_{2}I)c_{n1}=2w_{11}a_{n}+24w_{51}k(2F(b_{n},t_{n}+\frac{k}{2})-F(U_{n},t_{n}))$.

For the main scheme we make use of the decomposition (20), (22)–(24) to
obtain,

$\displaystyle U_{n+1}$
$\displaystyle=R_{2,2}(kA)U_{n}+P_{1}(kA)F(U_{n},t_{n})+8k\mathrm{Re}(w_{31}(kA-
c_{1}I)^{-1})G(a_{n},b_{n},t_{n}+\frac{k}{2})$
$\displaystyle+2k\mathrm{Re}(w_{41}(kA-c_{1}I)^{-1})F(c_{n},t_{n}+k)$
$\displaystyle=\left(I+2\mathrm{Re}\left(w_{11}(kA-
c_{1}I)^{-1}\right)\right)U_{n}+2k\mathrm{Re}\left(w_{21}(kA-
c_{1}I)^{-1}\right)F(U_{n},t_{n})$ $\displaystyle+2k\mathrm{Re}\left[(kA-
c_{1}I)^{-1}(4w_{31}G(a_{n},b_{n},t_{n}+\frac{k}{2})+w_{41}F(c_{n},t_{n}+k))\right]$
$\displaystyle=U_{n}+2\mathrm{Re}(U_{n1})\,,$ (30)

where $(kA-
c_{1}I)U_{n1}=w_{11}U_{n}+kw_{21}F(U_{n},t_{n})+4w_{31}kG(a_{n},b_{n},t_{n}+\frac{k}{2})+w_{41}kF(c_{n},t_{n}+k))$.
Putting all the results together, the ETDRK4P22 scheme can be implemented in
serial with the following steps:

  1. 1.

Solve: $(kA-c_{2}I)a_{n1}=2w_{11}u_{n}+24w_{51}kF(u_{n},t_{n})$

  2. 2.

Set: $a_{n}=u_{n}+2\mathrm{Re}(a_{n1})$

  3. 3.

Solve: $(kA-c_{2}I)b_{n1}=2w_{11}u_{n}+24w_{51}kF(a_{n},t_{n}+\frac{k}{2})$

  4. 4.

Set: $b_{n}=u_{n}+2\mathrm{Re}(b_{n1})$

  5. 5.

Solve: $(kA-
c_{2}I)c_{n1}=2w_{11}a_{n}+24w_{51}k[2F(b_{n},t_{n}+\frac{k}{2})-F(u_{n},t_{n})]$

  6. 6.

Set: $c_{n}=a_{n}+2\mathrm{Re}(c_{n1})$

  7. 7.

Solve: $(kA-
c_{1}I)u_{n1}=w_{11}u_{n}+w_{21}kF(u_{n},t_{n})+4w_{31}kG(a_{n},b_{n},t_{n}+\frac{k}{2})+w_{41}kF(c_{n},t_{n}+k)$.

  8. 8.

Set: $u_{n+1}=u_{n}+2\mathrm{Re}(u_{n1})$.

####  2.2.5 Implementation of ETDRK4P22-IF scheme

Following from the derivation above, we lay out the implementation of the
dimensionally split scheme begining from the intermediate solutions $a_{n}$,
$b_{n}$ and $c_{n}$. For $a_{n}$ we make use of (27) to obtain,

$\displaystyle a_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA_{1})\left[\tilde{R}_{2,2}(kA_{2})U_{n}+\tilde{P}(kA_{2})F(U_{n},t_{n})\right]$
$\displaystyle=\tilde{R}_{2,2}(kA_{1})(u_{n}+2\mathrm{Re}(a_{n1}))$ where
$(kA_{2}-c_{2}I)a_{n1}=2w_{11}U_{n}+k24w_{51}F(U_{n},t_{n})$. Setting
$a_{n2}=U_{n}+2\mathrm{Re}(a_{n1})$ and applying the partial fraction
decomposition of $\tilde{R}_{2,2}(kA_{1})$, we have $\displaystyle a_{n}$
$\displaystyle=\left(I+4\mathrm{Re}\left(w_{11}(kA_{1}-c_{2}I)^{-1}\right)\right)a_{n2}$
$\displaystyle=a_{n2}+2\mathrm{Re}((kA_{1}-c_{2}I)^{-1}(2w_{11}a_{n2}))$
$\displaystyle=a_{n2}+2\mathrm{Re}(a_{n3}),$

with $(kA_{1}-c_{2}I)a_{n3}=2w_{11}a_{n2}$. Similarly for $b_{n}$ we use (28)
to obtain

$\displaystyle b_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA_{1})\tilde{R}_{2,2}(kA_{2})U_{n}+\tilde{P}(kA_{2})F(a_{n},t_{n}+\frac{k}{2})$
$\displaystyle=\tilde{R}_{2,2}(kA_{1})(U_{n}+2\mathrm{Re}(b_{n1}))+2\mathrm{Re}(b_{n2})),$
where $(kA_{2}-c_{2}I)b_{n1}=2w_{11}U_{n}$ and
$(kA_{2}-c_{2}I)b_{n2}=24kw_{51}F(a_{n},t_{n}+\frac{k}{2})$. Setting
$b_{n3}=U_{n}+2\mathrm{Re}(b_{n1})$ and applying the partial fraction
decomposition of $\tilde{R}_{2,2}(kA_{1})$ we see that $\displaystyle b_{n}$
$\displaystyle=\left(I+4\mathrm{Re}\left(w_{11}(kA_{1}-c_{2}I)^{-1}\right)\right)b_{n3}+2\mathrm{Re}(b_{n2})$
$\displaystyle=b_{n3}+2\mathrm{Re}((kA_{1}-c_{2}I)^{-1}(2w_{11}b_{n3}))+2\mathrm{Re}(b_{n2})$
$\displaystyle=b_{n3}+2\mathrm{Re}(b_{n4})+2\mathrm{Re}(b_{n2}),$

with $(kA_{1}-c_{2}I)b_{n4}=2w_{11}b_{n3}$. Similarly from (29),

$\displaystyle c_{n}$
$\displaystyle=\tilde{R}_{2,2}(kA_{1})[\tilde{R}_{2,2}(kA_{2})a_{n}+2\tilde{P}(kA_{2})F(b_{n},t_{n}+\frac{k}{2})]-R_{2,2}(kA_{1})\tilde{P}(kA_{2})F(U_{n},t_{n})]$
$\displaystyle=\tilde{R}_{2,2}(kA_{1})c_{n1}^{*}-R_{2,2}(kA_{1})c_{n2}^{*},$

where
$c_{n1}^{*}=\tilde{R}_{2,2}(kA_{2})a_{n}+2\tilde{P}(kA_{2})F(b_{n},t_{n}+\frac{k}{2})$
and $c_{n2}^{*}=\tilde{P}(kA_{2})F(U_{n},t_{n})$. Using the respective partial
fraction decomposition we get,

$c_{n}=c_{n1}^{*}+2\mathrm{Re}(c_{n3})-[c_{n2}^{*}+2\mathrm{Re}(c_{n4})]$

where $(kA_{1}-c_{2}I)c_{n3}=2w_{11}c_{n1}^{*}$ and
$(kA_{1}-c_{1}I)c_{n4}=w_{11}c_{n2}^{*}$. Now it remains to find $c_{n1}^{*}$
and $c_{n2}^{*}$. We know that

$c_{n1}^{*}=\tilde{R}_{2,2}(kA_{2})a_{n}+2\tilde{P}(kA_{2})F(b_{n},t_{n}+\frac{k}{2})=a_{n}+2\mathrm{Re}(c_{n1})$

with $(kA_{2}-c_{2}I)c_{n1}=2w_{11}a_{n}+48w_{51}kF(b_{n},t_{n}+\frac{k}{2})$.
Similarly we obtain,

$c_{n2}^{*}=\tilde{P}(kA_{2})F(U_{n},t_{n})=2\mathrm{Re}(c_{n2})\,,$

where $(kA_{2}-c_{2}I)c_{n2}=24w_{51}kF(U_{n},t_{n})$. Now for the main
scheme, set $U_{n1}^{*}=R_{2,2}(kA_{2})U_{n}+P_{1}(kA_{2})F(U_{n},t_{n})$,
$U_{n2}^{*}=2P_{2}(kA_{2})G(a_{n},b_{n},t_{n}+\frac{k}{2})$ and
$U_{n3}^{*}=P_{3}(kA_{2})F(c_{n},t_{n}+k)$, then we can express $U_{n+1}$ as

$\displaystyle U_{n+1}$
$\displaystyle=R_{2,2}(kA_{1})U_{n1}^{*}+\tilde{R}_{2,2}(kA_{1})U_{n2}^{*}+U_{n3}^{*}$
$\displaystyle=\left(I+2\mathrm{Re}\left(w_{11}(kA_{1}-c_{1}I)^{-1}\right)\right)U_{n1}^{*}+\left(I+4\mathrm{Re}\left(w_{11}(kA_{1}-c_{2}I)^{-1}\right)\right)U_{n2}^{*}+U_{n3}^{*}$
$\displaystyle=U_{n1}^{*}+U_{n2}^{*}+2\mathrm{Re}((kA_{1}-c_{1}I)^{-1}(w_{11}U_{n1}^{*}))+2\mathrm{Re}((kA_{1}-c_{2}I)^{-1}(2w_{11}U_{n2}^{*}))+U_{n3}^{*}$
$\displaystyle=U_{n1}^{*}+U_{n2}^{*}+2\mathrm{Re}(U_{n4})+2\mathrm{Re}(U_{n5})+U_{n3}^{*}\,,$

where $(kA_{1}-c_{1}I)U_{n4}=w_{11}U_{n1}^{*}$ and
$(kA_{1}-c_{2}I)U_{n5}=2w_{11}U_{n2}^{*}$. Now from (30)
$U_{n1}^{*}=U_{n}+2\mathrm{Re}(U_{n1})$ where
$(kA_{2}-c_{1}I)U_{n1}=w_{11}U_{n}+kw_{21}F(U_{n},t_{n}).$ Similarly,
$U_{n2}^{*}=2\mathrm{Re}(U_{n2})$ where
$(kA_{2}-c_{1}I)U_{n2}=4w_{31}kG(a_{n},b_{n},t_{n}+\frac{k}{2}).$ Finally,
$U_{n3}^{*}=2\mathrm{Re}(U_{n3})$ where
$(kA_{2}-c_{1}I)U_{n3}=w_{41}kF(c_{n},t_{n}+k)$.

#### Algorithm for 2D ETDRK4P22-IF

The full ETDRK4P22-IF scheme can be implemented in parallel with up to three
processors in the following steps.

  1. 1.

Solve: $(kA_{2}-c_{2}I)a_{n1}=2w_{11}U_{n}+24w_{51}kF(U_{n},t_{n})$

  2. 2.

Set: $a_{n2}=U_{n}+2\mathrm{Re}(a_{n1})$

  3. 3.

Solve: $(kA_{1}-c_{2}I)a_{n3}=2w_{11}a_{n2}$

  4. 4.

Set: $a_{n}=a_{n2}+2\mathrm{Re}(a_{n3})$.

  5. 5.

Solve Processor 1: $(kA_{2}-c_{2}I)b_{n1}=2w_{11}U_{n}$

  6. 6.

Solve Processor 2: $(kA_{2}-c_{2}I)b_{n2}=24w_{51}kF(a_{n},t_{n}+\frac{k}{2})$

  7. 7.

Set: $b_{n3}=U_{n}+2\mathrm{Re}(b_{n1})$

  8. 8.

Solve: $(kA_{1}-c_{2}I)b_{n4}=2w_{11}b_{n3}$

  9. 9.

Set: $b_{n}=b_{n3}+2\mathrm{Re}(b_{n4})+2\mathrm{Re}(b_{n2})$.

  10. 10.

Solve Processor 1:
$(kA_{2}-c_{2}I)c_{n1}=2w_{11}a_{n}+48w_{51}kF(b_{n},t_{n}+\frac{k}{2})$

  11. 11.

Solve Processor 2: $(kA_{2}-c_{2}I)c_{n2}=24w_{51}kF(U_{n},t_{n})$

  12. 12.

Set: $c_{n1}^{*}=a_{n}+2\mathrm{Re}(c_{n1})$ and
$c_{n2}^{*}=2\mathrm{Re}(c_{n2})$

  13. 13.

Solve Processor 1: $(kA_{1}-c_{2}I)c_{n3}=2w_{11}c_{n1}^{*}$

  14. 14.

Solve Processor 2: $(kA_{1}-c_{1}I)c_{n4}=w_{11}c_{n2}^{*}$

  15. 15.

Set: $c_{n}=c_{n1}^{*}+2\mathrm{Re}(c_{n3})-[c_{n2}^{*}+2\mathrm{Re}(c_{n4})]$

  16. 16.

Solve Processor 1: $(kA_{2}-c_{1}I)U_{n1}=w_{11}U_{n}+w_{21}kF(U_{n},t_{n})$

  17. 17.

Solve Processor 2:
$(kA_{2}-c_{1}I)U_{n2}=4w_{31}kG(a_{n},b_{n},t_{n}+\frac{k}{2})$

  18. 18.

Solve Processor 3: $(kA_{2}-c_{1}I)U_{n3}=w_{41}kF(c_{n},t_{n}+k)$

  19. 19.

Set: $u_{n1}^{*}=u_{n}+2\mathrm{Re}(U_{n1})$,
$U_{n2}^{*}=2\mathrm{Re}(U_{n2})$ $U_{n3}^{*}=2\mathrm{Re}(U_{n3})$

  20. 20.

Solve Processor 1: $(kA_{1}-c_{1}I)U_{n4}=w_{11}U_{n1}^{*}$

  21. 21.

Solve Processor 2: $(kA_{1}-c_{2}I)U_{n5}=2w_{11}U_{n2}^{*}$

  22. 22.

Set:
$U_{n+1}=U_{n1}^{*}+U_{n2}^{*}+U_{n3}^{*}+2\mathrm{Re}(U_{n4})+2\mathrm{Re}(U_{n5})$.

For non-smooth problems the following third order scheme, developed in [16],
can be used for presmoothing. We show in our numerical experiments that
applying a few steps of this presmoothing scheme followed by ETDRK4P22-IF
preserves the fourth order accuracy of ETDRK4P22-IF.

  1. 1.

Solve Processor 1: $(kA-f_{1}I)a_{n1}=2s_{11}u_{n}+ks_{51}F(u_{n},t_{n})$

  2. 2.

Solve Processor 2: $(kA-f_{2}I)a_{n2}=2s_{12}u_{n}+ks_{52}F(u_{n},t_{n})$

  3. 3.

Set $a_{n}=a_{n1}+2\mathrm{Re}(a_{n2})$.

  4. 4.

Solve Processor 1: $(kA-
f_{1}I)b_{n1}=2s_{11}u_{n}+ks_{51}F(a_{n},t_{n}+\frac{k}{2})$

  5. 5.

Solve Processor 2: $(kA-
f_{2}I)b_{n2}=2s_{12}u_{n}+ks_{52}F(a_{n},t_{n}+\frac{k}{2})$

  6. 6.

Set $b_{n}=b_{n1}+2\mathrm{Re}(b_{n2})$.

  7. 7.

Solve Processor 1: $(kA-
f_{1}I)c_{n1}=2s_{11}a_{n}+ks_{51}(2F(b_{n},t_{n}+\frac{k}{2})-F(u_{n},t_{n}))$

  8. 8.

Solve Processor 2: $(kA-
f_{2}I)c_{n2}=2s_{12}a_{n}+ks_{52}(2F(b_{n},t_{n}+\frac{k}{2})-F(u_{n},t_{n}))$.

  9. 9.

Set $c_{n}=c_{n1}+2\mathrm{Re}(c_{n2})$.

  10. 10.

Solve Processor 1: $(kA-
e_{1}I)u_{n1}=s_{11}u_{n}+ks_{21}F(u_{n},t_{n})+2ks_{31}G(a_{n},b_{n},t_{n}+\frac{k}{2})+ks_{41}F(c_{n},t_{n}+k)$.

  11. 11.

Solve Processor 2 $(kA-
e_{2}I)u_{n2}=s_{12}u_{n}+ks_{22}F(u_{n},t_{n})+2ks_{32}G(a_{n},b_{n},t_{n}+\frac{k}{2}))+ks_{42}F(c_{n},t_{n}+k)$.

  12. 12.

Set $u_{n+1}=u_{n1}+2\mathrm{Re}(u_{n2})$.

The poles and weights are given by

$\displaystyle f_{1}$ $\displaystyle=-3.19214327596664;\
f_{2}=-1.40392836201668+3.61467898890404\mathrm{i}$ $\displaystyle e_{1}$
$\displaystyle=0.5f_{1};\ e2=0.5f_{2};\
s_{11}=\frac{6}{|e_{1}-e_{2}|^{2}};s_{12}=\frac{-3\mathrm{i}}{(\mathrm{Im}(e_{2}))(e_{2}-e_{1})}$
$\displaystyle s_{21}$
$\displaystyle=\frac{1-e_{1}}{|e_{1}-e_{2}|^{2}};s_{22}=\frac{\mathrm{i}(e_{2}-1)}{2(\mathrm{Im}(e_{2}))(e_{2}-e_{1})};\
s_{31}=\frac{1+e_{1}}{|e_{1}-e_{2}|^{2}};s_{32}=\frac{-\mathrm{i}(e_{2}+1)}{2(\mathrm{Im}(e_{2}))(e_{2}-e_{1})}$
$\displaystyle s_{41}$
$\displaystyle=\frac{1+e_{1}^{2}}{|e_{1}-e_{2}|^{2}};s_{42}=\frac{-\mathrm{i}(e_{2}^{2}+1)}{2(\mathrm{Im}(e_{2}))(e_{2}-e_{1})};\
s_{51}=\frac{24+6f_{1}+f_{1}^{2}}{|f_{1}-f_{2}|^{2}};s_{52}=\frac{-\mathrm{i}(24+6f_{2}+f_{2}^{2})}{2(\mathrm{Im}(f_{2}))(f_{2}-f_{1})}.$

##  3 Results and Discussion

We investigate the order of accuracy and efficiency of the proposed scheme
ETDRK4P22-IF, which applies integrating factor-type dimensional splitting to
improve the run time of the existing fourth-order ETDRK4P22 scheme developed
in [16]. We also compare the accuracy of the new scheme with the existing
second-order ETDRDP-IF developed in [19] with dimensional splitting as well as
the fourth-order semi-implicit backward differentiation scheme (SBDF4)
developed in [23]. SBDF4 is an IMEX scheme which employs a fourth-order
backward differentiation formula to discretize the linear part of (2) and a
fourth order Adam’s Bashford scheme to discretize the non-linear part. The
resulting scheme is of the form

$\begin{split}(25I+13kA)U_{n+1}&=48U_{n}-36U_{n-1}+16U_{n-2}-3U_{n-3}\\\
&+k\left(48F(U_{n},t_{n})-72kF(U_{n-1},t_{n-1})\right.\\\
&+\left.48kF(U_{n-2},t_{n-2})-12kF(U_{n-3},t_{n-3})\right).\end{split}$ (31)

The initial stages of the scheme $U_{1},U_{2},U_{3}$ are computed using the
first-order SBDF scheme

$(I+k_{0}A)U_{n+1}=U_{n}+k_{0}F(U_{n},t_{n})$

with a very small time step of $k_{0}=\frac{k}{2000}$, as was previously
suggested in [23, 22].

We evaluate the performance of the schemes on problems in one and two spatial
dimensions as well as for scalar and systems of reaction-diffusion equations
with Dirichlet and Neumann boundary conditions. The error, $E(k)$ at time step
$k$, is measured with the $L_{\infty}$ norm and the estimated rate of
convergence is calculated with the formula

$p=\frac{\log[E(k)/E(k/2)]}{\log(2)}$ (32)

For problems with exact solution, $E(k)=\|U(k)-\hat{U}(k)\|$, where
$U(k),\hat{U}(k)$ represent the exact and approximate grid solutions at the
final time $T$ using a temporal step size of $k$. All spatial derivatives are
approximated with a fourth-order finite difference scheme on a mesh with size
$h$ which enabled us to derive the correct order of convergence of the fully
discrete scheme by setting $k=h$. For problems where the exact solution is not
known, we only evaluated the convergence of the time discretization by fixing
the spatial mesh size, $h$. The errors are computed by using the numerical
solution on the next fine mesh as the reference solution. For a grid
refinement study with temporal step sizes $k,k/2,k/4$, the formula

$\tilde{p}=\frac{\log[\tilde{E}(k)/\tilde{E}(k/2)]}{\log(2)}$ (33)

is used, where $\tilde{E}(k)=\|\hat{U}(k)-\hat{U}(k/2)\|$ has been shown in
[24] to determine the correct rate of convergence of the scheme for
sufficiently small $k$. Each numerical scheme evaluated was used to calculate
its own reference solution. Numerical experiments were run in MATLAB R2020b on
a THINKMATE Work Station with Intel (R) Core(TM) i9-10900k CPU @ 3.70 GHz with
128 GB RAM. All codes used to generate convergence tables and figures have
been made publicly available at https://github.com/easantea/ETDRK4P22-IF.

###  3.1 Model problem with Dirichlet boundary

In order to investigate the order of convergence of the proposed scheme, we
consider the following two-dimensional reaction-diffusion equation with exact
solution

$\frac{\partial u}{\partial t}=\Delta u-u\,,\quad-\frac{\pi}{2}<x\,,\quad
y<\frac{\pi}{2}\,,\quad t\in[0,T]\,.$ (34)

with homogeneous Dirichlet boundary conditions and initial condition

$u(x,y,0)=\cos(x)\cos(y).$

The exact solution is $u(x,y,t)=\mathrm{e}^{-3t}\cos(x)\cos(y)$.

Table 1: Convergence table showing fourth-order convergence of the new
ETDRK4P22-IF compared with the existing ETDRK4P22 scheme for the model test
problem with exact solution at final time $T=1$

|  | ETDRK4P22-IF | ETDRK4P22
---|---|---|---
$k$ | $h$ | Error | $p$ | Cpu Time (sec) | Error | $p$ | Cpu Time (sec)
0.10 | 0.08 | 1.639$\times 10^{-7}$ | – | 0.06 | 9.069$\times 10^{-7}$ | – | 0.02
$0.05$ | 0.04 | 1.0805 $\times 10^{-8}$ | 3.92 | 0.08 | 5.6131$\times 10^{-8}$ | 4.01 | 0.35
$0.025$ | 0.02 | 6.958$\times 10^{-10}$ | 3.96 | 0.65 | 3.496$\times 10^{-9}$ | 4.01 | 5.21
$0.0125$ | 0.010 | 4.456$\times 10^{-11}$ | 3.96 | 3.72 | 2.1391$\times 10^{-10}$ | 4.03 | 79.70

We discretize the spatial derivative in each dimension using the following
fourth-order central difference approximation described in Section 2.1, where
a fourth degree polynomial extrapolation is used to compute the derivatives
near the boundary nodes ($j=0,j=m+1$). Each dimension of our spatial domain is
descritized with $m+2$ nodes with $h=\frac{\pi}{m+1}$. Linear systems in the
numerical scheme are solved using Gaussian elimination with LU-decomposition.
From Table 1 and Fig. 1A we see that proposed dimensional splitting scheme
ETDRK4P22-IF is fourth-order accurate. For a serial implementation of both
schemes, ETDRK4P22-IF is 20 times faster than ETDRK4P22 at the finest spatial
mesh size ($h=0.010$) and temporal step size ($k=0.0125$) (Table 1).
Interestingly, the new ETDRK4P22-IF is more accurate than the existing
ETDRK4P22 scheme across all temporal and spatial step sizes examined (Fig 1B,
Table 1). Whereas ETDRK4P22-IF requires more CPU time per time step compared
to the second-order ETDRDP-IF, the errors for ETDRK4P22-IF are several orders
of magnitude smaller than those for ETDRDP-IF (Fig. 1B, Table A.1).

For all time steps examined, the SBDF4 scheme is less accurate and requires
substantially greater CPU time compared with ETDRK4P22-IF (Fig. 1B, Table
A.1). Most of this time is spent computing the starting values of the
multistep method. Increasing the time step used in computing the starting
values decreased the CPU time but also worsened the accuracy of the scheme.

Figure 1: Convergence and efficiency plots for Example 3.1. A) Fourth-order
convergence and B) Efficiency of the new fourth-order ETDRK4P22-IF scheme
compared with the existing fourth-order ETDRK4P22, second-order ETDRDP-IF and
the fourth order SBDF4 schemes.

###  3.2 Model problem with Neumann boundary

Next, we investigate the order of convergence of the proposed scheme for a
scalar reaction-diffusion problem with homogeneous Neumann boundary
conditions. The model problem is given by

$\frac{\partial u}{\partial t}=\Delta u-u\,,\quad-\pi<x\,,\quad y<\pi\,,\quad
t\in[0,T]\,.$ (35)

with homogeneous Neumann boundary conditions and initial condition

$u(x,y,0)=\cos(x)\cos(y)\,{\color[rgb]{1,0,0}.}$

The exact solution is $u(x,y,t)=\mathrm{e}^{-3t}\cos(x)\cos(y)$.

Table 2: Convergence table showing fourth-order convergence of the new
ETDRK4P22-IF scheme compared with the existing ETDRK4P22 scheme for the model
test problem in Example 3.2 with $T=1$.

|  | ETDRK4P22-IF | ETDRK4P22
---|---|---|---
$k$ | $h$ | Error | $p$ | Cpu Time (sec) | Error | $p$ | Cpu Time (sec)
0.10 | 0.1 | 1.0836$\times 10^{-5}$ | – | 0.003 | 1.1580$\times 10^{-5}$ | – | 0.002
$0.05$ | 0.05 | 6.8127 $\times 10^{-7}$ | 3.99 | 0.012 | 7.2661$\times 10^{-7}$ | 3.99 | 0.042
$0.025$ | 0.025 | 4.2638$\times 10^{-8}$ | 4.00 | 0.081 | 4.5439$\times 10^{-8}$ | 4.00 | 0.732
$0.0125$ | 0.0124 | 2.6657$\times 10^{-9}$ | 4.00 | 1.961 | 2.8397$\times 10^{-9}$ | 4.00 | 10.703

Figure 2: Convergence and efficiency plots for Example 3.2. A) Exact solution,
B) Numerical solution and log-log plots showing C) Fourth-order convergence
and D) Efficiency of the new fourth-order ETDRK4P22-IF scheme compared with
the existing fourth-order ETDRK4P22, second-order ETDRDP-IF and fourth order
SBDF4 schemes.

We partition the spatial domain into $m+2$ nodes with $h=\frac{2\pi}{m+1}$.
The second derivative in each dimension is discretized using the fourth-order
central difference approximation described in Section 2.1. The resulting
linear systems in the numerical scheme is solved using Gaussian elimination
with LU-decomposition. The exact and numerical solution profiles are shown in
Fig. 2(A-B). The proposed scheme converges with fourth-order accuracy (Table
2, Fig 2C). In a serial implementation, ETDRK4P22-IF is 5 times faster than
ETDRK4P22 at the finest spatial mesh size ($h=0.0124$) and temporal step size
($k=0.0125$) (Table 2). The accuracy of the proposed scheme is equivalent to
the existing fourth-order scheme with reduced computational time being the
main source of efficiency (Fig 2C, D). Interestingly, the speed of
ETDRK4P22-IF is comparable to the second-order dimensional splitting scheme
ETDRDP-IF, but its errors are three orders of magnitude smaller for each time
step (Fig 2D, Table A.2). Similarly, SBDF4 is much slower and less accurate
than the proposed scheme at all time steps (Fig 2D, Table A.2).

###  3.3 Enzyme Kinetics

We consider the two-dimensional analogue of the one-dimensional reaction-
diffusion equation with Michaelis-Menten enzyme kinetics [25, 26, 27]

$\frac{\partial u}{\partial t}=d\Delta u-\frac{u}{(1+u)}\,,\quad
0<x\,,\,y<1\,,\quad t>0$ (36)

with homogeneous Dirichlet boundary conditions and initial condition

$u(x,y,0)=\sin(\pi x)\sin(\pi y)\,,\quad 0\leq x\,,\,y\leq 1\,.$

We use this problem to evaluate the performance of the ETDRK4P22-IF scheme for
non-linear scalar RDEs with Dirichlet boundary conditions and compare its
performance to the existing ETDRK4P22 scheme and ETDRDP-IF schemes. The
spatial derivative is discretized as described in Example 3.3. We note that
this problem has no exact solution and so all errors shown are course-to-fine
grid errors as described in the introduction of this section. Our empirical
convergence studies show that the scheme is fourth-order accurate with better
accuracy than the existing fourth-order ETDRK4P22 scheme (Table 3, Fig 3A) and
the second-order ETDRDP-IF scheme (Fig. 3A, Table A.3). Whereas the new scheme
is more computational efficient (Table 2, Fig. 3B) than the existing ETDRK4P22
scheme, the savings in computational time is not as substantial as was
observed in the previous problem. This is most likely because our convergence
analysis is performed at a fixed spatial mesh size of $h=0.05$. At this
resolution, the computational cost in solving the linear systems in the scheme
is relatively small, thus providing little advantage for the dimensional
splitting. Whereas ETDRK4P22-IF is generally more accurate than ETDRK4P22, for
relatively small problems (less than $20\times 20$ system), it provides little
to no speed up in CPU time at this scale. Observe that for all time steps
SBDF4 has a substantially high computational cost compared to the other
schemes. In fact, there appears to be very little difference in the CPU time
as we decrease the time step (Fig. 3B, Table A.3). This is due to the very
small time step needed to compute the starting values of the scheme, as
pointed out in the previous examples. Thus, the one-step nature of fourth
order ETD schemes make them more computationally efficient than fourth-order
IMEX schemes.

Table 3: Convergence table showing fourth-order convergence of the new
ETDRK4P22-IF compared with the existing ETDRK4P22 scheme for the enzyme
kinetics problem.

|  | ETDRK4P22-IF | ETDRK4P22
---|---|---|---
$k$ | $h$ | Error | $\tilde{p}$ | Cpu Time (sec) | Error | $\tilde{p}$ | Cpu Time (sec)
0.10 | 0.05 | 4.2433$\times 10^{-7}$ | – | 0.003 | 1.9274$\times 10^{-6}$ | – | 0.004
0.05 | 0.05 | 7.2737 $\times 10^{-9}$ | 5.87 | 0.007 | 1.1628$\times 10^{-7}$ | 4.05 | 0.008
0.025 | 0.05 | 4.666$\times 10^{-10}$ | 3.96 | 0.012 | 7.1638$\times 10^{-9}$ | 4.02 | 0.014
0.0125 | 0.05 | 3.0407$\times 10^{-11}$ | 3.94 | 0.023 | 4.4488$\times 10^{-10}$ | 4.01 | 0.028

Figure 3: Convergence and efficiency plots for Example 3.3. Log-log plots
showing A) convergence and B) efficiency of ETDRK4P22-IF for the enzyme
kinetics problem. Errors are generated between successive numerical solutions
with varying temporal step size on a fixed spatial mesh size of $h=0.05$. The
diffusion coefficient is $d=0.25$.

###  3.4 Non-smooth enzyme kinetics with mismatched initial and boundary data

Here, we evaluate the performance of ETDRK4P22-IF for solving non-smooth
problems by presmoothing with the lower order scheme presented in Section
2.2.5. Such pre-smoothing techniques have been used to improve the damping
properties of the Crank-Nicolson scheme [28]. The test problem is similar to
the model presented in Example 3.3 with a diffusion coefficient $d=1$ and the
following initial condition

$u(x,y,0)=1.$

This initial condition is discontinuous at the boundary where homogeneous
Dirichlet boundary conditions are applied. Such mismatch of initial and
boundary data are known to generate spurious oscillations which, if not
properly damped at the initial stages of the simulation, can significantly
affect the accuracy of solutions.

Table 4: Convergence of ETDRK4P22-IF on the nonsmooth enzyme kinetics problem
with and without initial smoothing steps. ETDRK4P22-IFs denotes the scheme
with three presmoothing steps

|  | ETDRK4P22-IF | ETDRK4P22-IFs
---|---|---|---
$k$ | $h$ | Error | $\tilde{p}$ | Cpu Time (sec) | Error | $\tilde{p}$ | Cpu Time (sec)
0.10 | 0.05 | 6.1306$\times 10^{-3}$ | – | 0.002 | 1.0894$\times 10^{-9}$ | – | 0.002
0.05 | 0.05 | 2.0160 $\times 10^{-5}$ | 8.25 | 0.003 | 9.9321$\times 10^{-11}$ | 3.46 | 0.004
0.025 | 0.05 | 7.2147$\times 10^{-11}$ | 18.09 | 0.006 | 8.5536$\times 10^{-12}$ | 3.54 | 0.007
0.0125 | 0.05 | 4.7483$\times 10^{-15}$ | 13.89 | 0.013 | 6.2814$\times 10^{-13}$ | 3.77 | 0.013

Without smoothing, the ETDRK4P22-IF scheme fails to damp out spurious
oscillations for relatively large time steps (Fig. 4B, $k=0.1$) and struggles
to maintain fourth-order accuracy (Table 4). However, when a sufficiently
small time step is used ($k=0.0125$), the scheme is able to recover the
expected accuracy (Fig. 4C,D). Performing just three smoothing steps [28] is
sufficient to recover the desired solution profile at the relatively course
time step $k=0.1$ (Fig. 4A) and preserves the fourth-order accuracy of the
scheme (Table 4). The dimensionally split scheme with smoothing
(ETDRK4P22-IFs) is also more efficient than the unsplit scheme (ETDRK4P22)
with or without presmoothing for all time steps examined (Fig. 4D). Again, the
SBDF4 scheme is less efficient than all the ETD schemes considered (Fig. 4B,
Table A.4).

Figure 4: Numerical solution for Example 3.3 with mismatched initial and
boundary data. A) Solution profile using ETDRK4P22-IF with three steps of
presmoothing. B) Solution profile using ETDRK4P22-IF without presmoothing. C)
Convergence Plots. D) Efficiency Plots. All numerical solutions are generated
on a spatial grid with mesh size $h=0.05$. A time step of $k=0.1$ was used for
the solutions profiles in A and B. The diffusion coefficient for the model is
taken as $d=1$ in all cases.

###  3.5 Brusselator in 2D

Finally, we examine the performance of ETDRK4P22-IF for simulating a system of
non-linear reaction-diffusion equations, with the Brusselator model [29] as a
test case. The model is given by

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=\epsilon_{1}\Delta
u+A+u^{2}v-(B+1)u\,,$ $\displaystyle\frac{\partial v}{\partial t}$
$\displaystyle=\epsilon_{2}\Delta v+Bu-u^{2}v\,,$

in 2D, where the diffusion coefficients are $\epsilon_{1}=\epsilon_{2}=2\cdotp
10^{-3}$, and the chemical parameters $A=1$, $B=3.4$. At the boundary of the
domain homogeneous Neumann conditions are imposed.

The initial conditions are given by

$u(x,y,0)=\frac{1}{2}+y\,,\quad v(x,y,0)=1+5x\,.$

Spatial derivatives are discretized as described in Example 3.2. Again, this
problem has no exact solution and so all errors shown are course-to-fine grid
errors as described in the introduction of this section. The numerical
solution for each component is shown in Fig 5A,B. Solution profiles are
similar to those produced in [29]. Our empirical convergence studies show that
the scheme is fourth-order accurate with similar accuracy to the existing
fourth-order ETDRK4P22 scheme (Table 5, Fig 5C). However, the new scheme is
more computationally efficient than ETDRK4P22, with seven times speed up in
CPU time at the spatial mesh size of $h=0.0125$ and temporal step size of
$k=0.00625$ (Table 5, Fig. 5D). As expected, ETDRK4P22-IF is more accurate
than the existing second-order ETDRDP-IF scheme, providing errors that are
about four orders of magnitude smaller (Fig. 5B, Table A.5). ETDRK4P22-IF is
also more efficient than SBDF4 (Fig. 5B, Table A.5) .

Table 5: Convergence table showing fourth-order convergence of the new
ETDRK4P22-IF scheme compared with the existing ETDRK4P22 scheme for
Brusselator Example 3.5 with $T=2$.

|  | ETDRK4P22-IF | ETDRK4P22
---|---|---|---
$k$ | $h$ | Error | $\tilde{p}$ | Cpu Time (sec) | Error | $\tilde{p}$ | Cpu Time (sec)
0.05 | 0.0125 | 3.1532$\times 10^{-4}$ | – | 0.160 | 3.1384$\times 10^{-4}$ | – | 1.486
$0.025$ | 0.0125 | 1.7359$\times 10^{-5}$ | 4.18 | 0.303 | 1.7592$\times 10^{-5}$ | 4.16 | 2.948
$0.0125$ | 0.0125 | 1.0814 $\times 10^{-6}$ | 4.00 | 0.603 | 1.1859 $\times 10^{-6}$ | 3.89 | 5.879
$0.00625$ | 0.0125 | 6.7987$\times 10^{-8}$ | 3.99 | 1.172 | 7.6968 $\times 10^{-8}$ | 3.95 | 11.836

Figure 5: Convergence and efficiency plots for Example 3.5. A) Numerical
solution $U$, B) Numerical solution $V$, log-log plots showing C) Fourth-order
convergence, and D) Efficiency of the new fourth-order ETDRK4P22-IF scheme
compared with the existing fourth-order ETDRK4P22, second-order ETDRDP-IF and
fourth order SBDF4 schemes. Numerical solutions are generated using $40$ nodes
per spatial dimension and time step of $0.05$ with final time $T=2$.

##  4 Conclusions

In this article, we have developed a fourth-order ETD scheme with dimensional
splitting that uses a Padé (2,2) rational function to approximate the matrix
exponential. The order of convergence of the scheme has been validated using
two linear problems with exact solution and two non-linear problems with
Dirichlet and Neumann boundary conditions in 2D. In all cases the new scheme
is more computational efficient than competing schemes. We expect the same
savings to persist for problems with periodic boundary conditions with
additional savings introduced through parallelization of the new scheme.
Addition of up to three pre-smoothing steps using a third-order L-stable
scheme effectively damps out spurious oscillations arising from mismatched
initial and boundary data and preserves the fourth-order accuracy of the
scheme. Future work will examine 3D applications and investigate the
assumptions on the differential operator and non-linear source term that
ensure fourth-order accuracy and positivity of the numerical solution.

## Declaration of Interest Statement

None.

## Author Contributions Statement

E.O. Asante-Asamani: Conception and design, Analysis and interpretation of
data, Writing original draft; B.A. Wade: Analysis and interpretation of data,
Writing-review & editing. A. Kleefeld: Analysis and interpretation of data,
Writing-review & editing.

## References

  * [1] R. C. Mittal, R. Goel, N. Ahlawat, An Efficient Numerical Simulation of a Reaction-Diffusion Malaria Infection Model using B-spline Collocation, Chaos. Soliton. Fract. 143 (2021) 110566.
  * [2] S. Kondo, T. Miura, Reaction-diffusion model as a framework for understanding biological pattern formation, Science 329 (2010) 1616.  doi:DOI:10.1126/science.1179047.
  * [3] A. N. Landge, B. M. Jordan, X. Diego, P. Müller, Pattern formation mechanisms of self-organizing reaction-diffusion systems, Dev. Biol. 460 (1) (2020) 2–11.
  * [4] T. Ghafouri, Z. G. Bafghi, N. Nouri, N. Manavizadeh, Numerical solution of the Schrödinger equation in nanoscale side-contacted FED applying the finite-difference method, Results. Phys. 19 (2020) 103502.
  * [5] N. Sundararajan, S. Sankaran, Groundwater modeling of Musi basin Hyderabad, India: a case study, Appl. Water Sci. 10 (1) (2020) 1–22.
  * [6] N. Rambeerich, D. Y. Tangman, A. Gopaul, M. Bhuruth, Exponential time integration for fast finite element solutions of some financial engineering problems, J. Comput. Appl. Math. 224 (2) (2009) 668–678.
  * [7] W. Hundsdorfer, J. G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer-Verlag Berlin, 2003.
  * [8] A. Iserles, A first course in the numerical analysis of differential equations, no. 44, Cambridge University Press, 2009.
  * [9] L. F. Shampine, Numerical Solution of Ordinary Differential Equations, Taylor and Francis, 1994.
  * [10] J. Martín-Vaquero, A. Kleefeld, Eserk5: A fifth-order extrapolated stabilized explicit runge–kutta method, J.Comput. Appl. Math. 356 (2019) 22–36.
  * [11] J. Martín-Vaquero, A. Kleefeld, Solving nonlinear parabolic pdes in several dimensions: parallelized eserk codes, J. Comput. Phys. 423 (2020) 109771\.
  * [12] S. M. Cox, P. C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys. 176 (2002) 430–455.
  * [13] B. Kleefeld, A. Q. M. Khaliq, B. A. Wade, An ETD Crank-Nicolson method for reaction-diffusion systems, Numer. Methods Partial Differ. Equ. 28 (2012) 1309–1335.
  * [14] M. Yousuf, A. Q. M. Khaliq, B. Kleefeld, The numerical approximation of nonlinear Black–Scholes model for exotic path–dependent American options with transaction cost, Int. J. Comput. Math 89 (9) (2012) 1239–1254.
  * [15] M. Yousuf, Efficient L-stable method for parabolic problems with application to pricing American options under stochastic volatility, Appl. Math. Comput. 213 (1) (2009) 121–136.
  * [16] A. Q. M. Khaliq, J. Martín-Vaquero, B. A. Wade, M. Yousuf, Smoothing schemes for reaction-diffusion systems with nonsmooth data, J. Comput. Appl. Math. 223 (1) (2009) 374–386.
  * [17] E. O. Asante-Asamani, A. Q. M. Khaliq, B. A. Wade, A real distinct poles exponential time differencing scheme for reaction–diffusion systems, J. Comput. Appl. Math. 299 (2016) 24–34.
  * [18] E. O. Asante-Asamani, B. A. Wade, A dimensional splitting of ETD schemes for reaction-diffusion systems, Commun. Comput. Phys. 19 (5) (2016) 1343–1356.
  * [19] E. O. Asante-Asamani, A. Kleefeld, B. A. Wade, A second-order exponential time differencing scheme for non-linear reaction-diffusion systems with dimensional splitting, J. Comput. Phys. 415 (2020) 109490.
  * [20] J. H. Mathews, K. D. Fink, et al., Numerical Methods Using MATLAB, 4e, Vol. 4, Pearson Prentice Hall, 2004.
  * [21] F. Gibou, R. Fedkiw, A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem, J. Comput. Phys. 202 (2) (2005) 577–601.
  * [22] A. Kassam, L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput. 26 (4) (2005) 1214–1233.
  * [23] U. M. Ascher, S. J. Ruuth, B. T. R. Wetton, Implicit-explicit methods for time-dependent partial differential equations, SIAM J. Numer. Anal. 32 (3) (1995) 797–823.
  * [24] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, SIAM, 2007.
  * [25] J. Müller, C. Kuttler, Methods and Models in Mathematical Biology, Springer, 2015.
  * [26] H. P. Bhatt, A. Q. M. Khaliq, The locally extrapolated exponential time differencing LOD scheme for multidimensional reaction-diffusion systems, J. Comput. Appl. Math. 285 (2015) 256–278.  doi:10.1016/j.cam.2015.02.017.
  * [27] Y. Cherruault, M. Choubane, J. M. Valleton, J. C. Vincent, Stability and asymptotic behavior of a numerical solution corresponding to a diffusion-reaction equation solved by a finite difference scheme (Crank-Nicolson), Comput. Math. Appl. 20 (11) (1990) 37–46.
  * [28] B. A. Wade, A. Q. M. Khaliq, M. Yousuf, J. Vigo-Aguiar, R. Deininger, On smoothing of the Crank-Nicolson scheme and higher order schemes for pricing barrier options, J.Comput. Appl. Math. 204 (1) (2007) 144–158.
  * [29] P. A. Zegeling, H. P. Kok, Adaptive moving mesh computations for reaction-diffusion systems, J. Comput. Appl. Math. 168 (2004) 519–528.

## Appendices

##  Appendix A Supplementary tables

Table A.1: Convergence table showing fourth-order accuracy of SBDF4 scheme and
second-order accuracy of ETDRDP-IF for Example 3.1

|  | ETDRDP-IF | SBDF4
---|---|---|---
$k$ | $h$ | Error | $p$ | Cpu Time (sec) | Error | $p$ | Cpu Time (sec)
0.1000 | 0.0785 | 3.3043$\times 10^{-4}$ | – | 0.008 | 2.2150$\times 10^{-4}$ | – | 1.145
0.0500 | 0.0393 | 8.4134$\times 10^{-5}$ | 1.97 | 0.018 | 1.2419$\times 10^{-5}$ | 4.16 | 14.823
0.0250 | 0.0196 | 2.1234 $\times 10^{-5}$ | 1.99 | 0.127 | 7.752 $\times 10^{-7}$ | 4.00 | 115.230
0.0125 | 0.0098 | 5.3332$\times 10^{-6}$ | 1.99 | 1.167 | 6.1782$\times 10^{-8}$ | 3.65 | 1209.314

Table A.2: Convergence table showing fourth-order accuracy of SBDF4 scheme and
second-order accuracy of ETDRDP-IF for Example 3.2

|  | ETDRDP-IF | SBDF4
---|---|---|---
$k$ | $h$ | Error | $p$ | Cpu Time (sec) | Error | $p$ | Cpu Time (sec)
0.1000 | 0.1496 | 4.4405$\times 10^{-4}$ | – | 0.002 | 2.3248$\times 10^{-4}$ | – | 0.189
0.0500 | 0.0766 | 1.0839$\times 10^{-4}$ | 2.03 | 0.010 | 1.3094$\times 10^{-5}$ | 4.15 | 1.198
0.0250 | 0.0388 | 2.6811 $\times 10^{-5}$ | 2.02 | 0.073 | 8.1722 $\times 10^{-7}$ | 4.00 | 14.316
0.0125 | 0.0195 | 6.6694$\times 10^{-6}$ | 2.01 | 1.809 | 6.4410$\times 10^{-8}$ | 3.67 | 114.896

Table A.3: Convergence table showing fourth-order accuracy of SBDF4 scheme and
second-order accuracy of ETDRDP-IF for Example 3.3

|  | ETDRDP-IF | SBDF4
---|---|---|---
$k$ | $h$ | Error | $\tilde{p}$ | Cpu Time (sec) | Error | $\tilde{p}$ | Cpu Time (sec)
0.1000 | 0.1653 | 9.4569$\times 10^{-5}$ | – | 0.0006 | 3.3077$\times 10^{-4}$ | – | 0.1511
0.0500 | 0.1653 | 2.4719$\times 10^{-5}$ | 1.94 | 0.0009 | 1.6554$\times 10^{-5}$ | 4.32 | 0.1437
0.0250 | 0.1653 | 6.3326$\times 10^{-6}$ | 1.96 | 0.0017 | 8.4646 $\times 10^{-7}$ | 4.29 | 0.1464
0.0125 | 0.1653 | 1.6036$\times 10^{-6}$ | 1.98 | 0.0037 | 4.2354$\times 10^{-8}$ | 4.32 | 0.1560

Table A.4: Convergence table for ETDRK4P22 and ETDRK4P22s schemes applied to
Example 3.4

|  | ETDRK4P22 | ETDRK4P22s
---|---|---|---
$k$ | $h$ | Error | $\tilde{p}$ | Cpu Time (sec) | Error | $\tilde{p}$ | Cpu Time (sec)
0.1000 | 0.1653 | 1.8184$\times 10^{-1}$ | – | 0.002 | 2.5622$\times 10^{-9}$ | – | 0.003
0.0500 | 0.1653 | 1.0646$\times 10^{-2}$ | 2.09 | 0.004 | 1.4589$\times 10^{-10}$ | 4.13 | 0.005
0.0250 | 0.1653 | 1.7765$\times 10^{-5}$ | 9.23 | 0.008 | 1.124 $\times 10^{-11}$ | 3.70 | 0.012
0.0125 | 0.1653 | 1.3239$\times 10^{-12}$ | 23.68 | 0.016 | 7.9819$\times 10^{-13}$ | 3.82 | 0.016

Table A.5: Convergence table showing fourth-order accuracy of SBDF4 scheme and
second-order accuracy of ETDRDP-IF for Example 3.5

|  | ETDRDP-IF | SBDF4
---|---|---|---
$k$ | $h$ | Error | $\tilde{p}$ | Cpu Time (sec) | Error | $\tilde{p}$ | Cpu Time (sec)
0.0500 | 0.0125 | 1.8794$\times 10^{-2}$ | – | 0.046 | 6.7785$\times 10^{-2}$ | – | 31.672
0.0250 | 0.0125 | 5.3233$\times 10^{-3}$ | 1.82 | 0.082 | 1.4339$\times 10^{-3}$ | 5.56 | 31.391
0.0125 | 0.0125 | 1.4270$\times 10^{-3}$ | 1.90 | 0.160 | 1.0627 $\times 10^{-4}$ | 3.75 | 33.662
0.0063 | 0.0125 | 3.7001$\times 10^{-4}$ | 1.95 | 0.310 | 6.8328$\times 10^{-6}$ | 3.96 | 34.430