# Matrix-oriented FEM formulation for stationary and time-dependent PDEs on
x-normal domains

Massimo Frittelli, Ivonne Sgura††footnotemark:  University of Salento,
Department of Mathematics and Physics “E. De Giorgi”, Via per Arnesano, 73100
Lecce, Italy. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

When numerical solution of elliptic and parabolic partial differential
equations is required to be highly accurate in space, the discrete problem
usually takes the form of large-scale and sparse linear systems. In this work,
as an alternative, for spatial discretization we provide a Matrix-Oriented
formulation of the classical Finite Element Method, called MO-FEM, of
arbitrary order $k\in\mathbb{N}$. On structured 2D domains (e.g. squares or
rectangles) the discrete problem is then reformulated as a Sylvester matrix
equation, that we solve by the _reduced approach_ in the associated spectral
space. On a quite general class of domains, namely _normal domains_ , and even
on special surfaces, the MO-FEM yields a _multiterm Sylvester matrix equation_
where the additional terms account for the geometric contribution of the
domain shape. In particular, we obtain a sequence of these equations after
time discretization of parabolic problems by the IMEX Euler method. We apply
the matrix-oriented form of the Preconditioned Conjugate Gradient (MO-PCG)
method to solve each multiterm Sylvester equation for MO-FEM of degree
$k=1,\dots,4$ and for the lumped $\mathbb{P}_{1}$ case. We choose a matrix-
oriented preconditioner with a single-term form that captures the spectral
properties of the whole multiterm Sylvester operator.
For several numerical examples, we show a gain in computational time and
memory occupation wrt the classical vector approach solving large sparse
linear systems by a direct method or by the vector PCG with same
preconditioning. As an application, we show the advantages of the MO-FEM-PCG
to approximate Turing patterns with high spatial resolution in a reaction-
diffusion PDE system for battery modeling.

## Keywords

Finite elements, Sylvester matrix equations, Reaction-diffusion, Turing
pattern, Preconditioned Conjugate Gradients

## Mathematics Subject Classification

65F45, 65M60, 65N30

##  1 Introduction

We are interested in the discretisation of (i) elliptic PDEs of the form

$-\Delta u+\gamma u=f(x,y),\qquad(x,y)\in\Omega\subset\mathbb{R}^{2},$ (1)

where $\gamma\geq 0$, (ii) parabolic PDEs of the form

$u_{t}-d_{u}\Delta u=f(x,y,t),\qquad(x,y,t)\in\Omega\times[0,T],$ (2)

with $d_{u}>0$ being a diffusion coefficient, and (iii) reaction-diffusion
systems (RDS) of the form

$\begin{cases}u_{t}-d_{u}\Delta u=f(u,v);\\\ v_{t}-d_{v}\Delta
v=g(u,v),\end{cases}\qquad(x,y,t)\in\Omega\times[0,T],$ (3)

with $d_{u},d_{v}>0$ being diffusion coeffcients. Problems (1)-(3) are endowed
with either homogeneous Dirichlet or Neumann boundary conditions, problems
(2)-(3) are endowed with suitable initial conditions. The relevance of the PDE
problems (1)-(3) is well-known, as they find numerous applications across all
fields of science. We stress that the RDS (3) is the playground of Turing’s
theory of morphogenesis [32], which encompasses extremely diverse applications
such as biological patterning [2], biomembrane modelling [11], tumour growth
[7], metal dealloying [10], financial risk management [4], oscillating
chemical reactions [33] and the recent applications to metal electrodeposition
[5] which we will consider in the present work. We focus on the approximation
of Turing patterns, because, from a computational point of view, this is a
challenging task since fine meshes are required in space to capture the
morphological class of the pattern itself (spots, labyrinths, etc) that must
be attained as steady state of the PDE dynamics for long time of integration.

Among the existing methods for the spatial discretisation of problems (1)-(3)
we mention finite differences [17], finite elements [16], spectral methods
[7], kernel methods [18] and many more. It is well known that numerical
methods typically approximate the elliptic problem (1) through an algebraic
system in vector form:

$A\boldsymbol{\xi}=\boldsymbol{b},$ (4)

with the vector $\boldsymbol{\xi}$ containing the coefficients of the
expansion of the numerical solution in a given discrete function basis, the
matrix $A$ approximating the operator $\mathcal{L}(u)=-\Delta u+\gamma u$ and
the vector $\boldsymbol{b}$ approximating the right-hand side of (1). For the
time-dependent problems (2)-(3), a common general approach is the so-called
method of lines (MOL), which consists of discretising the spatial variables
with a spatial method of choice, thereby producing a continuous-in-time ODE
system. For problem (2), in a general setting including several spatial
methods, the resulting spatially discrete problem takes the form of the
following ODE system in vector form

$\displaystyle
M\dot{\boldsymbol{\xi}}+A_{u}\boldsymbol{\xi}=\boldsymbol{b}(t),\qquad\boldsymbol{\xi}(0)=\boldsymbol{\xi}_{0},\qquad
t\in[0,T],$ (5)

with the vector $\boldsymbol{\xi}=\boldsymbol{\xi}(t)$ containing the time-
dependent coefficients of the expansion of $u$, the matrix $A_{u}$
approximating the operator $\mathcal{L}(u):=-d_{u}\Delta u$, the vector
$\boldsymbol{b}(t)$ approximating the right-hand-side of (2), and the matrix
$M$ depends on the spatial methods (e.g. the identity matrix for finite
differences, mass matrix for finite elements). Similarly, the spatially
discrete formulation of the RDS (3) becomes the following (possibly nonlinear)
ODE system in vector form

$\displaystyle\begin{cases}M\dot{\boldsymbol{\xi}}+A_{u}\boldsymbol{\xi}=\boldsymbol{b}_{1}(\boldsymbol{\xi},\boldsymbol{\eta});\\\
M\dot{\boldsymbol{\eta}}+A_{v}\boldsymbol{\eta}=\boldsymbol{b}_{2}(\boldsymbol{\xi},\boldsymbol{\eta});\\\
\boldsymbol{\xi}(0)=\boldsymbol{\xi}_{0},\quad\boldsymbol{\eta}(0)=\boldsymbol{\eta}_{0},\end{cases}\qquad
t\in[0,T].$ (6)

The main computational challenge of (4)-(6) is dimensionality. If the discrete
function space has dimension $d$, the matrices $A,A_{u},A_{v},M$ appearing in
(4)-(6) are of size $d\times d$. There are special cases, as detailed below,
where such matrices possess a general Kronecker structure with
$n\in\mathbb{N}$ terms, e.g.

$A=\sum_{i=1}^{n}R_{i}\otimes L_{i},$ (7)

with $\otimes$ denoting the Kronecker product and $R_{i},L_{i}$ being matrices
of lower dimension, e.g. $\sqrt{d}\times\sqrt{d}$, see for instance [25, 28,
26, 23, 9, 15]. In such cases, since

$A\boldsymbol{\xi}=\texttt{vec}\left(\sum_{i=1}^{n}L_{i}^{T}UR_{i}\right),$
(8)

where $U$ is such that $\texttt{vec}(U)=\boldsymbol{\xi}$, problem (4) can be
reformulated as the following linear algebraic matrix equation

$\sum_{i=1}^{n}L_{i}^{T}UR_{i}=B,$ (9)

where $B$ is such that $\texttt{vec}(B)=\boldsymbol{b}$. Problem (9) is called
a _multiterm Sylvester equation_ , see [31]. The solution of general multiterm
Sylvester equations is mostly uncharted territory, as discussed in [31, 30]
and references therein.
A special case of (9) worth mentioning is the two-term case $n=2$, when (9)
specialises to a generalised Sylvester equation

$L_{1}^{T}UR_{1}+L_{2}^{T}UR_{2}=B,$ (10)

and closed-form algorithms are available, such as the Bartels-Stewart
algorithm [3] or its improvement proposed by Golub and others [13]. If
$L_{i},R_{i}$ further fulfil suitable assumptions, even more efficient closed-
form algorithms are available, based on spectral decomposition, see for
instance [9], as we will also discuss in the next sections.
For the time-dependent problems (5)-(6), if the matrices $A_{u},A_{v},M$
possess a Kronecker decomposition similar to (7), problems (5)-(6) can be
reformulated as _matrix ODE systems_. In this work, we apply the Implicit-
Explicit (IMEX) Euler scheme directly to (6) in vector form, then we will
consider its MO counterpart. This will yield a sequence of multiterm Sylvester
matrix equations as detailed in Section 6.

Matrix formulations of spatial methods for PDEs were successfully carried out
in some notable cases. A class of elliptic problems with convection, posed on
rectangular or parallelepypedal domains, was discretised via central finite
differences in matrix-oriented form in [25] and the discrete problem takes the
form of a multiterm Sylvester equation. The methodology was then extended to
address more general polygonal domains, see [15]. Elliptic anisotropic PDEs
with stochastic terms were approximated via Galerkin method in matrix-oriented
form in [26], the discrete problem is a multiterm Sylvester equation.
Isogeometric analysis was successfully applied to various elliptic problems,
see for instance [28, 23, 1]. On square domains, the discrete problem is a
generalized Sylvester equation, see [28]. On more general domains defined
through splines or NURBS, a Sylvester form can still be achieved if using
suitable low-rank approximations of kernels, see [23]. The work in [9]
addresses time-dependent problems, specifically the heat equation and RDSs, on
rectangular domains, where the spatial discretisation is carried out via
central finite differences in matrix-oriented form and the discrete problem
takes the form of a two-term Sylvester equation, which lends itself to an
extremely efficient numerical treatment based on spectral decomposition.

In the present work we consider both elliptic PDEs of the form (1) and
parabolic PDE problems such as the semilinear heat equation (2) and RDSs (3),
posed on a class of two-dimensional spatial domains known as _normal domains_.
For such PDE problems we propose a Matrix-Oriented Finite Element Method for
the spatial discretisation, that we will define as MO-FEM. The proposed
framework advances the existing theory on matrix-oriented spatial
discretisation of PDEs in several directions, as listed below.

  * •

To the best of the author’s knowledge, the present work provides the first MO
formulation of the finite element method for elliptic and parabolic PDEs. The
proposed theory is general and applies to a large class of basis functions,
such as Lagrangian $\mathbb{P}_{k}$ basis functions, $k\in\mathbb{N}$, thereby
covering arbitrarily high-order convergence in space. Special focus is given
to the practical special case of lumped $\mathbb{P}_{1}$ finite elements.

  * •

On rectangular domains, where Cartesian-structured mesh are immediate to
construct, the discrete Laplacian takes the form of a two-term Sylvester
equation, in analogy with matrix-oriented finite differences [9] or
isogeometric analysis [28]. Moreover, since central finite differences are
equivalent to lumped $\mathbb{P}_{1}$ finite elements, the matrix-oriented FD
discretization of the considered PDE problems can be considered as a special
case of the proposed theory. In this case we solve the two term Sylvester
equation by the _reduced approach_ in the associated spectral space, see
Section 4.

  * •

Thanks to a suitable coordinate transformation, the proposed theory applies to
domains more general than rectangles, namely _normal domains_. To the best of
the authors’ knowledge, the first applications of a matrix-oriented method to
non-rectangular domains are (i) the work in [23], where isogeometric analysis
with suitable low rank approximation of kernels is applied to elliptic
problems and (ii) the work in [15], where a matrix-oriented finite difference
scheme with conformal mappings is applied to elliptic problems on _polygonal
domains_. On normal domains, the proposed approach adopts a curved mesh that
matches the (possibly curved) boundary exactly. Hence, the proposed approach
combines the low dimensionality of a matrix-oriented approach with the absence
of geometric error. Since normal domains can be wrapped around a cylinder, the
proposed method also applies to spatial domains that are special surfaces,
namely cylinders with curvilinear edges. In all these cases, he discrete
problem takes the form of a multiterm Sylvester equation, where the additional
terms account for domain shape.
In this case, for the numerical approximation of general multiterm Sylvester
equations, we propose an iterative method: a matrix-oriented preconditioned
conjugate gradient method (MO-PCG) that always converges for the considered
PDE problems,thanks to the involved differential operators being coercive and
self-adjoint, see [27, Section 6.7].

  * •

We provide numerical experiments that demonstrate that (i) for elliptic
problems on square domains, both the reduced approach and the MO-PCG exhibit
optimal spatial convergence, (ii) for elliptic and parabolic problems on
$x$-normal domains, MO-PCG exhibits optimal convergence in space (and time, if
combined with IMEX Euler), and (iii) both the reduced approach and the MO-PCG
provide a significant gain in terms of computational time and memory storage
in comparison to the standard vector form, solved both via a direct solver and
vector PCG.

  * •

Concerning reaction-diffusion systems, we show that the MO-FEM allows for
accurate simulation of Turing patterns -obtained as asymptotic solutions- that
might be prohibitive, in terms of time and memory, through standard finite
elements in vector form, because fine spatial grids are required to capture
the features of pattern morphology, after long-time integration, see [9]. In
particular, here we solve a RDS of interest for battery modeling on some
$x$-normal domains, cap and jar shaped, and cylindrical surfaces with
curvilinear boundaries. It is worth noting that the right-hand sides of such
RDSs are not low-rank, then the solution of the multiterm Sylvester equations
cannot be approximated by Krylov methods such those in [30].

The paper is organized as follows. In Section 2, we introduce preliminary
definitions and results, we elaborate on the classes of spatial domains to be
considered, and we introduce curvilinear Cartesian-structured meshes. In
Section 3, we define a general finite element method and we derive a Kronecker
decomposition of the discrete differential operators, thereby considering as a
practical variation lumped $\mathbb{P}_{1}$ finite elements.

In Section 4, we introduce the MO formulation for these FEMs. on different
kind of spatial 2D domains. Specifically,we discretise the stationary PDE
problem (1) on a square domain and we solve the corresponding Sylvester
equation by a spectral (_reduced_) approach, also in the case of Lumped FEM.
In Section 5, we present the discretisation of the elliptic problem (1) on
_$x$ -normal domains_ and we present the solution of the corresponding
multiterm Sylvester equations by the proposed matrix-oriented PCG method for
FEM of orders $k=1,\dots,4$ in space. A comparison with the vector PCG is
provided in the numerical examples.

In Section 6, we extend the proposed MO-FEM approach to the semilinear heat
equation (2) and we apply the MO-PCG method to solve the sequence of Sylvester
multiterm equations by the application of the IMEX-Euler method in time.
Specifically, in Section 6.1, we present the convergence results on a cap-
shaped domain and the computational performance in terms of execution time. In
Section 7 we present the numerical simulations of electrochemical patterns
arising in batteries for different choice of the parameters in the reaction-
diffusion system (3) corresponding to the DIB model [19, 29] yielding Turing
patterns with spots-worms and holes. We show that these patterns can be seen
as PDE solutions on cylindrical surfaces. In Section 8, we provide some
concluding remarks and highlight future research directions.

##  2 Normal domains and curvilinear structured meshes

In this section we introduce the classes of domains to be considered in this
work and we construct suitable curvilinear meshes whose nodes possess a
Cartesian ordering and that match curved boundaries exactly. We also introduce
preliminary definitions and results that we will adopt in the construction of
the proposed matrix-oriented finite element method.

###  2.1 Square domains

To make the reader familiar with the proposed setting, we start by
discretising square domains. On the one-dimensional unitary domain:
$K:=[0,1]$, we consider the equally spaced mesh $K_{h}$ with $(N+1)$ nodes,
with $N\in\mathbb{N}$. For $m=0,\dots,N$, we define the $m$-th node as
$x_{m}=m/N$. Each element of $K_{h}$ is of the form $E_{m}=[x_{m},x_{m+1}]$
for some $m=0,\dots,N-1$.
On the square two-dimensional domain: $\Omega:=[0,1]^{2}$ we consider the
Cartesian mesh $\Omega_{h}$ with $(N+1)\times(N+1)$ nodes. For
$m,n=0,\dots,N$, we define the node
$\boldsymbol{x}_{mn}:=(x_{m},y_{n})=(m/N,n/N)$. Each element of $\Omega_{h}$
is of the form $Q_{mn}=[x_{m},x_{m+1}]\times[y_{n},y_{n+1}]$ for some
$m,n=0,\dots,N-1$, see Fig. 1(a) for an illustration.
The above construction has a tensor structure: $\Omega=K\times K$,
$\Omega_{h}=K_{h}\times K_{h}$ and $Q_{mn}=E_{m}\times E_{n}$.

###  2.2 $x$-normal domains

We now introduce the more general class of $x$-normal domains, we discretise
such domains through a curvilinear mesh and we provide some related results.
Consider a smooth Cartesian curve $x=L(y)$ for $y\in[0,1]$ such that $L(y)>0$
for all $y\in[0,1]$. Consider the following _$x$ -normal domain_

$\Omega^{L}:=\\{(x,y)\in\mathbb{R}^{2}|y\in[0,1],\ 0\leq x\leq L(x)\\}.$ (11)

The reference domain $\Omega=[0,1]^{2}$ and the $x$-normal domain $\Omega^{L}$
are linked by the diffeomorphism $\eta:\Omega\rightarrow\Omega^{L}$ defined as
follows

$(x^{L},y^{L})=\eta(x,y)=(xL(y),y),\qquad\forall(x,y)\in\Omega.$ (12)

The idea of mapping a class of domains onto the reference square in order to
exploit the tensor structure of the mesh is reminiscent, for example, of the
work in [15]. Let $\boldsymbol{x}^{L}_{mn}:=\eta(\boldsymbol{x}_{mn})$,
$m,n=0,\dots,N$ be the transformed nodes and let
$\Omega^{L}_{h}:=\eta(\Omega_{h})$ be the transformed (curvilinear) mesh. For
each element $Q\in\Omega_{h}$ let $T:=\eta(Q)$ be the corresponding
transformed element, see Fig. 1(b) for an illustration. The Jacobian
$J_{\eta}$ of $\eta$ is given by

$J_{\eta}(x,y)=\begin{pmatrix}L(y)&xL^{\prime}(y)\\\ 0&1\end{pmatrix},$ (13)

so that

$\det J_{\eta}(x,y)=L(y).$ (14)

The inverse transformation $\eta^{-1}:\Omega^{L}\rightarrow\Omega$ is given by

$(x,y)=\eta^{-1}(x^{L},y^{L})=\left(\frac{x^{L}}{L(y^{L})},y^{L}\right),$ (15)

so the Jacobian $J_{\eta^{-1}}$ of $\eta^{-1}$ is given by

$J_{\eta^{-1}}(x,y)=\begin{pmatrix}\frac{1}{L(y^{L})}&-\frac{x^{L}L^{\prime}(y)}{L^{2}(y)}\\\
0&1\end{pmatrix}.$ (16)

Switching back to the original coordinates, we have

$J_{\eta^{-1}}\circ\eta(x,y)=\begin{pmatrix}\frac{1}{L(y)}&-\frac{xL^{\prime}(y)}{L(y)}\\\
0&1\end{pmatrix},$ (17)

which implies that

$H(x,y):=(J_{\eta^{-1}}\circ\eta(x,y))(J_{\eta^{-1}}\circ\eta(x,y))^{T}=\begin{pmatrix}\frac{1}{L^{2}(y)}+\frac{x^{2}L^{\prime
2}(y)}{L^{2}(y)}&-\frac{xL^{\prime}(y)}{L(y)}\\\
-\frac{xL^{\prime}(y)}{L(y)}&1\end{pmatrix}.$ (18)

We finally define the matrix

$\widehat{H}(x,y):=H(x,y)\det
J_{\eta}(x,y)=\begin{pmatrix}\frac{1}{L(y)}+\frac{x^{2}L^{\prime
2}(y)}{L(y)}&-xL^{\prime}(y)\\\ -xL^{\prime}(y)&L(y)\end{pmatrix},$ (19)

which is symmetric and uniformly positive definite on $\Omega$.

######  Remark 1 (Symmetric $x$-normal domains)

The proposed theory still holds true on _symmetric $x$-normal domains_ of the
form

$\Omega^{S}:=\\{(x,y)\in\mathbb{R}^{2}|y\in[0,1],\ |x|\leq S(y)\\},$ (20)

where $S(y)>0$, $y\in[0,1]$ is a smooth function. In this case, equations
(12)-(19) hold true by setting $\Omega=[-1/2,1/2]\times[0,1]$ and
$L(y)=2S(y)$. A symmetric $x$-normal domain $\Omega^{S}$ of the form (20) with
its discretisation $\Omega_{h}^{S}$ is shown in Fig. 1(c). It is worth
remarking that the proposed theory can be easily extended to _non-symmetric_
$x$-normal domains, we do not consider this case for ease of presentation.

######  Remark 2 (Curvilinear cylindrical surfaces)

Every $x$-normal domain (11) satisfying $L(0)=L(1)$ or symmetric $x$-normal
domain (20) satisfying $S(0)=S(1)$ can be wrapped around a cylinder through
the transformation $\sigma:\Omega^{L}\rightarrow\mathbb{R}^{3}$ defined by

$\sigma(x,y)=\left(x,\frac{\sin 2\pi y}{2\pi},\frac{\cos 2\pi
y}{2\pi}\right),\qquad(x,y)\in\Omega^{L}.$ (21)

We can thus define the curvilinear cylinder $\Gamma$ and its curved mesh
$\Gamma_{h}$ as

$\Gamma:=\sigma(\Omega^{L});\qquad\Gamma_{h}:=\sigma(\Omega^{L}_{h}),$ (22)

respectively, see Fig. 2 for an illustration. We will show that the matrix
approach proposed in the next section can be applied also to special surface
PDEs, see for instance [12].

######  Remark 3 (More general domains)

The choice of normal domains is justified by the property that each entry of
the matrix $\widehat{H}(x,y)$ defined in (19) is a finite sum of separable
terms. The proposed approach still applies if the coordinate transformation
$\eta$ defined in (12) is such that the corresponding $\widehat{H}(x,y)$
retains this property.

$Q_{00}$$Q_{10}$$Q_{01}$$Q_{11}$$\boldsymbol{x}_{00}$$\boldsymbol{x}_{10}$$\boldsymbol{x}_{20}$$\boldsymbol{x}_{01}$$\boldsymbol{x}_{11}$$\boldsymbol{x}_{21}$$\boldsymbol{x}_{02}$$\boldsymbol{x}_{12}$$\boldsymbol{x}_{22}$
(a) Unit square $\Omega$ and its approximation $\Omega_{h}$ for $N=2$.

$T_{00}$$T_{10}$$T_{01}$$T_{11}$$\boldsymbol{x}_{00}^{L}$$\boldsymbol{x}_{10}^{L}$$\boldsymbol{x}_{20}^{L}$$\boldsymbol{x}_{01}^{L}$$\boldsymbol{x}_{11}^{L}$$\boldsymbol{x}_{21}^{L}$$\boldsymbol{x}_{02}^{L}$$\boldsymbol{x}_{12}^{L}$$\boldsymbol{x}_{22}^{L}$
(b) An $x$-normal domain $\Omega^{L}$ with $L(y)=2+\cos 2\pi y$ and its
approximation $\Omega_{h}^{L}$ for $N=2$.

$T_{00}$$T_{10}$$T_{01}$$T_{11}$$\boldsymbol{x}_{00}^{L}$$\boldsymbol{x}_{10}^{L}$$\boldsymbol{x}_{20}^{L}$$\boldsymbol{x}_{01}^{L}$$\boldsymbol{x}_{11}^{L}$$\boldsymbol{x}_{21}^{L}$$\boldsymbol{x}_{02}^{L}$$\boldsymbol{x}_{12}^{L}$$\boldsymbol{x}_{22}^{L}$
(c) A symmetric $x$-normal domain $\Omega^{S}$ with $S(y)=1+\frac{1}{2}\cos
2\pi y$ and its approximation $\Omega_{h}^{S}$ for $N=2$.

Figure 1: Pictorial illustration of the classes of spatial domains considered
in this work, together with their curvilinear meshes.

(a) Symmetric $x$-normal domain $\Omega^{S}$ of the form (20) with
$S(y)=1+\frac{1}{2}\sin 2\pi y$ with its mesh $\Omega^{S}_{h}$.

(b) Curvilinear cylinder $\Gamma$ defined in (22) with its mesh $\Gamma_{h}$
corresponding to the 2D x-normal domain in (a).

Figure 2: An $x$-normal domain $\Omega^{S}$ and its mesh $\Omega_{h}^{S}$ are
transformed into the curvilinear cylinder $\Gamma$ and its mesh $\Gamma_{h}$,
respectively. The transformation $\sigma$ in (21) joins the top- and bottom
edges of $\Omega_{S}$, highlighted in red.

##  3 Matrix-oriented formulation of FEM

We are now ready to (i) formulate a finite element method on the (possibly
curvilinear) meshes introduced in the previous section and (ii) derive a
Kronecker decomposition of the discrete differential operators involved in
problems (1)-(3). To make the reader familiar with the notations and results,
we start with the case of square domains.

###  3.1 Square domains

Let $\\{\psi_{i}\\}_{i=1}^{N+1}$ be any finite element basis functions on the
one-dimensional mesh $I_{h}$, e.g. piecewise Lagrange polynomials of any fixed
degree. The corresponding stiffness- and mass matrices
$A,M\in\mathbb{R}^{(N+1)\times(N+1)}$ in 1D are defined by

$a_{ij}:=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{j}}}{\mathrm{d}x}(x)\mathrm{d}x;\qquad
m_{ij}:=\int_{0}^{1}\psi_{i}(x)\psi_{j}(x)\mathrm{d}x,$ (23)

for $i,j=1,\dots,N+1$, respectively. On $\Omega_{h}$ we choose the tensor-
product local Lagrange basis $\\{\phi_{j}\\}_{j=1}^{(N+1)^{2}}$ defined as

$\phi_{i+(N+1)j}(x,y):=\psi_{i}(x)\psi_{j}(y),\qquad\forall\ i,j=1,\dots,N+1.$
(24)

In general we will write $\phi_{i}(x,y)=\phi_{i_{x}}(x)\phi_{i_{y}}(y)$ for
all $i=1,\dots,(N+1)^{2}$, meaning that every 2D basis function can be
uniquely decomposed as a product of 1D basis functions. The stiffness- and
mass matrices
$\widetilde{A},\widetilde{M}\in\mathbb{R}^{(N+1)^{2}\times(N+1)^{2}}$ in 2D
are defined by

$\widetilde{a}_{ij}:=\iint_{\Omega}\nabla\phi_{i}(x,y)\cdot\nabla\phi_{j}(x,y)\mathrm{d}x\mathrm{d}y;\qquad\widetilde{m}_{ij}:=\iint_{\Omega}\phi_{i}(x,y)\phi_{j}(x,y)\mathrm{d}x\mathrm{d}y,$
(25)

for $i,j=1,\dots,(N+1)^{2}$, respectively. The stiffness- and mass matrices
$\widetilde{A}$ and $\widetilde{M}$ in 2D fulfil the following relations:

$\begin{split}\widetilde{m}_{ij}=&\int_{0}^{1}\int_{0}^{1}\psi_{i_{x}}(x)\psi_{i_{y}}(y)\psi_{j_{x}}(x)\psi_{j_{y}}(y)\mathrm{d}x\mathrm{d}y\\\
=&\int_{0}^{1}\psi_{i_{x}}(x)\psi_{j_{x}}(x)\mathrm{d}x\int_{0}^{1}\psi_{i_{y}}(y)\psi_{j_{y}}(y)\mathrm{d}y=m_{i_{x}j_{x}}m_{i_{y}j_{y}};\end{split}$
(26)
$\begin{split}\widetilde{a}_{ij}&=\int_{0}^{1}\int_{0}^{1}\nabla\left(\psi_{i_{x}}(x)\psi_{i_{y}}(y)\right)\cdot\nabla\left(\psi_{j_{x}}(x)\psi_{j_{y}}(y)\right)\mathrm{d}x\mathrm{d}y\\\
&=\int_{0}^{1}\int_{0}^{1}\frac{\partial}{\partial
x}\left(\psi_{i_{x}}(x)\psi_{i_{y}}(y)\right)\frac{\partial}{\partial
x}\left(\psi_{j_{x}}(x)\psi_{j_{y}}(y)\right)\mathrm{d}x\mathrm{d}y\\\
&+\int_{0}^{1}\int_{0}^{1}\frac{\partial}{\partial
y}\left(\psi_{i_{x}}(x)\psi_{i_{y}}(y)\right)\frac{\partial}{\partial
y}\left(\psi_{j_{x}}(x)\psi_{j_{y}}(y)\right)\mathrm{d}x\mathrm{d}y\\\
&=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i_{x}}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{j_{x}}}}{\mathrm{d}x}(x)\mathrm{d}x\int_{0}^{1}\psi_{i_{y}}(y)\psi_{j_{y}}(y)\mathrm{d}y\\\
&+\int_{0}^{1}\psi_{i_{x}}(x)\psi_{j_{x}}(x)\mathrm{d}x\int_{0}^{1}\frac{\mathrm{d}{\psi_{i_{y}}}}{\mathrm{d}y}(y)\frac{\mathrm{d}{\psi_{j_{y}}}}{\mathrm{d}y}(y)\mathrm{d}y=a_{i_{x},j_{x}}m_{i_{y}j_{y}}+m_{i_{x},j_{x}}a_{i_{y}j_{y}}.\end{split}$
(27)

Equations (26)-(27) translate to the following matrix identities:

$\displaystyle\widetilde{M}=M\otimes M;\qquad\widetilde{A}=A\otimes M+M\otimes
A,$ (28)

where $\otimes$ denotes the Kronecker product.

###  3.2 $x$-normal domains

Let us now consider the more general case of $x$-normal and symmetric
$x$-normal domains. We introduce the novel 1D matrices
$B_{1},B_{2},C_{1},C_{2},M_{1},M_{2},M_{3}\in\mathbb{R}^{(N+1)\times(N+1)}$,
defined as follows:

$\displaystyle
b_{ij}^{1}:=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{j}}}{\mathrm{d}x}(x)L(x)\mathrm{d}x;\qquad
b_{ij}^{2}:=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{j}}}{\mathrm{d}x}(x)x^{2}\mathrm{d}x;$
(29) $\displaystyle
c_{ij}^{1}:=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\psi_{j}(x)x\mathrm{d}x;\qquad\hskip
22.76219ptc_{ij}^{2}:=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\psi_{j}(x)L^{\prime}(x)\mathrm{d}x;$
(30)
$\displaystyle\begin{split}&m_{ij}^{1}:=\int_{0}^{1}\psi_{i}(x)\psi_{j}(x)\frac{1}{L(x)}\mathrm{d}x;\qquad\hskip
8.53581ptm_{ij}^{2}:=\int_{0}^{1}\psi_{i}(x)\psi_{j}(x)\frac{L^{\prime
2}(x)}{L(x)}\mathrm{d}x;\\\
&m_{ij}^{3}:=\int_{0}^{1}\psi_{i}(x)\psi_{j}(x)L(x)\mathrm{d}x,\end{split}$
(31)

for $i,j=1,\dots,(N+1)$. On $\Omega_{h}^{L}$ we choose the transformed
Lagrange basis $\\{\phi_{i}^{L}\\}_{i=1}^{(N+1)^{2}}$ where
$\phi_{i}^{L}(x,y):=\phi_{i}(\eta^{-1}(x,y))$. The stiffness- and mass
matrices $\widetilde{A},\widetilde{M}\in\mathbb{R}^{(N+1)^{2}\times(N+1)^{2}}$
in 2D are defined by

$\widetilde{a}_{ij}:=\iint_{\Omega^{L}}\nabla\phi_{i}^{L}(x,y)\cdot\nabla\phi_{j}^{L}(x,y)\mathrm{d}x\mathrm{d}y;\qquad\widetilde{m}_{ij}:=\iint_{\Omega^{L}}\phi_{i}^{L}(x,y)\phi_{j}^{L}(x,y)\mathrm{d}x\mathrm{d}y,$
(32)

for $i,j=1,\dots,(N+1)^{2}$, respectively. By using (19), the stiffness and
mass matrices $\widetilde{A}$ and $\widetilde{M}$ in 2D fulfil the following
relations:

$\begin{split}\widetilde{m}_{ij}&=\iint_{\Omega^{L}}\phi_{i}^{L}(x,y)\phi_{j}^{L}(x,y)\mathrm{d}x\mathrm{d}y=\iint_{\Omega}\phi_{i}(x,y)\phi_{j}(x,y)L(y)\mathrm{d}x\mathrm{d}y\\\
&=\int_{0}^{1}\psi_{i_{x}}(x)\psi_{j_{x}}(x)\mathrm{d}x\int_{0}^{1}\psi_{i_{y}}(y)\psi_{j_{y}}(y)L(y)\mathrm{d}y=m_{i_{x}j_{x}}m_{i_{y}j_{y}}^{3};\end{split}$
(33)
$\begin{split}\widetilde{a}_{ij}&=\iint_{\Omega^{L}}\nabla\phi_{i}^{L}(x,y)\cdot\nabla\phi_{j}^{L}(x,y)\mathrm{d}x\mathrm{d}y\\\
&=\iint_{\Omega^{L}}\left(\nabla\phi_{i}(\eta^{-1}(x,y))J_{\eta^{-1}}(x,y)\right)\cdot\left(\nabla\phi_{j}(\eta^{-1}(x,y))J_{\eta^{-1}}(x,y)\right)\mathrm{d}x\mathrm{d}y\\\
&=\iint_{\Omega}\nabla\phi_{i}(x,y)^{T}H(x,y)\nabla\phi_{j}(x,y)\det
J_{\eta}(x,y)\mathrm{d}x\mathrm{d}y\\\
&=\int_{0}^{1}\frac{\mathrm{d}{\psi_{ix}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{jx}}}{\mathrm{d}x}(x)\mathrm{d}x\int_{0}^{1}\psi_{i_{y}}(y)\psi_{j_{y}}(y)\frac{1}{L(y)}\mathrm{d}y\\\
&+\int_{0}^{1}\psi_{ix}(x)\psi_{jx}(x)\mathrm{d}x\int_{0}^{1}\frac{\mathrm{d}{\psi_{i_{y}}}}{\mathrm{d}y}(y)\frac{\mathrm{d}{\psi_{j_{y}}}}{\mathrm{d}y}(y)L(y)\mathrm{d}y\\\
&+\int_{0}^{1}\frac{\mathrm{d}{\psi_{ix}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{jx}}}{\mathrm{d}x}(x)x^{2}\mathrm{d}x\int_{0}^{1}\psi_{i_{y}}(y)\psi_{j_{y}}(y)\frac{L^{\prime
2}(y)}{L(y)}\mathrm{d}y\\\
&-\int_{0}^{1}\frac{\mathrm{d}{\psi_{ix}}}{\mathrm{d}x}(x)\psi_{jx}(x)x\mathrm{d}x\int_{0}^{1}\psi_{i_{y}}(y)\frac{\mathrm{d}{\psi_{j_{y}}}}{\mathrm{d}y}(y)L^{\prime}(y)\mathrm{d}y\\\
&-\int_{0}^{1}\psi_{ix}(x)\frac{\mathrm{d}{\psi_{jx}}}{\mathrm{d}x}(x)x\mathrm{d}x\int_{0}^{1}\frac{\mathrm{d}{\psi_{i_{y}}}}{\mathrm{d}y}(y)\psi_{j_{y}}(y)L^{\prime}(y)\mathrm{d}y\\\
&=a_{i_{x}j_{x}}m_{i_{y}j_{y}}^{1}+m_{i_{x}j_{x}}b_{i_{y}j_{y}}^{1}+b_{i_{x}j_{x}}^{2}m_{i_{y}j_{y}}^{2}-c_{i_{x}j_{x}}^{1}c_{j_{y}i_{y}}^{2}-c_{j_{x}i_{x}}^{1}c_{i_{y}j_{y}}^{2}.\end{split}$
(34)

Equations (33)-(34) translate to the following matrix identities:

$\displaystyle\widetilde{M}=M\otimes M_{3};\qquad\widetilde{A}=A\otimes
M_{1}+M\otimes B_{1}+B_{2}\otimes M_{2}-C_{1}\otimes
C_{2}^{T}-C_{1}^{T}\otimes C_{2},$ (35)

where $A$ and $M$ are the standard stiffness- and mass matrices in 1D defined
in (23) and $B_{1},B_{2},C_{1},C_{2}$, $M_{1},M_{2},M_{3},$ are novel matrices
defined in (29)-(31). Observe that on square domains, i.e. when $L(y)=1$, we
have $M_{1}=M_{3}=M$, $B_{1}=A$, $M_{2}=C_{2}=0$, so the matrix identities
(35) encompass the special case (28) of square domains.

###  3.3 $x$-normal domains: $\mathbb{P}_{1}$ elements with lumping

The Kronecker decompositions (35) are absolutely general and encompass
arbitrary choices of the FEM basis function, thereby including the case of
Lagrangian spatial methods of any order. However, for the sake of
practicality, we consider the special case of Lagrangian $\mathbb{P}_{1}$
finite elements _with mass and stiffness lumping_ , see [24]. We show that
such special case significantly simplifies the matrix identities (35). This
has two advantages: (i) existing finite element codes can be adapted to the
proposed approach with minor modifications and (ii) the resulting numerical
schemes, which take the form of matrix equations, become much easier to solve.
In the remainder of this section, we specialise $\\{\psi_{i}\\}_{i=1}^{N+1}$
to be the standard $\mathbb{P}_{1}$ Lagrangian (also called _pyramidal_) basis
functions. We consider the (diagonal) lumped mass matrix
$M_{0}\in\mathbb{R}^{(N+1)\times(N+1)}$ in 1D defined by

$m_{ij}^{0}:=\int_{0}^{1}I_{h}(\psi_{i}(x)\psi_{j}(x))\mathrm{d}x=\delta_{ij}\int_{0}^{1}\psi_{i}(x)\mathrm{d}x,\qquad
i,j=1,\dots,N+1,$ (36)

with $I_{h}$ being the element-wise interpolant operator [24] and
$\delta_{ij}$ being the Kronecker symbol. We also consider the (tri-diagonal)
convection matrix $C\in\mathbb{R}^{(N+1)\times(N+1)}$ defined by

$c_{ij}:=\int_{0}^{1}\frac{\mathrm{d}{\psi_{i}(x)}}{\mathrm{d}x}\psi_{j}(x)\mathrm{d}x,\qquad
i,j=1,\dots,N+1.$ (37)

Then we consider the (tri-diagonal) _modified stiffness matrices_
$A_{1},A_{2}\in\mathbb{R}^{(N+1)\times(N+1)}$ defined by

$\displaystyle
a_{ij}^{1}:=\int_{0}^{1}I_{h}\left(\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{j}}}{\mathrm{d}x}(x)L(x)\right)\mathrm{d}x,\qquad
a_{ij}^{2}:=\int_{0}^{1}I_{h}\left(\frac{\mathrm{d}{\psi_{i}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{j}}}{\mathrm{d}x}(x)x^{2}\right)\mathrm{d}x,$
(38)

for $i,j=1,\dots,N+1$. We finally define the following auxiliary diagonal
matrices $D_{1},D_{2},D_{3}\in\mathbb{R}^{(N+1)\times(N+1)}$ defined as
follows

$d_{ij}^{1}=\delta_{ij}L(x_{i}),\qquad
d_{ij}^{2}=\delta_{ij}L^{\prime}(x_{i}),\qquad
d_{ij}^{3}=\delta_{ij}x_{i},\qquad i,j=1,\dots,N+1.$ (39)

By combining (19) and (34) we get

$\widetilde{a}_{ij}=\iint_{\Omega}\left(\widehat{H}(x,y)\nabla\phi_{i}(x,y)\right)\cdot\nabla\phi_{j}(x,y)\mathrm{d}x\mathrm{d}y.$
(40)

The matrices defined by (33) and (40) can be understood as the mass and
stiffness matrices of the anisotropic elliptic equation

$-\nabla(\widehat{H}(x,y)\nabla
u(x,y))+L(y)u(x,y)=f(x,y),\qquad(x,y)\in\Omega,$ (41)

with $\widehat{H}$ as defined in (19). Hence, by following [24], the lumped
counterparts $\widehat{M}$ of $\widetilde{M}$ and $\widehat{A}$ of
$\widetilde{A}$ are defined by

$\displaystyle\begin{split}\widehat{m}_{ij}&=\iint_{\Omega}I_{h}\left(\phi_{i}(x,y)\phi_{j}(x,y)L(y)\right)\mathrm{d}x\mathrm{d}y;\\\
\widehat{a}_{ij}&=\iint_{\Omega}I_{h}\left(\left(\widehat{H}(x,y)\nabla\phi_{i}(x,y)\right)\cdot\nabla\phi_{j}(x,y)\right)\mathrm{d}x\mathrm{d}y,\end{split}$
(42)

for $i,j=1,\dots,(N+1)^{2}$, where $I_{h}$ is the element-wise, vector-valued
interpolant operator, see [12]. So, the tensorial decompositions of (42) are
carried out as follows

$\begin{split}\widehat{m}_{ij}&=\int_{0}^{1}I_{h}(\psi_{i_{x}}(x)\psi_{j_{x}}(x))\mathrm{d}x\int_{0}^{1}I_{h}(\psi_{i_{y}}(y)\psi_{j_{y}}(y)L(y))\mathrm{d}y=m^{0}_{i_{x}j_{x}}d^{1}_{i_{y}j_{y}}m^{0}_{i_{y}j_{y}};\end{split}$
(43)
$\begin{split}\widehat{a}_{ij}&=\int_{0}^{1}\frac{\mathrm{d}{\psi_{ix}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{jx}}}{\mathrm{d}x}(x)\mathrm{d}x\int_{0}^{1}I_{h}\left(\psi_{i_{y}}(y)\psi_{j_{y}}(y)\frac{1}{L(y)}\right)\mathrm{d}y\\\
&+\int_{0}^{1}I_{h}\left(\psi_{ix}(x)\psi_{jx}(x)\right)\mathrm{d}x\int_{0}^{1}I_{h}\left(\frac{\mathrm{d}{\psi_{i_{y}}}}{\mathrm{d}y}(y)\frac{\mathrm{d}{\psi_{j_{y}}}}{\mathrm{d}y}(y)L(y)\right)\mathrm{d}y\\\
&+\int_{0}^{1}I_{h}\left(\frac{\mathrm{d}{\psi_{ix}}}{\mathrm{d}x}(x)\frac{\mathrm{d}{\psi_{jx}}}{\mathrm{d}x}(x)x^{2}\right)\mathrm{d}x\int_{0}^{1}I_{h}\left(\psi_{i_{y}}(y)\psi_{j_{y}}(y)\frac{L^{\prime
2}(y)}{L(y)}\right)\mathrm{d}y\\\
&-\int_{0}^{1}I_{h}\left(\frac{\mathrm{d}{\psi_{ix}}}{\mathrm{d}x}(x)\psi_{jx}(x)x\right)\mathrm{d}x\int_{0}^{1}I_{h}\left(\psi_{i_{y}}(y)\frac{\mathrm{d}{\psi_{j_{y}}}}{\mathrm{d}y}(y)L^{\prime}(y)\right)\mathrm{d}y\\\
&-\int_{0}^{1}I_{h}\left(\psi_{ix}(x)\frac{\mathrm{d}{\psi_{jx}}}{\mathrm{d}x}(x)x\right)\mathrm{d}x\int_{0}^{1}I_{h}\left(\frac{\mathrm{d}{\psi_{i_{y}}}}{\mathrm{d}y}(y)\psi_{j_{y}}(y)L^{\prime}(y)\right)\mathrm{d}y\\\
&=\frac{a_{i_{x}j_{x}}m_{i_{y}j_{y}}^{0}}{d^{1}_{i_{y}i_{y}}}+m^{0}_{i_{x}j_{x}}a_{i_{y}j_{y}}^{1}+\frac{a_{i_{x}j_{x}}^{2}m_{i_{y}j_{y}}^{0}(d^{2}_{i_{y}i_{y}})^{2}}{d^{1}_{i_{y}i_{y}}}-c_{i_{x}j_{x}}d^{3}_{j_{x}j_{x}}c_{j_{y}i_{y}}d^{2}_{i_{y}i_{y}}-c_{j_{x}i_{x}}d^{3}_{i_{x}i_{x}}c_{i_{y}j_{y}}d^{2}_{j_{y}j_{y}}.\end{split}$
(44)

Equations (43)-(44) translate to the following matrix identities:

$\displaystyle\begin{split}\widehat{M}&=M_{0}\otimes M_{0}D_{1};\\\
\widehat{A}&=A\otimes D_{1}^{-1}M_{0}+M_{0}\otimes A_{1}+A_{2}\otimes
D_{2}^{2}D_{1}^{-1}M_{0}-CD_{3}\otimes D_{2}C^{T}-D_{3}C^{T}\otimes
CD_{2}.\end{split}$ (45)

As mentioned earlier, the matrix relations (45) are simpler than (35). This is
because (i) the matrices $D_{1},D_{2},D_{3},M_{0}$ are now diagonal and (ii)
the only non-diagonal matrices involved, i.e. $A,A_{1},A_{2},C$ are now
tridiagonal.

######  Remark 4 (Homogeneous Dirichlet boundary conditions)

In the presence of homogeneous Dirichlet boundary conditions, it is well-known
that the boundary basis functions must be eliminated from the basis, see [16],
hence all the matrices involved in (45) must be trimmed by eliminating all
boundary entries. Hence, the dimension of such matrices drops from
$(N+1)\times(N+1)$ to $(N-1)\times(N-1)$ and the following additional
properties hold true.

  * •

$M_{0}$ is a multiple of the identity: for $m_{0}:=\frac{1}{N}$ it holds that

$M_{0}=m_{0}I.$ (46)

  * •

The stiffness matrix $A$ is Toeplitz.

  * •

Thanks to the symmetries of the $\psi_{i}$’s, $C$ is Toeplitz and skew-
symmetric, i.e.

$C^{T}=-C.$ (47)

Because the dimension of the matrices depends on the kind of boundary
conditions, see Remark 4, we set

$q:=\begin{cases}N+1\qquad\text{for Neumann boundary conditions}\\\
N-1\qquad\text{for Dirichlet boundary conditions}.\end{cases}$ (48)

In the presence of homogeneous Dirichlet boundary conditions, thanks to (46)
and (47), relations (45) become

$\displaystyle\begin{split}&\widehat{M}=m_{0}^{2}I\otimes D_{1};\\\
&\widehat{A}=m_{0}(A\otimes D_{1}^{-1}+I\otimes A_{1}+A_{2}\otimes
D_{2}^{2}D_{1}^{-1})+CD_{3}\otimes D_{2}C+D_{3}C\otimes CD_{2}.\end{split}$
(49)

######  Remark 5 (Recap on lumped $\mathbb{P}_{1}$ matrix properties)

We recap here the properties of the matrices appearing in (45)and (49):

  * •

$M_{0}$ and $D_{1}$ are diagonal and positive definite. In the Dirichlet case
$M_{0}$ is multiple of the identity;

  * •

$D_{2}$ is diagonal and it is non-singular only when the curve $x=L(y)$ is
strictly monotone;

  * •

$D_{3}$ is diagonal and is singular when $N$ is odd, or tends to being
singular when $N$ is even and approaches infinity;

  * •

$A,A_{1}$ are tridiagonal, moreover they are positive definite in the
Dirichlet case and semidefinite in the Neumann case (one null eigenvalue). $A$
is symmetric and in the Dirichlet case it is Toeplitz;

  * •

$A_{2}$ is tridiagonal non-symmetric, it is singular for $N$ odd, or tends to
being singular for $N$ even approaching infinity;

  * •

$C=\texttt{tridiag}([-1,0,1])$, except for boundary entries in the Neumann
case (it is the matrix of centered first derivatives). According to $N$ being
even or odd, and depending on the boundary conditions, $C$ is singular or
tends to being singular as $N\rightarrow+\infty$.

##  4 Stationary PDEs and the Sylvester equation

In this section we consider the Poisson equation on the unit square
$\Omega=[0,1]^{2}$:

$-\Delta u(\boldsymbol{x})+\gamma
u(\boldsymbol{x})=f(\boldsymbol{x}),\qquad\boldsymbol{x}\in\Omega,$ (50)

where $\gamma\geq 0$ for the case of Dirichlet boundary conditions and
$\gamma>0$ for the case of Neumann boundary conditions. Since we encompass
general boundary condition, we use the notation in (48) for the dimension of
matrices. The general FEM discretisation of problem (50) in vector form is
then

$\widetilde{A}\boldsymbol{\xi}+\gamma\widetilde{M}\boldsymbol{\xi}=\widetilde{M}\boldsymbol{f},$
(51)

with $\boldsymbol{\xi}$ being the nodal vector of the numerical solution and
$\boldsymbol{f}$ being the corresponding nodal vector of $f$. By using (28),
the linear system (51) becomes the following generalised Sylvester matrix
equation:

$\left(A+\frac{\gamma}{2}M\right)UM+MU\left(A+\frac{\gamma}{2}M\right)=MFM,$
(52)

where $U$ and $F$ are such that $\texttt{vec}(U)=\boldsymbol{\xi}$ and
$\texttt{vec}(F)=\boldsymbol{f}$. Since the mass matrix $M$ is positive-
definite, we can pre- and post-multiply both sides of (52) by $M^{-1}$:

$\left(M^{-1}A+\frac{\gamma}{2}I\right)U+U\left(AM^{-1}+\frac{\gamma}{2}I\right)=F.$
(53)

Even if $M^{-1}A$ and $AM^{-1}$ are not symmetric, they are diagonalizable
nonetheless because $A,M$ are both symmetric and $M$ is positive definite.
Hence, also $Z_{1}:=M^{-1}A+\frac{\gamma}{2}I$ and
$Z_{2}:=AM^{-1}+\frac{\gamma}{2}I$ are diagonalizable. Then, (54) becomes

$Z_{1}U+UZ_{2}=F,$ (54)

and we can diagonalise $Z_{1}$, $Z_{2}$ as follows

$\displaystyle Z_{1}=X^{(1)}\Lambda^{(1)}{X^{(1)}}^{-1};\qquad
Z_{2}=X^{(2)}\Lambda^{(2)}{X^{(2)}}^{-1},$ (55)

with $\Lambda^{(k)}\in\mathbb{R}^{q\times q}$, $k=1,2$, being the diagonal
matrices containing the eigenvalues $\lambda_{i}^{(k)}$, $i=1,\dots,q$ of
$Z_{1}$ and $Z_{2}$, respectively. Except special cases, the $\Lambda^{(k)}$’s
and $X^{(k)}$’s must be computed numerically. Now, by setting
$\widehat{U}:={X^{(1)}}^{-1}UX^{(2)}$ and
$\widehat{F}:={X^{(1)}}^{-1}FX^{(2)}$, (54) becomes

$\Lambda^{(1)}\widehat{U}+\widehat{U}\Lambda^{(2)}=\widehat{F},$ (56)

which can be solved as follows. Let $L\in\mathbb{R}^{q\times q}$ be the matrix
defined by

$\ell_{ij}=\frac{1}{\lambda_{i}+\lambda_{j}},\qquad i,j=1,\dots,q.$ (57)

As shown in [9], the solution to (56) in the spectral space can be expressed
as:

$\widehat{U}=L\circ\widehat{F},$ (58)

with $\circ$ denoting the Hadamard product. The original variable $U$ is thus
given by

$U=X^{(1)}(L\circ\widehat{F}){X^{(2)}}^{-1}.$ (59)

We call this technique _FEM reduced method_.

######  Remark 6 (Memory performance of the reduced approach)

Because $U$, $F$, $L$, $X^{(1)}$, $X^{(2)}$ are the only full matrices
involved in the computation, the memory occupation of the proposed approach is
$5N^{2}+O(kN)$ floating point numbers for any $k$. On the other hand, since
the large matrices $\widetilde{A}$ and $\widetilde{M}$ are
$(2k+1)^{2}$-diagonal, where $k$ is the polynomial order of the method, and
$\boldsymbol{f}$ and $\boldsymbol{\xi}$ are both full vectors, the vector
formulation (51) has a memory occupation of $((2k+1)^{2}+2)N^{2}+O(kN)$
floating point numbers.

###  4.1 Special case: Lagrangian $\mathbb{P}_{1}$ elements with Dirichlet
boundary conditions

In the special case of Lagrangian $\mathbb{P}_{1}$ elements with Dirichlet
boundary conditions, we are able to solve the matrix equation (54) in closed
form without computing spectral decompositions numerically. In fact, since $M$
and $A$ are both symmetric, positive definite, tridiagonal Toeplitz matrices,
there exist $\alpha,\beta>0$ such that $M=\alpha A+\beta I$. In this specific
case, we have

$\alpha=-\frac{1}{6N^{2}},\qquad\beta=\frac{1}{N}.$ (60)

This implies that (i) $M^{-1}$ and $A$ commute, (ii) $M^{-1}$ and $A$ share
the same eigenvectors and (iii) $M^{-1}A=AM^{-1}$ is diagonalizable and shares
the same eigenvectors of $M^{-1}$ and $A$. We can thus write

$\displaystyle Z_{1}=Z_{2}=X\Lambda X^{-1},$ (61)

with $\Lambda\in\mathbb{R}^{(N-1)\times(N-1)}$ being the diagonal matrix
containing the eigenvalues $\lambda_{i}$, $i=1,\dots,N-1$ of $Z_{1}=Z_{2}$.
Now, by setting $\widehat{U}:=X^{-1}UX$ and $\widehat{F}:=X^{-1}FX$, (54)
becomes

$\Lambda\widehat{U}+\widehat{U}\Lambda=\widehat{F},$ (62)

in the spectral space. We are left to show that $X$ and $\Lambda$ are known in
closed form. In fact, as shown in [8], the eigenvalues and eigenvectors of $A$
are given by

$\displaystyle\lambda^{A}_{i}=\frac{2}{N}\left(1-\cos\frac{i\pi}{N}\right),\qquad
i=1,\dots,N-1;$ (63)
$\displaystyle(\boldsymbol{v}^{A}_{i})_{j}=\sin\frac{ij\pi}{N},\qquad
i,j=1,\dots,N-1.$ (64)

Hence, by using (60) the entries of $\Lambda$ are given by

$\lambda_{i}=\frac{(12N^{2}-\gamma)\lambda^{A}_{i}+6N\gamma}{12N-2\lambda_{i}^{A}},\quad
i=1,\dots,N-1,$ (65)

while the entries of $X$ (common basis of eigenvectors of $A$, $Z_{1}$ and
$Z_{2}$) are given in (64). Hence, $U$ is given by (59) with
$\Lambda^{(1)}=\Lambda^{(2)}=\Lambda$ defined by (65) and $X^{(1)}=X^{(2)}=X$
defined by (64).

###  4.2 Numerical example

We consider the following Poisson equation with Dirichlet boundary conditions
on the square $\Omega=[0,1]^{2}$:

$\begin{cases}-&\Delta u(x,y)=8\pi^{2}\sin(2\pi x)\sin(2\pi
y),\qquad(x,y)\in\Omega;\\\
&u(x,y)=0,\qquad(x,y)\in\partial\Omega,\end{cases}$ (66)

whose exact solution is $u(x,y)=\sin(2\pi x)\sin(2\pi y)$. We apply Lagrangian
$\mathbb{P}_{k}$ elements, $k=1,\dots,4$, and we compare the performances of
the Kronecker (vector) formulation (51), solved through the direct solver
mldivide of MATLAB (known also as "backslash" $\backslash$ ), and the
_reduced_ approach (56). For $k=1$, we further compare the aforementioned
methods with the _reduced_ method in closed form (62). Here we consider a
sequence of $7$ meshes $\Gamma_{i}$ with $N_{i}=24\cdot 2^{i}$ for all
$i=0,\dots,6$. Such $N_{i}$’s are compatible with $\mathbb{P}_{k}$ elements
for all $k=1,\dots,4$. The numerical results are shown in Fig. 3. On the
finest mesh ($N=1536$), the reduced approach (56) for $k=1$ is approximately
$1.38$ times quicker than the direct (vector) method, $1.58$ times quicker for
$k=2$, $1.70$ times quicker for $k=3$, and $1.96$ times quicker for $k=4$. The
_reduced_ approach in closed form (62) (only $k=1$) is $21.73$ times quicker
than the vector formulation. Furthermore, on equal meshes, the methods produce
the same solutions up to rounding errors (they are equivalent to each other)
and exhibit optimal convergence in space, with the case $k=2$ being
superconvergent (fourth order), as we can see in Fig. 3, right plot.
This and all the following experiments are carried out in MATLAB R2019a on a
HP Z230 Tower Workstation with Intel Core i7-440 CPU and 16GB RAM. The timings
were just taken once and not averaged with several measurements.

Figure 3: Poisson equation (66) on the square $\Omega=[0,1]^{2}$. FEM of order
$k=1,2,3,4$: comparison between the Kronecker (vector) approach (51) (dashed
lines), the MO-FEM solved by the _reduced_ approach (56) (continuous lines)
and, only for $k=1$, the reduced approach in closed form (62) (thick
continuous line). Left plot: computational times. Right plot: convergence
behavior. Dotted lines indicate slopes 2,4,5.

##  5 Stationary PDEs on $x$-normal domains and multiterm Sylvester equations

We consider the following stationary PDE problem on an $x$-normal domain
$\Omega^{L}$:

$-\Delta u(\boldsymbol{x})+\gamma
u(\boldsymbol{x})=f(\boldsymbol{x}),\qquad\boldsymbol{x}\in\Omega^{L},$ (67)

with $\gamma\geq 0$ in the presence of zero Dirichlet boundary conditions and
$\gamma>0$ in the presence of zero Neumann boundary conditions. The FEM
discretisation of problem (67) in vector form is then

$\widetilde{A}\boldsymbol{\xi}+\gamma\widetilde{M}\boldsymbol{\xi}=\widetilde{M}\boldsymbol{f},$
(68)

with $\widetilde{A}$ and $\widetilde{M}$ as defined in (32). The lumped
counterpart of (68) is

$\widehat{A}\boldsymbol{\xi}+\gamma\widehat{M}\boldsymbol{\xi}=\widehat{M}\boldsymbol{f},$
(69)

with $\widehat{M}$ and $\widehat{A}$ in (42). Here,
$\boldsymbol{\xi}\in\mathbb{R}^{q^{2}}$ is the nodal vector of the numerical
solution and $\boldsymbol{f}\in\mathbb{R}^{q^{2}}$ being the corresponding
nodal vector of $f$. By using (35), the linear system (68) becomes the
following multiterm Sylvester matrix equation:

$M_{1}UA+(B_{1}+\gamma
M_{3})UM+M_{2}UB_{2}-C_{2}UC_{1}-C_{2}^{T}UC_{1}^{T}=M_{3}FM.$ (70)

Similarly, by using (49) the “lumped” linear system (69) translates to:

$\begin{split}M_{0}D_{1}^{-1}UA+(A_{1}+\gamma
M_{0}D_{1})UM_{0}+M_{0}D_{1}^{-1}D_{2}^{2}UA_{2}-CD_{2}UCD_{3}&-D_{2}C^{T}UD_{3}C^{T}\\\
&=M_{0}D_{1}FM_{0}.\end{split}$ (71)

Thanks to (47), in the case of zero Dirichlet boundary conditions, (71) can be
simplified as:

$\begin{split}m_{0}(D_{1}^{-1}UA+(A_{1}+\gamma
D_{1})U+D_{2}^{2}D_{1}^{-1}UA_{2})-(D_{2}CUD_{3}C+CD_{2}UCD_{3})=m_{0}^{2}D_{1}F.\end{split}$
(72)

We now propose an iterative strategy that can be applied for the solution of
the multiterm matrix equations (70), (71) and (72). For the systems (68) and
(69), the well known Preconditioned Conjugate Gradient method (PCG) [14] would
be a suitable choice since the matrices $\widetilde{A}+\gamma\widetilde{M}$
and $\widehat{A}+\gamma\widehat{M}$ are symmetric and positive definite. Here,
we propose its matrix oriented version to solve the corresponding Sylvester
multiterm equations (70) and (72). To this end, we define the following matrix
operators:

$\displaystyle\widetilde{\mathcal{L}}(U):=M_{1}UA+(B_{1}+\gamma
M_{3})UM+M_{2}UB_{2}-C_{2}UC_{1}-C_{2}^{T}UC_{1}^{T};$ (73)
$\displaystyle\widetilde{\mathcal{R}}(U):=M_{3}FM;$ (74)
$\displaystyle\widehat{\mathcal{L}}(U):=M_{0}D_{1}^{-1}UA+(A_{1}+\gamma
M_{0}D_{1})UM_{0}+M_{0}D_{1}^{-1}D_{2}^{2}UA_{2}-CD_{2}UCD_{3}-D_{2}C^{T}UD_{3}C^{T};$
(75) $\displaystyle\widehat{\mathcal{R}}(U):=M_{0}D_{1}FM_{0};$ (76)

for all $U\in\mathbb{R}^{q\times q}$. We also consider suitable
preconditioning operators
$\widetilde{\mathcal{P}},\widehat{\mathcal{P}}:\mathbb{R}^{q\times
q}\rightarrow\mathbb{R}^{q\times q}$ for the FEM and lumped FEM, respectively,
whose choice will be discussed later. With these settings, the matrix-oriented
formulation of the PCG for systems (68) and (69), that we define as MO-PCG, is
given by

$\displaystyle\begin{cases}R^{(0)}=\mathcal{R}(U^{(0)})-\mathcal{L}(U^{(0)});\\\
Z^{(0)}=\mathcal{P}^{-1}(R^{(0)});\\\ Q^{(0)}=Z^{(0)};\end{cases}$ (77)
$\displaystyle\begin{cases}\alpha^{(s)}=\dfrac{\texttt{sum}(\texttt{sum}(Q^{(s)}\circ
R^{(s)}))}{\texttt{sum}(\texttt{sum}(\mathcal{L}(Q^{(s)})\circ R^{(s)}))};\\\
U^{(s+1)}=U^{(s)}+\alpha^{(s)}Q^{(s)};\\\
R^{(s+1)}=R^{(s)}-\alpha^{(s)}\mathcal{L}(Q^{(s)});\\\
Z^{(s+1)}=\mathcal{P}^{-1}(R^{(s+1)});\\\
\beta^{(s)}=\dfrac{\texttt{sum}(\texttt{sum}(\mathcal{L}(Q^{(s)})\circ
Z^{(s+1)}))}{\texttt{sum}(\texttt{sum}(\mathcal{L}(Q^{(s)})\circ
R^{(s)}))};\\\ Q^{(s+1)}=Q^{(s)}-\beta^{(s)}Q^{(s)},\end{cases}\qquad s\geq
0,$ (78)

where the $\ \widetilde{}\ $ and $\ \widehat{}\ $ are omitted for ease of
presentation. For the initial guess $U^{(0)}$ we can choose for instance
$U^{(0)}=\boldsymbol{0}$.
For illustrative purposes, we consider the following stopping criterion. If
$err_{v}$ is the absolute error obtained by the direct method solving the
linear systems (vector formulation) (68) or (69), we stop the MO-PCG
iterations when the increment fulfils

$\|U^{(s+1)}-U^{(s)}\|_{F}\leq 0.05\|err_{v}\|_{F}$

where $\|\cdot\|_{F}$ denotes the Frobenius norm.

#### Memory performance of the PCG approach

Because $U^{(s)}$, $R^{(s)}$, $Z^{(s)}$, $Q^{(s)}$, $F$ are the only full
matrices involved in the computation, the memory occupation of the proposed
approach is $5N^{2}+O(kN)$, with $N$ being the grid size and $k$ is the FEM
polynomial order On the other hand, the vector - Kronecker formulation (51)
has a memory occupation of $((2k+1)^{2}+2)N^{2}+O(kN)$ for matrix storage,
where $k$ is the polynomial order of the method, as discussed in Remark 6.

#### Matrix PCG vs. classical PCG

The large linear system (68) or (69) could be solved via PCG in classical
vector formulation. However, in the experiments (carried out in MATLAB R2019a
on a HP Z230 Tower Workstation with Intel Core i7-4770 CPU and 16 GB RAM), we
find that the matrix-oriented PCG is significantly faster than its classical
vector counterpart, even if the methods are equivalent. Some details are
provided in the test of Section 5.2.

#### Choice of the preconditioners

The choice of fast and efficient matrix-oriented preconditioners for problems
(70)-(71) is an open problem. Here we consider the case of Dirichlet boundary
conditions and we derive preconditioners that experimentally prove more
efficient than the identity operator, i.e. no preconditioner. To this end, we
first construct suitable preconditioners for the large linear systems
(68)-(69), which we will use to derive matrix-oriented preconditioners for the
multiterm Sylvester problems (70)-(72). For the systems (68)-(69) in matrix
form, we could consider the _ideal_ preconditioners

$\displaystyle\widetilde{P}:=A\otimes M_{1}+M\otimes B_{1}+B_{2}\otimes
M_{2};\qquad\widehat{P}:=A\otimes D_{1}^{-1}M_{0}+M_{0}\otimes
A_{1}+A_{2}\otimes D_{2}^{2}D_{1}^{-1}M_{0},$ (79)

which contain the discrete operators for the second-order derivative terms of
the stiffness matrices $\widetilde{A}$ and $\widehat{A}$, respectively.
Nevertheless, in matrix form the preconditioning operations
$P^{-1}\texttt{vec}(U)$, with $P$ equal $\widetilde{P}$ or $\widehat{P}$,
imply again the solution of another multiterm Sylvester equation at each
iteration of the PCG due to the presence of more Kronecker products in (79).
In order to avoid this, we approximate in a spectral sense the preconditioners
$\widetilde{P}$ and $\widehat{P}$ in (79), respectively, such that a single
Kronecker product is present. After several experiments, the choice
$\widetilde{P}_{1}:=(A+B_{2})\otimes B_{1}$ and
$\widehat{P}_{1}:=(A+A_{2})\otimes A_{1}$ proved the best among the tested
ones. In matrix form, these preconditioners translate to the following
operators:

$\displaystyle\widetilde{\mathcal{P}}(U):=B_{1}U(A+B_{2});\qquad\widehat{\mathcal{P}}(U):=A_{1}U(A+A_{2}).$
(80)

that satisfy
$\texttt{vec}(\widetilde{\mathcal{P}}(U))=\widetilde{P}_{1}\texttt{vec}(U)$
and $\texttt{vec}(\widehat{\mathcal{P}}(U))=\widehat{P}_{1}\texttt{vec}(U)$
for all $U\in\mathbb{R}^{q\times q}$. The respective inverse operators are
given by

$\widetilde{\mathcal{P}}^{-1}(U)=B_{1}^{-1}U(A+B_{2})^{-1};\qquad\widehat{\mathcal{P}}^{-1}(U)=A_{1}^{-1}U(A+A_{2})^{-1}.$

In the special case of square domains, the choices (80) reduce to

$\displaystyle\widetilde{\mathcal{P}}(U)=\widehat{\mathcal{P}}(U)=AUA,$ (81)

which will prove particularly efficient in the following experiments.

Figure 4: Poisson equation (66) on the square $\Omega=[0,1]^{2}$. FEM of order
$k=1,2,3,4$: comparison between the "vector" formulation (51) (dashed lines)
solved by the direct method and the MO-FEM solved by the MO-PCG (continuous
lines). Left plot: computational times. Right plot: Convergence behaviour of
the two approaches for all $k$. Dotted lines indicate slopes 2,4 and 5. For
$k=1,3,4$, the convergence is optimal of order $k+1$, for $k=2$, we observe
superconvergence of order 4. PCG is cheaper in execution time and for all $k$
and for all $N$ stops after two iterations, see also Table 1.

###  5.1 Example 1: Poisson equation on the square, $\mathbb{P}_{k}$ elements

We consider the Poisson equation (66) on the square $\Omega=[0,1]^{2}$. We
consider Lagrangian $\mathbb{P}_{k}$ finite elements for $k=1,2,3,4$, and we
compare the classical vector approach that solves the Kronecker form (51)
through the direct solver mldivide of MATLAB, and the MO-PCG approach
(77)-(78) using (81) as preconditioner. Observe that, since the domain is the
unit square, the multiterm Sylvester formulation (70) reduces to the classical
Sylvester equation (52) and then in this case the MO-PCG can be regarded also
as an alternative to the _reduced_ approach presented in Section 4. We
therefore consider the same sequence of seven meshes considered in Section 4.2
with $N=24\cdot 2^{i}$, $i=0,\dots,6$. On equal meshes, the MO-PCG with
preconditioning in (81) is quicker than the direct method for the vector
approach and the gap increases with $N$, as we can see in Fig. 4, left plot.
The detailed time comparisons on the finest mesh are shown in Table 1. For all
$k$, the methods exhibit optimal convergence ($(k+1)$-th order), with the case
$k=2$ being superconvergent (fourth order) as we can see in Fig. 4, right
plot. On the finest mesh ($N=1536,k=4$) errors near the machine precision
affect the convergence order, with MO-PCG being more accurate. For all $k$ and
for all $N$, PCG terminates with two iterations. In conclusion, the MO-PCG
method outperforms also the _reduced matrix_ approach (56) when the eigenvalue
decompositions (55) are computed numerically.

Table 1: Poisson equation (66) on the square $\Omega=[0,1]^{2}$: computational times on the finest mesh ($N=1536$) for all $k=1,2,3,4$ (see Fig.4, left plot) and the respective time ratios between the direct method for the vector form (51) and the MO-PCG method (77)-(78) with preconditioning in (81). $k$ |  Time (s) Direct method (vector form) |  Time (s) Matrix PCG |  Time ratio (Direct/PCG)
---|---|---|---
$1$ | 9.0273 | 0.2992 | 30.1726
$2$ | 9.4791 | 0.4255 | 22.2760
$3$ | 9.3019 | 0.4354 | 21.3665
$4$ | 10.1502 | 0.4438 | 22.8726

###  5.2 Example 2: Poisson equation on curved domain

Figure 5: Poisson equation (83) on the cap-shaped domain $\Omega^{S}$ defined
in (82). Continuous lines indicate MO-PCG, while dashed lines indicate the
vector (Kronecker) formulations solved by the direct method. In terms of
computational time, MO-PCG is competitive only for $k=1$, without lumping
(upper left panel). Both approaches exhibit optimal convergence for all FEM
order $k$, with the case $k=2$ being superconvergent, round-off error becomes
dominant on the finest mesh for $k=4$ (upper-right and lower-left panels).
Number of iterations required by the MO-PCG increases with $k$ and $N$, lower-
right panel.

Consider the _cap-shaped_ symmetric $x$-normal domain

$\Omega^{S}:=\left\\{(x,y)\in\mathbb{R}^{2}\middle|0\leq y\leq 1,\ |x|\leq
1-\frac{y^{2}}{2}\right\\},$ (82)

(shown in the next Fig. 6). On $\Omega^{S}$ we consider the Poisson equation
with zero Dirichlet boundary conditions:

$\begin{cases}-&\Delta u(x,y)=f(x,y),\qquad(x,y)\in\Omega^{S};\\\
&u(x,y)=0,\qquad(x,y)\in\partial\Omega^{S},\end{cases}$ (83)

were $f(x,y)$ is chosen in such a way that the exact solution is

$u(x,y)=y(y-1)\left(-\frac{y^{2}}{2}+x+1\right)\left(\frac{y^{2}}{2}+x-1\right),\qquad(x,y)\in\Omega^{S},$

we omit the cumbersome expression of such $f(x,y)$. We consider the
$\mathbb{P}_{k}$ elements (68), $k=1,\dots,4$, and the lumped $\mathbb{P}_{1}$
elements in vector form (69) and we solve again these classical linear systems
via the MATLAB direct solver mldivide. Hence, we compare the perfomance in
terms of computational time with the matrix-oriented PCG method (77)-(78) with
preconditioners as in (80), respectively, that solve the multiterm Sylvester
equations arising by the MO-FEM.

By applying all methods on a sequence of $6$ meshes $\Omega_{i}^{S}$,
$i=0,\dots,5$ with $N_{i}=24\cdot 2^{i}$ for all $i=0,\dots,5$ we find optimal
quadratic convergence in $L^{2}(\Omega^{S})$ and almost equal errors, as
reported in Fig. 5, upper right plot. By comparing the computational times, as
shown in Fig. 5, upper left plot, the MO-PCG is competitive only for $k=1$ and
moderate meshsizes $N<768$. Note that the no lumped FEM is slightly better
than the lumped version. On finer meshes, say for $N\geq 800$, the number of
iterations required by MO-PCG increases dramatically for $k\geq 2$, so
affecting its global performance. We guess that different preconditioners
could improve this behaviour, but this study is outside the scope of the
present work. On the other hand, we find in this particular experiment that,
the PCG method in its classical vector form is much slower the our MO-PCG, as
discussed before. In fact, by applying the built-in function _pcg_ of Matlab
with the same preconditioner, we find that MO-PCG is much quicker and the gap
increases with $N$. For example, for $k=1$ without lumping, the speedup factor
is: $1.48$ for $N=24$ and $30.82$ for $N=768$. As a final remark, it is worth
noting that in any case, as $N$ and $k$ increase, MO-PCG becomes more
competitive than the Direct (Kronecker) solver in terms of memory occupation,
as explained before in more detail.

In the next section we will show that the MO-PCG approach will prove far more
convenient in the case of time-dependent PDEs.

##  6 Time-dependent PDEs

On a general $x$-normal domain $\Omega^{L}$ we consider the following
semilinear heat equation

$u_{t}-d_{u}\Delta u=f(u,\boldsymbol{x},t),\qquad\boldsymbol{x}\in\Omega^{L},\
t\in[0,T],$ (84)

with $d_{u}>0$, endowed with either homogeneous Dirichlet boundary conditions
$u_{|\partial\Omega^{S}}=0$ or homogeneous Neumann boundary conditions
$(\nabla u\cdot\boldsymbol{n})_{|\partial\Omega^{S}}=0$. The treatment of RDSs
of the form (3) is completely analogous and we omit the details. The
$\mathbb{P}_{k}$ and lumped $\mathbb{P}_{1}$ FEM spatial discretisations of
(84) in vector form are as follows

$\displaystyle\widetilde{M}\boldsymbol{\xi}_{t}+d_{u}\widetilde{A}\boldsymbol{\xi}=\widetilde{M}f(\boldsymbol{\xi}),$
(85)
$\displaystyle\widehat{M}\boldsymbol{\xi}_{t}+d_{u}\widehat{A}\boldsymbol{\xi}=\widehat{M}f(\boldsymbol{\xi}),$
(86)

respectively, where $\boldsymbol{\xi}(t)\in\mathbb{R}^{q^{2}}$, with $q$ as
defined in (48), is the time-dependent nodal vector of the spatially discrete
solution. A full discretisation can be obtained by applying the IMEX Euler
timestepping scheme with timestep $\tau>0$ to (85) and (86), which yields

$\displaystyle(I+d_{u}\tau\widetilde{M}^{-1}\widetilde{A})\boldsymbol{\xi}^{(k+1)}=\boldsymbol{\xi}^{(k)}+\tau\boldsymbol{f}^{(k)},\qquad
k=0,\dots,N_{T}-1,$ (87)
$\displaystyle(I+d_{u}\tau\widehat{M}^{-1}\widehat{A})\boldsymbol{\xi}^{(k+1)}=\boldsymbol{\xi}^{(k)}+\tau\boldsymbol{f}^{(k)},\qquad
k=0,\dots,N_{T}-1,$ (88)

respectively, where $N_{T}:=\left\lceil\frac{T}{\tau}\right\rceil$,
$\boldsymbol{\xi}^{(k)}$ is the nodal vector of the fully discrete solution at
time $t_{k}:=k\tau$ and $\boldsymbol{f}^{(k)}:=f(\boldsymbol{\xi}^{(k)})$, see
[12]. Following [9], the fully discrete scheme (88) can be further accelerated
through an a-priori LU factorisation of the coefficient matrices (in
combination with symamd reordering of such matrix to further increase
sparsity). In matrix-oriented form, (87) and (88) become

$\displaystyle\widetilde{\mathcal{L}}(U^{(k+1)})=\widetilde{\mathcal{R}}(U^{(k)}),\qquad
k=0,\dots,N_{T},$ (89)
$\displaystyle\widehat{\mathcal{L}}(U^{(k+1)})=\widehat{\mathcal{R}}(U^{(k)}),\qquad
k=0,\dots,N_{T},$ (90)

respectively, where $U^{(k)}$ is such that
$\texttt{vec}(U^{(k)})=\boldsymbol{\xi}^{(k)}$ and

$\displaystyle\widetilde{\mathcal{L}}(U):=M_{3}UM+d_{u}\tau\Big{(}M_{1}UA+B_{1}UM+M_{2}UB_{2}-C_{2}UC_{1}-C_{2}^{T}UC_{1}^{T}\Big{)};$
(91)
$\displaystyle\begin{split}&\widehat{\mathcal{L}}(U):=M_{0}D_{1}^{2}UM_{0}+d_{u}\tau\Big{(}(M_{0}UA+D_{1}A_{1}UM_{0}+M_{0}D_{2}^{2}UA_{2})\\\
&\hskip 28.45274pt-
D_{1}(CD_{2}UCD_{3}+(CD_{2})^{T}U(CD_{3})^{T})\Big{)};\end{split}$ (92)
$\displaystyle\widetilde{\mathcal{R}}(U):=M_{3}(U+\tau f(U))M;$ (93)
$\displaystyle\widehat{\mathcal{R}}(U):=M_{0}D_{1}^{2}(U+\tau f(U))M_{0}.$
(94)

Each iteration of the fully discrete scheme (89) or (90) is a multiterm
Sylvester equation that can be solved through the matrix-oriented PCG method
(77)-(78) where $\mathcal{L}$ and $\mathcal{R}$ are chosen accordingly. For
the preconditioning operator $\mathcal{P}$, observe in (89) that, in the limit
$\frac{\tau}{N^{2}}\rightarrow 0$, the operator $\mathcal{L}(U)$ reduces to
$M_{3}UM$. Hence, $\mathcal{P}(U):=M_{3}UM$ would make for a reasonable
preconditioner. However, we experimentally found that by adding suitable
corrections we obtain a more accurate preconditioner that retains a single-
term form and it is given by:

$\displaystyle\widetilde{\mathcal{P}}^{-1}(U)=(M_{3}+d_{u}\tau
B_{1})^{-1}U(M+d_{u}\tau(A+B_{2}))^{-1};$ (95)
$\displaystyle\widehat{\mathcal{P}}^{-1}(U)=(M_{0}D_{1}+d_{u}\tau
A_{1})^{-1}U(M_{0}+d_{u}\tau(A+A_{2}))^{-1};$ (96)

We use the following stop criterion: at each timestep, we stop the iterations
of (78) when the truncated solution $U^{(s)}$ fulfils
$\|R^{(s)}\|\leq\tau\|R^{(0)}\|$, where $R^{(s)}$ is the residual defined by
$R^{(s)}:=\mathcal{L}(U^{(s)})-\mathcal{R}(U^{(s)})$, with the corresponding
operators as defined in (91)-(94). 111This stopping criterion guarantees
optimal convergence in space and time, and is justified as follows. At each
timestep, the initial guess $U^{(0)}$ is an $O(\tau)$-accurate approximation
of the exact solution $\widehat{U}$ of (90) which in turn contains an
$O(\tau^{2})$ discretisation error, since IMEX-Euler is first-order accurate.
Hence, to preserve the accuracy of the method, the truncated solution
$U^{(s)}$ of (77)-(78) must be an $O(\tau^{2})$-accurate approximation of
$\widehat{U}$ as well. Consequently, $U^{(s)}$ must approximate $\widehat{U}$
better than $U^{(0)}$ by $O(\tau)$ times. In terms of residuals, $\|R^{(s)}\|$
must be $O(\tau)$ times $\|R^{(0)}\|$ in (77)-(78).

###  6.1 Numerical Example: Semilinear heat equation on $x$-normal domain

We consider the following heat equation with zero Dirichlet boundary
conditions on the the cap-shaped domain $\Omega^{S}$ defined in (82) :

$\begin{cases}&u_{t}-d_{u}\Delta u=f(x,y,t),\qquad(x,y)\in\Omega^{S},\quad
t\in[0,1];\\\ &u(x,y,t)=0,\qquad(x,y)\in\partial\Omega^{S},\quad t\in[0,1];\\\
&u(x,y,0)=y(y-1)\left(-\frac{y^{2}}{2}+x+1\right)\left(\frac{y^{2}}{2}+x-1\right),\qquad(x,y)\in\Omega^{S},\end{cases}$
(97)

were $d_{u}=0.1$ and $f(x,y,t)$ is chosen in such a way that the exact
solution is $u(x,y,t)=u(x,y,0)\exp(t)$, we omit the cumbersome expression of
such $f(x,y,t)$. We consider both $\mathbb{P}_{k}$, $k=1,2,3,4$, and lumped
$\mathbb{P}_{1}$ elements. Also in this case, we solve the vector formulations
(87) (for $\mathbb{P}_{k}$ elements) and (88) (for lumped $\mathbb{P}_{1}$
elements) via the MATLAB direct solver mldivide, with only one preliminary LU-
decomposition. We compare these results with the matrix-oriented PCG approach
(90) with preconditioner (95). We present a Test 1, to study the convergence
of the two approaches and a Test 2 to highlight the computational advantages
in time of the MO-PCG approach.

Table 2: Test 1 \- Semilinear heat equation (97) on the cap-shaped domain: $\mathbb{P}_{k}$ finite elements, $k=1,2,3,4$ and lumped $\mathbb{P}_{1}$ elements. We show the convergence rates in space and time obtained by solving the vector formulations (87) and (88) by the direct method and the corresponding MO-FEM by the MO-PCG method (89). The convergence rates are optimal for all $k=1,2,3,4$. $k$ |  Convergence rate Vector method |  Convergence rate Matrix PCG
---|---|---
$1$ lumped | 1.9953 | 2.0509
$1$ | 1.9945 | 1.9946
$2$ | 2.9945 | 2.9947
$3$ | 3.9941 | 3.9943
$4$ | 4.9939 | 4.9942

In Test 1, for both methods we consider the second and third mesh
$\Omega_{i}^{S}$, $i=1,\dots,2$ of Experiment 5.2 for $N=48,96$.
Correspondingly, for each $i=1,2$ and $k=1,2,3,4$ we choose
$\tau_{i,k}=0.01\cdot 2^{(k+1)(i-1)}$ and $\texttt{tol}=1$. This choice of
timesteps allows to highlight optimal convergence in $L^{2}(\Omega^{S})$ norm
(i.e. $(k+1)$-th order in space and first order in time). We have confined
this test to these $N$ values because the timestep $\tau_{i,k}$ would become
too small for larger values of $N$ when $k>1$. The results are shown in Table
2.

In Test 2, we consider $N=480,960,1920$, fixed $\tau=1$e-2. This test is more
representative of typical user-case scenarios where high spatial resolution is
required. The obtained results indicate a significant advantage of MO-PCG and
are shown in Table 3 for all $k$ and for all $N$. At each timestep, MO-PCG
converges with just one iteration, except with lumped $\mathbb{P}_{1}$
elements, where up to two PCG iterations per timestep are required.

Table 3: Test 2 \- Semilinear heat equation (97) on the cap-shaped domain (82), solved through $\mathbb{P}_{k}$ finite elements, $k=1,2,3,4$ and lumped $\mathbb{P}_{1}$ elements. For each meshsize $N$, we apply the IMEX Euler method with timestep $\tau=0.01$. In all cases, the MO-PCG approach is quicker than the vector (direct) approach, especially for $k=1$ and $k=4$. The gap increases with $N$, as shown by comparing the time ratios. MO-PCG always converges with one single iteration, except with lumped $\mathbb{P}_{1}$ elements, where two iterations are required. $N$ | $k$ |  Time (s) Vector method |  Time (s) Matrix PCG |  Time ratio (Vector/PCG) |  Iterations PCG
---|---|---|---|---|---
$480$ |  $1$ lumped | 16.79 | 4.089 | 4.106 | 2
$1$ | 10.10 | 4.313 | 4.228 | 1
$2$ | 10.21 | 4.571 | 2.342 | 1
$3$ | 13.42 | 4.886 | 2.233 | 1
$4$ | 17.30 | 4.094 | 2.746 | 1
$960$ |  $1$ lumped | 129.1 | 18.58 | 6.949 | 2
$1$ | 131.7 | 16.41 | 8.023 | 1
$2$ | 54.98 | 19.10 | 2.879 | 1
$3$ | 58.79 | 20.25 | 2.902 | 1
$4$ | 89.18 | 21.62 | 4.126 | 1
$1920$ |  $1$ lumped | 1269 | 77.50 | 16.37 | 2
$1$ | 4337 | 74.35 | 58.33 | 1
$2$ | 497.9 | 87.01 | 5.722 | 1
$3$ | 473.0 | 81.17 | 5.828 | 1
$4$ | 3311 | 87.27 | 37.93 | 1

##  7 Applications to pattern formation in battery modeling

We now consider the following reaction-diffusion model in two variables
$\eta:\Omega\times[0,T]\rightarrow\mathbb{R}$ and
$\theta:\Omega\times[0,T]\rightarrow[0,1]$, endowed with zero Neumann boundary
conditions, on an arbitrary compact domain $\Omega\subset\mathbb{R}^{2}$:

$\begin{cases}\eta_{t}-\Delta\eta=\rho
f(\eta,\theta),\qquad(x,y,t)\in\Omega\times[0,T];\\\
\theta_{t}-d_{\theta}\Delta\theta=\rho
g(\eta,\theta),\qquad(x,y,t)\in\Omega\times[0,T];\\\
\nabla\eta\cdot\boldsymbol{n}=\nabla\theta\cdot\boldsymbol{n}=0,\qquad(x,y,t)\in\partial\Omega\times[0,T];\\\
\eta(x,y,0)=\eta_{0}(x,y),\quad\theta(x,y,0)=\theta_{0}(x,y),\qquad(x,y)\in\Omega,\end{cases}$
(98)

where $d_{\theta}>0$ is the diffusion coefficient, $\rho>0$ is a space-time
rescaling factor and the kinetics are

$\displaystyle f(\eta,\theta):=A_{1}(1-\theta)\eta-
A_{2}\eta^{3}-B(\theta-\alpha);$ (99) $\displaystyle
g(\eta,\theta):=C(1+k_{2}\eta)(1-\theta)[1-\gamma(1-\theta)]-D\theta(1+\gamma\theta)(1+k_{3}\eta),$
(100)

with $\alpha,\gamma,A_{1},A_{2},B,C,D,k_{2},k_{3}$ positive parameters.

The PDE system (98)-(100) is known as DIB model and has been introduced for
the first time in [5] to describe electrodeposition processes. Under suitable
choices of the parameters and of the domain $\Omega$, this model was shown to
possess a variety of spatially structured solutions, known as Turing patterns,
see for example [20, 21]. An interesting application in battery modeling is
reported in [19, 22]. Turing patterns are obtained as stationary solutions of
(98)-(100) and then their numerical approximation requires highly spatial
accuracy for longtime integration, this motivates the development of efficient
solvers. In this direction, a first work based on matrix oriented formulation
of (98)-(100) is [9] where finite differences and several time solvers have
been proposed on square domains. In [9], the Sylvester matrix equations
obtained at each time step have been approximated by the _reduced approach_ ,
similar to the one in Section 4. For example the IMEX Euler yielded the
_rEuler_ method, that revealed much more efficient than its classical vector
approach. On the other hand, domain geometry was also proven to play an
important role in pattern selection, as shown also in [19, 22], for this
reason efficient solvers that can be applied on domains as general as possible
are need. Towards this aim, here we propose the matrix oriented FEM spatial
approximation and the MO-PCG approach presented in Setion 6 with
preconditioner (95) to deal in particular with some x-normal domains. We will
present two kind of simulations, first on the cap shaped domain introduced in
(82) and then on the jar-domain shown in Fig. 2(a) that correspond to the
curvilinear cylinder in 2(b). In both cases we will consider domains of
increasing _effective domain size_ given by ${\cal A}=\rho|\Omega|$, where
$|\Omega|$ is the area of the domain in (98). In fact, as shown in [19, 22]
there exists a sufficiently large $\cal{A}^{*}$ such that for
$\cal{A}>\cal{A}^{*}$ the _intrinsic_ Turing pattern corresponding to the
fixed model parameters arises (see also [19] for more details), otherwise only
a portion of it can be approximated, giving rise to doubts about its
classification.

To solve on domains of large sizes, we exploit the meaning and the role of the
parameter $\rho$ in (98), as follows. By introducing new variables
$(\widetilde{x},\widetilde{y})=\sqrt{\rho}(x,y)$, the chain rule yields

$\displaystyle\frac{\partial\eta}{\partial
t}=\rho\frac{\partial\eta}{\partial\widetilde{t}},\qquad\nabla_{(x,y)}\eta=\sqrt{\rho}\nabla_{(\widetilde{x},\widetilde{y})}\eta,\qquad\Delta_{(x,y)}\eta=\rho\Delta_{(\widetilde{x},\widetilde{y})}\eta;$
(101) $\displaystyle\frac{\partial\theta}{\partial
t}=\rho\frac{\partial\theta}{\partial\widetilde{t}},\qquad\nabla_{(x,y)}\theta=\sqrt{\rho}\nabla_{(\widetilde{x},\widetilde{y})}\theta,\qquad\Delta_{(x,y)}\theta=\rho\Delta_{(\widetilde{x},\widetilde{y})}\theta;$
(102)

Hence, we define $\Omega_{\rho}:=\sqrt{\rho}\Omega$ and $T_{\rho}:=\rho T$.
Hence, $\rho$ acts as a rescaling parameter in space and time.

In the following simulations we always solve 110 in the reference domanis
$\Omega$ cap shaped e jar shaped with final time time $T$ and timespet $\tau$
such that $T_{\rho}=300$ and $\tau_{\rho}=5e-3$ which guarantees the stability
of the IMEX Euler method. The corresponding numerical solutions will be
plotted in the rescaled domain $\Omega_{\rho}$ for $T_{\rho}=300$. In all the
experiments fix the following model parameters:

$\alpha=0.5,\gamma=0.2,A_{1}=10,D=3.2727,k_{2}=2.5,k_{3}=1.5,d_{\theta}=20.$
(103)

The initial data are given by $\theta_{0}(x,y)=\theta_{e}+10^{-4}rand(x,y)$
and $\eta_{0}(x,y)=\eta_{e}+10^{-4}rand(x,y)$ and are small spatially random
perturbations of the homogeneous equilibrium $(\eta_{e},\theta_{e}):=(0,0.5)$.

###  7.1 Cap-shaped domain

In this example, we consider the cap-shaped domain (82) and we choose
$A_{2}=30,B=25,C=7$, that, according to the segmentation results in [29], can
yield mixed spots-worms Turing patterns at the steady state.

We solve the PDE RDS system with $\mathbb{P}_{k}$ elements ($k=1,2,3,4$) and
lumped $\mathbb{P}_{1}$ elements in space and by the IMEX Euler, as in the
previous section. We compare the vector approach solving the sequence of
linear systems in (87)-(88) by the direct method (that we will call "vector
method") with the MO-PCG approach (90), solving the multiterm Sylvester
equations arising at each time step for this choice of the domain.

We solve the DIB model (98) with different combinations of $\rho$, $N_{x}$ and
$N_{y}$ as listed in Table 4, that is for domains of larger area
$\mathcal{A}=|\Omega_{\rho}|$. In all the computations, we discretise the $x$
dimension with $N_{x}$ nodes and the $y$ dimension with $N_{y}$ nodes, with
$N_{x}=2N_{y}$, which reflects the aspect ratio of the domain. In Fig. 6, for
each simulation, we report the final patterns obtained by the vector approach
and by the MO-PCG, together with the respective increments
$\|\eta^{(k+1)}-\eta^{(k)}\|_{F}$ as a function of time, with $\|\cdot\|_{F}$
Frobenius norm. If such increment decreases over time and tends to an almost
small stationary value, we deduce that the numerical solution is converging to
a steady state.

For increasing values of the _effective domain size_ $\cal A$ the solution
morphology changes and a pattern with more structures is attained, as shows in
Fig.6,(a)–(c) corresponding to the values (a)–(c) in Table 4, respectively. In
simulation (a) a good pattern is attained by both methods, but its morphology
is not completely expressed. The vector and the matrix approach seem to be
equivalent in this case also in terms of computational times (see Table 4). To
capture the true Turing morphology a larger domain $\Omega$ and a
_sufficiently_ fine mesh is required, otherwise _phantom patterns_ could be
obtained. This is exactly what happens for the simulation in case (b),
corresponding to the second row of both Fig.6 and Table 4, where the same mesh
of case (a) yields a “pixelated” pattern. Hence, for the same domain, in
simulation (c) a finer mesh is used and both methods are finally able to
attain a “complete” pattern. Note that, in vector form at each time step we
solve a problem of dimension $N_{x}\cdot N_{y}=200\cdot 400=8\cdot 10^{4}$, by
the MO-PCG instead we solve a sequence of $N_{t}=6\cdot 10^{4}$ rectangular
multiterm Sylvester equations of size $400\times 200$. Moreover, this example
shows that the matrix-oriented PCG algorithm (77)-(78) can successfully solve
rectangular Sylvester equations. As shown in Table 4, the time ratios indicate
that the matrix PCG approach (89)-(90) tends to become quicker than the
vector-direct approach (88) with significant advantage only for
$\mathbb{P}_{4}$ and lumped and no lumped $\mathbb{P}_{1}$ elements. For
$k=2,3$, we guess that a different preconditioner could improve the results of
the MO-PCG method.

Table 4:  DIB model (98)-(100) with parameters (103) on the cap-shaped domain (82). Parameters for the simulations (a)-(b)-(c) reported in Figure 6 and performance comparison between the vector approach (87)-(88) based on LU decomposition and the matrix PCG method (90). When the time ratio $r_{t}$ is close to 1, the methods take approximately the same time; $r_{t}>1$ indicates that MO-PCG is quicker. The MO-PCG is less expensive as $N_{x}$ and $N_{y}$ increase. Simulation in (b) yields a _phantom pattern_. The last column shows the amount of iterations required by MO-PCG for each of the two PDEs of the model. | $\cal{A}$ | $\rho$ | $N_{x}$ | $N_{y}$ | $k$ |  Time (s) Vector method |  Time (s) Matrix PCG |  Time ratio $r_{t}$ |  Iterations PCG $(\eta,\theta)$
---|---|---|---|---|---|---|---|---|---
(a) | $2000/3$ | $400$ | $100$ | $50$ |  $1$ lumped | 57.18 | 93.76 | 0.6099 | (1,6)
$1$ | 62.76 | 94.60 | 0.6635 | (2,3)
$2$ | 75.60 | 182.0 | 0.4159 | (2,4)
$3$ | 226.6 | 193.4 | 1.172 | (2,3)
$4$ | 270.2 | 207.7 | 1.301 | (2,3)
(b) | $10000/3$ | $20000$ | $100$ | $50$ |  $1$ lumped | 58.08 | 47.53 | 1.222 | (1,2)
$1$ | 69.72 | 70.89 | 0.9835 | (1,3)
$2$ | 75.50 | 118.1 | 0.6394 | (1,3)
$3$ | 82.79 | 142.6 | 0.5804 | (1,3)
$4$ | 273.7 | 156.5 | 1.749 | (1,3)
(c) | $10000/3$ | $20000$ | $400$ | $200$ |  $1$ lumped | 1590 | 879.8 | 1.807 | (1,4)
$1$ | 1623 | 1249 | 1.299 | (2,3)
$2$ | 1809 | 1934 | 0.9352 | (1,3)
$3$ | 1973 | 2143 | 0.9209 | (2,3)
$4$ | 6323 | 2266 | 2.791 | (2,3)

Figure 6: DIB model (98)-(100) on the cap-shaped domain, lumped
$\mathbb{P}_{1}$ solutions. The solutions are plotted on the rescaled domain
$\Omega_{\rho}$. Each row corresponds to the $(\mathcal{A},\rho,N_{x},N_{y})$
combination in Table 4. For $\rho=400$ (a) only few structures arise in the
pattern. For $\rho=20000$: in (b) $N_{x}=200,N_{y}=100$ are not sufficient to
resolve the pattern structure, that instead is well-resolved for $N_{x}=400$
and $N_{y}=200$ in (c). Smaller values of $N_{x}$ and $N_{y}$ do not provide
sufficient spatial approximation on the larger domain and the pattern appear
grainy (b). The patterns obtained by the MO-PCG method and the vector (direct)
method are very similar, they are stationary solutions obtained at
$T_{\rho}=300$, are shown by the increment dynamics (left subplots). To
compare the execution times by the lumped $\mathbb{P}_{1}$ approximation see
Table 4.

###  7.2 Jar-shaped domain

Thanks to the results in the previous test, here we solve the model only with
$\mathbb{P}_{1}$ elements, both with and without lumping. To further explore
the robustness of the matrix PCG approach w.r.t. _domain complexity and mesh
distortion_ , we consider the _jar-shaped domain_ in Fig. 2(a). We fix the
parameters $A_{2}=1,B=30,C=3$ for the DIB model which are known to produce
Turing patterns with holes (also called reversed spots) [29] and again we
solve for increasing $\mathcal{A}$ with different combinations of $\rho$,
$N_{x}$ and $N_{y}$ as listed in Table 5. In all the computations, we consider
$N_{x}=3N_{y}$, which reflects the aspect ratio of the domain. The timestep
and the final time are as in the previous test on the cap-shaped domain. We
show the $\mathbb{P}_{1}$ solutions in Fig. 7.
As we can see in the figure, in case (a) on the smallest domain only few holes
arise in the pattern. In case (b) on the larger domain, more structures arise
in the pattern, but the mesh is too coarse and a _phantom pattern_ arises. In
case (c), the intrinsic Turing pattern is well-resolved for $N_{x}=600$ and
$N_{y}=200$. The MO-PCG and the vector solutions are very similar, but the MO
approach converges in significant less time (see Table 5).

We conclude by remarking that, since the jar-shaped domain in Fig. 2(a) can be
transformed to the cylinder $\Gamma$ shown in Fig. 2(b), then the solutions
shown in Fig. 7 can be interpreted, after the coordinate transformation (21),
as solutions to the _surface DIB model_ , that is (98)-(100) where the Laplace
operator $\Delta$ is replaced by the Laplace-Beltrami operator
$\Delta_{\Gamma}$ on the cylinder $\Gamma$. As an example, we report in Fig. 8
the solution in (a) wrapped on a curvilinear cylinder. The application of the
model on a cylindrical surface can be of applicative interest as shown in [6],
in which the authors consider the use of cylindrical Zn sponges as a means of
limiting the shape change and dendrite formation issues in Zn-based
rechargeable batteries.

Table 5: Reaction-diffusion model (98)-(100) with parameters as in Section 7.2 on the jar-shaped domain in Fig. 2(a): performance comparison between the vector method (87)-(88) with LU decomposition and the matrix PCG method (90). The vector PCG approach tends to become quicker than the vector approach as $N_{x}$ and $N_{y}$ increase. | $\mathcal{A}$ | $\rho$ | $N_{x}$ | $N_{y}$ | $k$ |  Time (s) Vector method |  Time (s) Matrix PCG | Time ratio |  Iterations PCG $(\eta,\theta)$
---|---|---|---|---|---|---|---|---|---
(a) | $800$ | $400$ | $150$ | $50$ |  $1$ lumped | 309.6 | 306.6 | 1.010 | (2,16)
$1$ | 328.7 | 278.4 | 1.181 | (3,7)
(b) | $40000$ | $20000$ | $150$ | $50$ |  $1$ lumped | 98.54 | 67.17 | 1.4670 | (1,2)
$1$ | 111.4 | 117.0 | 0.9528 | (1,2)
(c) | $40000$ | $20000$ | $600$ | $200$ |  $1$ lumped | 7995 | 2316.4 | 3.451 | (1,9)
$1$ | 8067 | 3028 | 2.664 | (2,7)

Figure 7: DIB model (98)-(100) on the jar-shaped domain, $\mathbb{P}_{1}$
solutions. Values of $(\rho,N_{x},N_{y})$ are given in Table 5. In case (a) on
the smallest domain only few holes arise in the pattern. In case (b), for
$\rho=20000$ on the larger domain the mesh is too coarse and a _phantom
pattern_ arises. In case (c), the intrinsic Turing pattern is well-resolved
for $N_{x}=600$ and $N_{y}=200$. MO-PCG and vector solutions are very similar,
but the MO approach converges in significant less time (see Table 5). Figure
8: Numerical solution of Fig. 7 (a) ($\rho=20000$, $N_{x}=600$, $N_{y}=200$)
mapped onto the cylinder $\Gamma$ in Fig. 2(b). It can be interpreted as the
solution of the DIB _surface reaction-diffusion model_ on $\Gamma$.

##  8 Conclusions

In this work we have provided a matrix-oriented formulation for Lagrangian
finite elements of arbitrarily high order $k\in\mathbb{N}$ on $x$-normal
domains. The proposed approach applies to both elliptic and parabolic PDE
problems. The discrete problems take the form of a matrix equation (or a
sequence of matrix Sylvester equations in the time-dependent case) of much
smaller dimension that is mathematically equivalent to the much larger
standard linear systems in Kronecker form. The proposed approach adopts a
curvilinear structured mesh that eliminates geometric boundary error.
Moreover, through a coordinate transformation, our approach applies also to
special surface domains, namely cylinders with arbitrary curved boundaries.
On square domains, the discrete problems take the form of a generalised two-
term Sylvester equation that we solve efficiently through a spectral approach
for all $k$. In this sense, our work extends the findings in [9], based on
classical finite differences, to the case of high order FEM in space.
On general $x$-normal domains, the discrete problem takes the form of (a
sequence of) _multiterm_ Sylvester equations which we solve through a matrix-
oriented PCG method with matrix-oriented preconditioner that is quick to
evaluate thanks to is single-term form. On one hand, such solver is always
quicker than the classical PCG in vector form. On the other hand, we show by
several numerical tests that it is quicker than the optimised direct solver
mldivide of MATLAB in the case of (i) time-dependent PDEs on general
$x$-normal domains and (ii) elliptic PDEs on square domains. In terms of
memory occupation, the matrix-oriented PCG method always improves on any
direct or iterative solver that relies on the full storage of the Kronecker
matrix, and the gap increases both with the number of gridpoints $N$ and the
polynomial order $k$ of the method. The only case in which we could not find
any speedup in the MO-PCG approach is that of elliptic problems on non-square
domains with $\mathbb{P}_{k}$ elements, $k\neq 1$. This opens the quest for
efficient preconditioners, which will be addressed in future studies.
Special consideration deserves the application to reaction-diffusion systems,
where the simulation of fine-grained patterns requires high spatial
resolution, which translates into computationally intense simulations both in
time and memory. Our experiments for the approximation of Turing patterns
arising in batteries as solution of the DIB morphochemical model provide
encouraging results in this direction and justify the matrix approach in terms
of execution times and storage. The best performance gains were observed with
$\mathbb{P}_{4}$ and lumped $\mathbb{P}_{1}$ elements. Also in this case, we
believe that further performance gains can be found through the development of
more efficient preconditioners and more efficient solvers for multiterm
Sylvester equations, such as a truncated PCG [30]. These aspects will be
addressed in future studies.

## Acknowledgments

The work of MF was funded by Regione Puglia (Italy) through the research
programme REFIN-Research for Innovation (protocol code 901D2CAA, project
number UNISAL026).
The work of IS has been funded by the MIUR (Italian Ministry of Education,
University and Research) project PRIN 2017, “Mathematics of active materials:
From mechanobiology to smart devices”, project no. 2017KL4EF3.
The work of MF and IS was performed under the auspices of GNCS-INdAM (Italian
National Group of Scientific Computing).

## References

  * [1] P Antolin, A Buffa, F Calabró, M Martinelli, and G Sangalli.  Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization.  Computer Methods in Applied Mechanics and Engineering, 285:817–828, 2015.  doi:10.1016/j.cma.2014.12.013.
  * [2] R Barreira, C M Elliott, and A Madzvamuse.  The surface finite element method for pattern formation on evolving biological surfaces.  Journal of Mathematical Biology, 63(6):1095–1119, 2011.  doi:10.1007/s00285-011-0401-0.
  * [3] R H Bartels and G W Stewart.  Solution of the matrix equation AX + XB = C.  Communications of the ACM, 15(9):820–826, 1972.  doi:10.1145/361573.361582.
  * [4] D Becherer, M Schweizer, et al.  Classical solutions to reaction–diffusion systems for hedging problems with interacting Itô and point processes.  Annals of Applied Probability, 15(2):1111–1144, 2005.  doi:10.1214/105051604000000846.
  * [5] B Bozzini, D Lacitignola, and I Sgura.  Spatio-temporal organization in alloy electrodeposition: a morphochemical mathematical model and its experimental validation.  Journal of Solid State Electrochemistry, 17(2):467–479, 2013.  doi:10.1007/s10008-012-1945-7.
  * [6] B Bozzini, C Mele, A Veneziano, N Sodini, G Lanzafame, A Taurino, and L Mancini.  Morphological evolution of Zn-sponge electrodes monitored by in situ X-ray computed microtomography.  ACS Applied Energy Materials, 3(5):4931–4940, 2020.  doi:10.1021/acsaem.0c00489.
  * [7] M A J Chaplain, M Ganesh, and I G Graham.  Spatio-temporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth.  Journal of Mathematical Biology, 42(5):387–423, 2001.  doi:10.1007/s002850000067.
  * [8] F Chung and S-T Yau.  Discrete Green’s functions.  Journal of Combinatorial Theory, Series A, 91(1-2):191–214, 2000\.  doi:10.1006/jcta.2000.3094.
  * [9] M C D’Autilia, I Sgura, and V Simoncini.  Matrix-oriented discretization methods for reaction–diffusion PDEs: Comparisons and applications.  Computers & Mathematics with Applications, 79(7):2067–2085, 2020\.  doi:10.1016/j.camwa.2019.10.020.
  * [10] C Eilks and C M Elliott.  Numerical simulation of dealloying by surface dissolution via the evolving surface finite element method.  Journal of Computational Physics, 227(23):9727–9741, 2008.  doi:10.1016/j.jcp.2008.07.023.
  * [11] C M Elliott and B Stinner.  Modeling and computation of two phase geometric biomembranes using surface finite elements.  Journal of Computational Physics, 229(18):6585–6612, 2010.  doi:10.1016/j.jcp.2010.05.014.
  * [12] M Frittelli, A Madzvamuse, I Sgura, and C Venkataraman.  Preserving invariance properties of reaction–diffusion systems on stationary surfaces.  IMA Journal of Numerical Analysis, 39(1):235–270, 2019.  doi:10.1093/imanum/drx058.
  * [13] G Golub, S Nash, and C Van Loan.  A Hessenberg-Schur method for the problem AX + XB = C.  IEEE Transactions on Automatic Control, 24(6):909–913, 1979.  doi:10.1109/tac.1979.1102170.
  * [14] G H Golub and C F Van Loan.  Matrix computations, volume 3.  JHU press, 2013.
  * [15] Yue Hao and Valeria Simoncini.  Matrix equation solving of PDEs in polygonal domains using conformal mappings.  Journal of Numerical Mathematics, 0(0), nov 2020.  doi:10.1515/jnma-2020-0035.
  * [16] T J R Hughes.  The finite element method: linear static and dynamic finite element analysis.  Courier Corporation, 2012.
  * [17] C Jordan and K Jordán.  Calculus of finite differences, volume 33.  American Mathematical Soc., 1965.
  * [18] E J Kansa.  Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics – II solutions to parabolic, hyperbolic and elliptic partial differential equations.  Computers & Mathematics with Applications, 19(8-9):147–161, 1990\.  doi:10.1016/0898-1221(90)90271-k.
  * [19] D Lacitignola, B Bozzini, M Frittelli, and I Sgura.  Turing pattern formation on the sphere for a morphochemical reaction-diffusion model for electrodeposition.  Communications in Nonlinear Science and Numerical Simulation, 48:484–508, 2017.  doi:10.1016/j.cnsns.2017.01.008.
  * [20] D Lacitignola, B Bozzini, and I Sgura.  Spatio-temporal organization in a morphochemical electrodeposition model: analysis and numerical simulation of spiral waves.  Acta Applicandae Mathematicae, 132(1):377–389, 2014.  doi:10.1007/s10440-014-9910-3.
  * [21] D Lacitignola, B Bozzini, and I Sgura.  Spatio-temporal organization in a morphochemical electrodeposition model: Hopf and Turing instabilities and their interplay.  European Journal of Applied Mathematics, 26(2):143–173, 2015.  doi:10.1017/s0956792514000370.
  * [22] D Lacitignola, I Sgura, B Bozzini, T Dobrovolska, and I Krastev.  Spiral waves on the sphere for an alloy electrodeposition model.  Communications in Nonlinear Science and Numerical Simulation, 79:104930, 2019.  doi:10.1016/j.cnsns.2019.104930.
  * [23] A Mantzaflaris, B Jüttler, B N Khoromskij, and U Langer.  Low rank tensor methods in Galerkin-based isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 316:1062–1085, 2017.  doi:10.1016/j.cma.2016.11.013.
  * [24] Y-Y Nie and V Thomée.  A lumped mass finite-element method with quadrature for a non-linear parabolic problem.  IMA Journal of Numerical Analysis, 5(4):371–396, 1985.  doi:10.1093/imanum/5.4.371.
  * [25] D Palitta and V Simoncini.  Matrix-equation-based strategies for convection–diffusion equations.  BIT Numerical Mathematics, 56(2):751–776, 2016.  doi:10.1007/s10543-015-0575-8.
  * [26] C E Powell, D Silvester, and V Simoncini.  An efficient reduced basis solver for stochastic Galerkin matrix equations.  SIAM Journal on Scientific Computing, 39(1):A141–A163, 2017.  doi:10.1137/15m1032399.
  * [27] Y Saad.  Iterative methods for sparse linear systems.  SIAM, 2003.
  * [28] G Sangalli and M Tani.  Isogeometric preconditioners based on fast solvers for the Sylvester equation.  SIAM Journal on Scientific Computing, 38(6):A3644–A3671, 2016.  doi:10.1137/16m1062788.
  * [29] I Sgura, A S Lawless, and B Bozzini.  Parameter estimation for a morphochemical reaction–diffusion model of electrochemical pattern formation.  Inverse Problems in Science and Engineering, 27(5):618–647, 2019\.  doi:10.1080/17415977.2018.1490278.
  * [30] S D Shank, V Simoncini, and D B Szyld.  Efficient low-rank solution of generalized Lyapunov equations.  Numerische Mathematik, 134(2):327–342, 2016.  doi:10.1007/s00211-015-0777-7.
  * [31] V Simoncini.  Computational methods for linear matrix equations.  SIAM Review, 58(3):377–441, 2016.  doi:10.1137/130912839.
  * [32] A M Turing.  The chemical basis of morphogenesis.  Bulletin of Mathematical Biology, 52(1):153–197, 1990.  doi:10.1093/oso/9780198250791.003.0022.
  * [33] V K Vanag.  Waves and patterns in reaction–diffusion systems. Belousov–Zhabotinsky reaction in water-in-oil microemulsions.  Physics-Uspekhi, 47(9):923, 2004.  doi:10.1070/pu2004v047n09abeh001742.