# A spectral Galerkin exponential Euler time-stepping scheme for parabolic
SPDEs on two-dimensional domains with a $\mathcal{C}^{2}$ boundary

Julian Clausnitzer
Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH, 52425 Jülich
<EMAIL_ADDRESS>
Andreas Kleefeld
Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH, 52425 Jülich
Faculty of Medical Engineering and Technomathematics, University of Applied
Sciences Aachen, 52428 Jülich
<EMAIL_ADDRESS>

###### Abstract

We consider the numerical approximation of second-order semi-linear parabolic
stochastic partial differential equations interpreted in the mild sense which
we solve on general two-dimensional domains with a $\mathcal{C}^{2}$ boundary
with homogeneous Dirichlet boundary conditions. The equations are driven by
Gaussian additive noise, and several Lipschitz-like conditions are imposed on
the nonlinear function. We discretize in space with a spectral Galerkin method
and in time using an explicit Euler-like scheme. For irregular shapes, the
necessary Dirichlet eigenvalues and eigenfunctions are obtained from a
boundary integral equation method. This yields a nonlinear eigenvalue problem,
which is discretized using a boundary element collocation method and is solved
with the Beyn contour integral algorithm. We present an error analysis as well
as numerical results on an exemplary asymmetric shape, and point out
limitations of the approach.

##  1 Introduction.

We consider the numerical approximation of semilinear parabolic stochastic
partial differential equations (SPDEs) of the form

$\displaystyle\mathrm{d}U(t)$ $\displaystyle=[\Delta
U(t)+F(U(t))]\mathrm{d}t+\mathrm{d}W^{Q}(t),$ $\displaystyle U(0)$
$\displaystyle=u_{0}\in H,\quad t\in[0,T],$

where $T>0$, $H\subset L^{2}(D)$ is a Hilbert space of functions satisfying a
homogeneous Dirichlet boundary condition on a domain $D\subset\mathbb{R}^{2}$
with a $\mathcal{C}^{2}$ boundary, $\Delta$ is the Laplace operator, $F$ is a
nonlinear continuous operator and $W^{Q}$ is a cylindrical $Q$-Wiener process
with covariance operator $Q$ (for more precise definitions, see §2). We make
the common assumption that the eigenfunctions of $Q$ and $\Delta$ coincide,
which allows for a series expansion of $W^{Q}$ in terms of eigenfunctions of
$\Delta$ (see Appendix A).
The numerical solution of SPDEs has attracted much attention, and many
different methods for temporal as well as spatial discretization have been
developed (for a comprehensive overview, see e.g. [38, Section 1.2]). However,
due to the irregular behavior of SPDEs, strong regularity conditions have to
be imposed on either the noise or the spatial domain. Some publications
focused on a one-dimensional spatial domain and potentially space-time white
noise [1, 2, 6, 9, 17, 22, 28, 29, 30, 42], while others have worked with
trace-class or smooth noise and a potentially two-dimensional domain [13,
6.24][24, 25, 26, 27, 36]. We intend to use the latter approach by considering
additive noise that has sufficiently quickly decaying modes, but in exchange
being able to consider a much larger class of general two-dimensional domains
$D\subset\mathbb{R}^{2}$ with a $\mathcal{C}^{2}$ boundary. This broadening of
the class of spatial domains has been pointed out in [38, p. 260] as a subject
of future research. Also, while there has been some theoretical treatment of
SPDEs on domains with a $\mathcal{C}^{2}$ boundary of dimension two or more
with homogeneous Dirichlet boundary conditions [8, 18, 19], to our knowledge
no solutions of SPDEs have been numerically approximated in such domains.
Spectral methods have been widely used for the solution of semilinear
parabolic SPDEs [6, 24, 28, 39] and are straightforward to implement in cases
where the base functions are known explicitly. If this is not the case, the
spectrum of the linear differential operator as well as its eigenfunctions
have to be numerically approximated. This can be achieved using boundary
element methods (BEM) resulting in a nonlinear eigenvalue problem. The
advantage of this method is that any bounded domain with a $\mathcal{C}^{2}$
boundary can be considered, and the eigenfunction can be computed pointwise at
any desired point in the interior of the domain, making it possible to store
the functions on a uniform grid. Another possibility to compute the base
functions are finite element methods (FEM), and for the Dirichlet Laplacian
this has been done in [40]. The difference to BEM is that instead of domains
with a $\mathcal{C}^{2}$ boundary, FEM can deal with polygons with a Lipschitz
boundary in $\mathbb{R}^{2}$. Since FEM requires the triangularization of the
domain, if the functions are to be stored on a uniform grid, an additional
error would be introduced by evaluating the functions on the triangularization
at the required grid points. In this work, we focus on BEM and domains with
$\mathcal{C}^{2}$ boundaries.
We give in §2 a precise mathematical framework, including the conditions
placed on each component of the equations which can be solved. In §3, the
discretization in space and time of the mild solutions is given as a Galerkin
projection in space and the exponential Euler scheme in time as described by
Jentzen and Kloeden [29]. The boundary integral equation method is recalled,
including the necessary discretization steps turning the Helmholtz equation
into a nonlinear eigenvalue problem. The algorithm developed by Beyn [5],
which allows for the solution of nonlinear eigenvalue problems, is recalled,
again also sketching the discretization steps necessary for the
implementation. Lastly, some remarks on assembling the first elements of an
orthonormal basis are made and it is shown how many boundary elements are
needed to compute the eigenfunctions to a given accuracy. In §4, we prove a
result on the strong error of the spectral Galerkin-exponential Euler scheme
with eigenvalues and eigenfunctions which are approximated to a given
accuracy. In §5, we demonstrate the convergence on an asymmetric shape and
present a few solutions for different nonlinear functions. The MATLAB program
code, including the data and programs to generate the plots shown in the
figures, is available at GitHub https://github.com/JulianSPDE/2D-SPDE.

##  2 Mathematical setting and assumptions.

Let $D\subset\mathbb{R}^{2}$ be a bounded domain whose boundary is of class
$\mathcal{C}^{2}$ and can be described by a parametrization $\gamma$ such that
$\|\gamma^{\prime}(t)\|_{2}>0$, $t\in[0,2\pi)$, $\gamma(0)=\gamma(2\pi)$ and
$\gamma$ is nonintersecting, i.e. $\gamma(s)\neq\gamma(t)$ for distinct
$s,t\in[0,2\pi)$. Let $T>0$ be fixed. We consider SPDEs of the form

$\displaystyle\begin{split}\mathrm{d}U(t)=[AU(t)+F(U(t))]\mathrm{d}t+\mathrm{d}W^{Q}(t),\quad
U(0)=u_{0},\quad t\in[0,T]\end{split}$ (1)

where $(U(t))_{t\in[0,T]}$ is a stochastic process taking values in a Hilbert
space $H$, $u_{0}\in H$, $A:D(A)\subset H\rightarrow H$ is a linear sectorial
operator, i.e. it generates an analytic semigroup
$(\mathrm{e}^{At})_{t\in[0,T]}$, and $F:H\rightarrow H$ is a nonlinear
operator. The process $(W^{Q}(t))_{t\in[0,T]}$ is a cylindrical $Q$-Wiener
process taking values in $H$, where $Q:H\rightarrow H$ is a trace-class
operator. We will consider the Laplacian $A=\Delta:D(\Delta)\subset
H\rightarrow H$ with homogeneous Dirichlet boundary conditions on the domain
$D$, so that we have $H=L^{2}(D)$ and $D(A)=H_{0}^{1}(D)\cap H^{2}(D)$. It is
also assumed that $\Delta$ and the covariance operator $Q$ commute, so that
they have a complete set consisting of the same orthonormal eigenfunctions
$e_{i}$, $i\in\mathbb{N}$. This assumption of a joint eigenbasis for $Q$ and
$\Delta$ has been treated in [7] and we show in Appendix A that this
assumption can be satisfied for more general two-dimensional shapes.
The SPDE (1) can be rigorously understood as an integral equation involving a
stochastic $\mathrm{d}W^{Q}$ integral [16, Chapter 4]. It is interpreted in
the mild sense [16, p. 161], meaning that an $H$-valued process
$(U(t))_{t\in[0,T]}$ is a solution of (1) if

$\displaystyle
U(t)=\mathrm{e}^{At}u_{0}+\int_{0}^{t}\mathrm{e}^{A(t-s)}F(U(s))\;\mathrm{d}s+\int_{0}^{t}\mathrm{e}^{A(t-s)}\;\mathrm{d}W^{Q}(s).$
(2)

We now give some details on the properties of each component of the SPDE (1).
The linear sectorial operator $A=\Delta$.
It is known that for Dirichlet, Neumann and mixed boundary conditions, the
Laplacian has a complete set of orthonormal smooth eigenfunctions $e_{i}$,
$i\in\mathbb{N}$. Also, the eigenvalues $-\lambda_{i}$ are all negative and
real having exactly one accumulation point at $-\infty$ (see e.g. [4, Appendix
W]; we use a notation where $\lambda_{i}>0$, $i\in\mathbb{N}$). Hence, it is
possible to give a spectral representation

$\displaystyle\textstyle\Delta f(x)=\sum_{i=1}^{\infty}-\lambda_{i}\langle
f,e_{i}\rangle_{H}e_{i}(x)$

for $f\in D(\Delta)$. We assume Dirichlet boundary conditions, so we call the
$\lambda_{i}$ Dirichlet eigenvalues and the $e_{i}$ Dirichlet eigenfunctions.
The nonlinear operator $F$.
In our setting, the operator $F$ on $H$ is a Nemytskii operator, i.e. there is
an associated function $u_{F}:\mathbb{R}\rightarrow\mathbb{R}$ such that the
operator $F$ can be defined as $[F(g)](x)=u_{F}(g(x))$ for $g\in H$. For $F$,
we adopt the assumption from [29, Assumption 2.4] and require that $F$ be two
times continuously Fréchet differentiable, and that its derivatives satisfy

$\displaystyle\|F^{\prime}(u)-F^{\prime}(v)\|$
$\displaystyle\leq\tilde{L}\|u-v\|,$ (3)
$\displaystyle\|(-\Delta)^{-r}F^{\prime}(u)(-\Delta)^{r}w\|$
$\displaystyle\leq\tilde{L}\|w\|$ (4)

for all $u,v\in H$ and $w\in D((-\Delta)^{r})$ for $r\in\\{0,\ 1/2,\ 1\\}$ and

$\displaystyle\|\Delta^{-1}F^{\prime\prime}(u)(v,w)\|\leq\tilde{L}\|(-\Delta)^{-1/2}v\|\cdot\|(-\Delta)^{-1/2}w\|$
(5)

for all $u,v,w\in H$ and $\tilde{L}>0$ a constant.
The Wiener process $W^{Q}$.
The process $W^{Q}=(W^{Q}(t))_{t\in[0,T]}$ is a $Q$-cylindrical Wiener process
[16, §4.1.2][37, §3.2.35], i.e. a collection

$\displaystyle\\{W_{u}^{Q}(t),\ u\in H,\ t\in[0,T]\\}$

of zero-mean Gaussian random variables that satisfy

$\displaystyle\mathbb{E}(W_{u}^{Q}(t)W_{v}^{Q}(s))=\langle
Qu,v\rangle_{H}\cdot\min(t,s)$

for $s,t\in[0,T]$, $u,v\in H$. Then, for an orthonormal basis
$\\{e^{Q}_{n}\\}_{n\in\mathbb{N}}$ of eigenfunctions of $Q$ with corresponding
eigenvalues $\\{q_{k}\\}_{k\in\mathbb{N}}$, the process $W^{Q}$ has the series
expansion in $H$ [16, Proposition 4.3]

$\displaystyle\textstyle
W^{Q}(t)=\sum_{n=1}^{\infty}\sqrt{q_{n}}\beta_{n}(t)e_{n}^{Q},$ (6)

where $\beta_{i}$, $i\in\mathbb{N}$, are independent real-valued standard
Brownian motions. We assume that $Q$ and the operator $\Delta$ both have a
joint eigenbasis consisting of the same functions, i.e. we assume
$\\{e_{n}^{Q}\\}_{n\in\mathbb{N}}=\\{e_{n}\\}_{n\in\mathbb{N}}$ (for details
on this assumption, see Appendix A). The stochastic integral with respect to
the cylindrical $Q$-Wiener process is defined in a similar way to the
classical Itô integral (for details, see [35, 2.5]).
In the literature on SPDEs, at least one of the two constraints have to be
made: Either a one-dimensional spatial domain is considered, where white noise
as well as trace-class noise can be considered, or the covariance operator $Q$
cannot be allowed to be $Q=I$, i.e. space-time white noise cannot be
considered. This is due to the fact that a solution of an SPDE on a
multidimensional domain (i.e. $D\subset\mathbb{R}^{n}$, $n\geq 2$) with space-
time white noise (i.e. $W^{Q}$ is a cylindrical Wiener process with $Q=I$)
does not exist even for $F=0$, since vital regularity properties of the
stochastic convolution

$\displaystyle\textstyle
W_{A}^{Q}(t)=\int_{0}^{t}\mathrm{e}^{A(t-s)}\;\mathrm{d}W^{Q}(s),\quad
t\in[0,T],$

in (2) are violated. See [16, Example 5.7] [41, p. 329] for discussions of
these failures in spatial dimensions of two and higher. The deciding factor is
the asymptotic behavior of Dirichlet eigenvalues, which has been described by
Weyl’s law [23, 4.1]: If $0<\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\dots$
are the Dirichlet eigenvalues of the negative Laplacian for an arbitrary
domain $D\subset\mathbb{R}^{d}$, then

$\displaystyle\textstyle\lambda_{n}\propto\frac{4\pi^{2}}{(\omega_{d}\mu_{d}(D))^{2/d}}n^{2/d},\quad
n\rightarrow\infty,$ (7)

where $\mu_{d}(D)$ is the $d$-dimensional Lebesgue measure of $D$ and
$\omega_{d}$ is the volume of the $d$-dimensional unit ball. The regularity
condition [29, Assumption 2.2] for the exponential Euler scheme which we will
use requires

$\displaystyle\sum_{n=1}^{\infty}(\lambda_{n})^{2\gamma-1}q_{n}<\infty$ (8)

for some $\gamma\in(0,1)$, which for space-time white noise, i.e. $q_{n}=1$
for all $n\in\mathbb{N}$, is not fulfilled in two dimensions. Instead, we
assume that $q_{n}=\mathcal{O}(n^{-\alpha})$ for some $\alpha>0$, in which
case (8) holds if $\gamma<\alpha/2$. In our numerical experiments, we will use
$q_{n}=n^{-2+\epsilon_{q}}$ for an arbitrarily small $\epsilon_{q}>0$ so that
$\gamma$ can be arbitrarily close to $1$ (adding $\epsilon_{q}$ is mostly a
formal requirement which is justified in Appendix A).

##  3 Numerical scheme.

In this work, we combine two different numerical methods: The spectral
Galerkin exponential Euler scheme is used in two spatial dimensions with
eigenvalues and corresponding eigenfunctions which are numerically
approximated using a boundary element collocation method in combination with
the Beyn algorithm.
We use the spectral Galerkin method for approximation in space and the
exponential Euler method for approximation in time to numerically solve the
SPDE (1) (Section 3.1). The spectral Galerkin scheme relies on projecting the
values of the process from $H$ to the $N$-dimensional subspace
$\text{span}(e_{1},\dots,e_{N})$, where the $e_{i}$ are the Dirichlet
eigenfunctions of $\Delta$ (spatial approximation, Sections 3.1.1, 3.4.2), and
then discretizing the mild solution (2) in time (temporal approximation,
Section 3.1.2). If the SPDE has to be solved on a bounded domain
$D\subset\mathbb{R}^{2}$, the Dirichlet eigenfunctions $e_{i}$ (and also the
eigenvalues) for the domain $D$ are needed. They will be computed using a
boundary integral equation method (Section 3.2) with the Beyn integral
algorithm (Section 3.3) as described in [31, 32]. Finally, we make a few
comments on assembling an accurate orthonormal basis of Dirichlet
eigenfunctions on a given shape (Section 3.4).

###  3.1 The spectral Galerkin exponential Euler scheme.

We recall the spectral Galerkin method and the exponential Euler scheme.

####  3.1.1 Approximation in space.

We use a projection

$\displaystyle P_{N}:H$ $\displaystyle\rightarrow
H_{N}:=\text{span}(e_{1},\dots,e_{N})$ $\displaystyle
f=\sum_{i=1}^{\infty}\langle f,e_{i}\rangle_{H}e_{i}$
$\displaystyle\mapsto\sum_{i=1}^{N}\langle f,e_{i}\rangle_{H}e_{i}$

with the inner product

$\displaystyle\langle f,e_{i}\rangle_{H}=\int_{D}f(x)\cdot
e_{i}(x)\;\mathrm{d}x$

and denote $U^{N}:=P_{N}U$, $F_{N}:=P_{N}F$, $A_{N}:=P_{N}A$,
$u_{0}^{N}:=P_{N}u_{0}$ and

$\displaystyle
W^{Q}_{N}:=P_{N}W^{Q}=\sum_{i=1}^{N}\sqrt{q_{i}}e_{i}\beta_{i}(t).$

The projected SPDE

$\displaystyle\mathrm{d}U^{N}(t)=(A_{N}U^{N}(t)+F_{N}(U^{N}(t)))\mathrm{d}t+\mathrm{d}W_{N}^{Q}(t)$

has the mild solution satisfying for $t\in[0,T]$

$\displaystyle
U^{N}(t)=\mathrm{e}^{A_{N}t}u_{0}^{N}+\int_{0}^{t}\mathrm{e}^{A_{N}(t-s)}F_{N}(U^{N}(s))\;\mathrm{d}s+\int_{0}^{t}\mathrm{e}^{A_{N}(t-s)}\;\mathrm{d}W_{N}^{Q}(s).$
(9)

####  3.1.2 Approximation in time.

We introduce the equidistant time steps $t_{i}=h\cdot i$, $i=0,\dots,M$, with
$h=T/M$, and denote $u_{k}^{N,M}=U^{N}(t_{k})$. Equation (9) can be rewritten
in one time step from $t_{k}$ to $t_{k+1}$ as

$\displaystyle
u_{k+1}^{N,M}=\mathrm{e}^{A_{N}h}u_{k}^{N,M}+\int_{t_{k}}^{t_{k+1}}\mathrm{e}^{A_{N}(t_{k+1}-s)}F_{N}(u^{N,M}(s))\mathrm{d}s+\int_{t_{k}}^{t_{k+1}}\mathrm{e}^{A_{N}(t_{k+1}-s)}\mathrm{d}W_{N}^{Q}(s).$

We now make the approximation $F_{N}(u^{N,M}(s))\approx
F_{N}(u^{N,M}(t_{k}))$. Let $\widetilde{U^{N}}$ be the process after this
approximation and let $v_{k}^{N,M}:=\widetilde{U^{N}}(t_{k})$. We denote the
inner products by

$\displaystyle v_{k,j}^{N,M}:=\langle v_{k}^{N,M},e_{j}\rangle_{H},\quad
F_{N}^{j}(v_{k}^{N,M}):=\langle F_{N}(v_{k}^{N,M}),e_{j}\rangle_{H}.$

Taking the inner product with $e_{j}$ on both sides and applying the above
approximation, we obtain for the approximate Fourier-Galerkin coefficients

$\displaystyle v^{N,M}_{k+1,j}$
$\displaystyle=\langle\mathrm{e}^{A_{N}h}v_{k}^{N,M},e_{j}\rangle+\langle\int_{t_{k}}^{t_{k+1}}\mathrm{e}^{A_{N}(t_{k+1}-s)}F_{N}(v_{k}^{N,M})\;\mathrm{d}s,e_{j}\rangle$
$\displaystyle\quad+\langle\int_{t_{k}}^{t_{k+1}}\mathrm{e}^{A_{N}(t_{k+1}-s)}\;\mathrm{d}W_{N}^{Q}(s),e_{j}\rangle$
$\displaystyle=\mathrm{e}^{-\lambda_{j}h}v_{k,j}^{N,M}+F_{N}^{j}(v_{k}^{N,M})\cdot\int_{t_{k}}^{t_{k+1}}\mathrm{e}^{-\lambda_{j}(t_{k+1}-s)}\;\mathrm{d}s$
$\displaystyle\quad+\int_{t_{k}}^{t_{k+1}}\mathrm{e}^{-\lambda_{j}(t_{k+1}-s)}\sqrt{q_{j}}\;\mathrm{d}\beta_{j}(s),$

using linearity of the inner product and the fact that
$A_{N}he_{i}=-\lambda_{i}he_{i}$ and hence
$\mathrm{e}^{A_{N}h}e_{i}=\mathrm{e}^{-\lambda_{i}h}e_{i}$ for $h>0$. The
$\mathrm{d}\beta_{j}(s)$ integral is normally distributed with mean zero and
its variance can be computed using the Itô isometry. We finally obtain

$\displaystyle
v_{k+1,j}^{N,M}=\mathrm{e}^{-\lambda_{j}h}v_{k,j}^{N,M}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}F_{N}^{j}(v_{k}^{N,M})+\left(\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})\right)^{1/2}R_{k}^{j},$
(10)

where the $R_{k}^{j}$ are independent identically distributed standard
normally distributed random variables, $j=1,\dots,N$, $k=0,\dots,M-1$. This is
the exponential Euler scheme, as described in [29, Section 3].

###  3.2 The boundary integral equation method.

The boundary integral equation method aims to solve boundary value problems
such as the Helmholtz equation with Dirichlet or Neumann boundary conditions
by writing the solution as an integral layer operator involving an unknown
density, approximating this density and subsequently reinserting it into the
integral to compute the solution.
Let $\Phi_{\kappa}(x,y)$ be the fundamental solution

$\displaystyle\Phi_{\kappa}(x,y)=\text{i}H_{0}^{(1)}(\kappa\|x-y\|)/4,\quad
x\neq y$ (11)

of the Helmholtz equation

$\displaystyle-\Delta u+\kappa^{2}u$ $\displaystyle=0\;\;\text{in}\ D,\qquad
u=0\;\;\text{on}\ \partial D,$ (12)

where $H_{0}^{(1)}$ denotes the first-kind Hankel function of order zero (see
[14, p. 66]) and let $\kappa\in\mathbb{C}\backslash\\{0\\}$ be fixed. For a
point $y$ on the boundary $\partial D$, we denote the exterior normal by
$\nu(y)$ and we use the double layer potential

$\displaystyle\text{DL}_{\kappa}[\psi](x)=\int_{\Gamma}\psi(y)\cdot\partial_{\nu(y)}\Phi_{\kappa}(x,y)\;\mathrm{d}s(y),\quad
x\in\mathbb{R}^{2}\backslash\partial D,$ (13)

where $\psi\in C(\partial D)$ is the unknown density function. We also define

$\displaystyle\text{D}_{\kappa}[\psi](x)=\int_{\Gamma}\psi(y)\cdot\partial_{\nu(y)}\Phi_{\kappa}(x,y)\;\mathrm{d}s(y),\quad
x\in\partial D.$

For $x\in\partial D$, it is known [14, Theorem 3.1] that the so-called jump
relations

$\displaystyle\begin{split}\lim_{h\rightarrow
0}&\int_{\Gamma}\psi(y)\partial_{\nu(y)}\Phi_{\kappa}(x\pm
h\nu(x),y)\;\mathrm{d}s(y)=\\\
&\int_{\Gamma}\psi(y)\partial_{\nu(y)}\Phi_{\kappa}(x,y)\;\mathrm{d}s(y)\pm\frac{1}{2}\psi(x)\end{split}$
(14)

hold. Since the fundamental solution $\Phi_{\kappa}$ satisfies the Helmholtz
equation (12), so does the double-layer potential (13). Conversely, every
solution of (12) can be written as a double-layer potential [14, p. 39]. We
use the double layer ansatz

$\displaystyle u(x)=\text{DL}_{\kappa}[\psi](x),\quad x\in D$ (15)

which, after letting $x$ tend towards the boundary $\partial D$, with the jump
relation (14) becomes

$\displaystyle-\frac{1}{2}\psi(x)+\text{D}_{\kappa}[\psi](x)=0,\quad
x\in\partial D.$ (16)

Writing $\mathrm{M}(\kappa):=-\frac{1}{2}\mathrm{I}+\mathrm{D}_{\kappa}$, (16)
becomes $\mathrm{M}(\kappa)\psi=0$, which is a nonlinear eigenvalue problem.
First, the boundary is subdivided into $n_{f}$ curved elements called
$\Delta_{1},\dots,\Delta_{n_{f}}$, and (16) becomes

$\displaystyle-\frac{1}{2}\psi(x)+\sum_{j=1}^{n_{f}}\int_{\Delta_{j}}\psi(y)\partial_{\nu(y)}\Phi_{\kappa}(x,y)\;\mathrm{d}s(y)=0,\quad
x\in\partial D.$ (17)

Let $m_{j}:[0,1]\rightarrow\Delta_{j}$, $j=1,\dots,n_{f}$, be bijective,
continuous maps. With a change of variables, (17) becomes

$\displaystyle-\frac{1}{2}\psi(x)+\sum_{j=1}^{n_{f}}\int_{0}^{1}\psi(m_{j}(s))\partial_{\nu(m_{j}(s))}\Phi_{\kappa}(x,m_{j}(s))\|\partial_{s}m_{j}(s)\|\;\mathrm{d}s(s)=0,\quad
x\in\partial D.$ (18)

We now denote by $L_{i}(s)$ the quadratic Lagrange polynomials and by
$\widehat{L_{i}}(s)$ the generalized Lagrange basis functions, $i=1,2,3$:

$\displaystyle L_{1}(s)=(1-s)(1-2s),\qquad$ $\displaystyle
L_{2}(s)=4s(1-s),\qquad\qquad\;\;\;$ $\displaystyle L_{3}(s)=s(2s-1),$
$\displaystyle\widehat{L_{1}}(s)=\frac{1-s-\alpha}{1-2\alpha}\frac{1-2s}{1-2\alpha},$
$\displaystyle\widehat{L_{2}}(s)=4\frac{s-\alpha}{1-2\alpha}\frac{1-s-\alpha}{1-2\alpha},$
$\displaystyle\widehat{L_{3}}(s)=\frac{s\alpha}{1-2\alpha}\frac{2s-1}{1-2\alpha}$

and we fix a parameter

$\displaystyle\alpha\in(0,1/2).$ (19)

Each $m_{j}(s)$ is approximated by a quadratic interpolation polynomial

$\displaystyle
m_{j}(s)\approx\tilde{m}_{j}(s)=\sum_{i=1}^{3}v^{(j)}_{i}L_{i}(s),$

where $v_{i}^{(j)}$, $i=1,3$ are the end points of the $j$-th boundary element
and $v_{2}^{(j)}$ is its midpoint. We define the collocation nodes
$\tilde{v}_{j,k}:=\tilde{m}_{j}(q_{k})$ for $j=1,\dots,n_{f}$, $k=1,2,3$ and
$q_{1}=\alpha,\ q_{2}=1/2,\ q_{3}=1-\alpha$. The unknown function
$\psi(\tilde{m}_{j}(s))$ is approximated on each piece $\Delta_{j}$ by a
quadratic interpolation polynomial

$\displaystyle\sum_{k=1}^{3}\psi(\tilde{m}_{j}(q_{k}))\widehat{L_{k}}(s)=\sum_{k=1}^{3}\psi(\tilde{v}_{j,k})\widehat{L_{k}}(s).$

Equation (18) then becomes

$\displaystyle-\frac{1}{2}\psi(x)+\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}\int_{0}^{1}\partial_{\nu(\tilde{m}_{j}(s))}\Phi_{\kappa}(x,\tilde{m}_{j}(s))\|\partial_{s}\tilde{m}_{j}(s)\|\widehat{L_{k}}(s)\;\mathrm{d}s(s)\psi(\tilde{v}_{j,k})=r(x),$

with a residue $r(x)$ which results from the various approximations. For the
boundary element collocation method, we require $r(\tilde{v}_{i,\ell})=0$ [3,
p. 50] and obtain the linear system of size $3n_{f}\times 3n_{f}$

$\displaystyle-\frac{1}{2}\psi(\tilde{v}_{i,\ell})+\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}a_{i,j,k,\ell}\psi(\tilde{v}_{j,k})=0$
(20)

with $i=1,\dots,n_{f}$, $\ell=1,2,3$ and with the integrals

$\displaystyle
a_{i,j,k,\ell}=\int_{0}^{1}\partial_{\nu(\tilde{m}_{j}(s))}\Phi_{\kappa}(\tilde{v}_{i,\ell},\tilde{m}_{j}(s))\|\partial_{s}\tilde{m}_{j}(s)\|\widehat{L_{k}}(s)\;\mathrm{d}s(s)$

which are approximated using a Gauß-Kronrod quadrature. We write (20)
abstractly as the nonlinear eigenvalue problem $M(\kappa)\phi=0$, where

$\displaystyle\phi=(\psi(v_{1,1}),\psi(v_{1,2}),\psi(v_{1,3}),\psi(v_{2,1}),\dots,\psi(v_{n_{f},3}))^{\top}.$
(21)

###  3.3 The Beyn integral algorithm.

For a given wavenumber $\kappa$, the solution of (12) can be found by
inserting the obtained density $\psi$ back into (15), so the eigenfunction
value can be computed at any given point inside $D$. But first, the nonlinear
eigenvalue problem $M(\kappa)\phi=0$ for the density must be solved. This is
accomplished with the Beyn algorithm [5, p. 3849], which makes use of
Keldysh’s theorem and is based on complex contour integrals of the resolvent
$M^{-1}(\kappa)$, whose poles are the Dirichlet eigenvalues.

######  Theorem 3.1 (Keldysh’s Theorem.).

Let $D\subset\mathbb{C}$ be a domain, and let $v^{*}=\overline{v^{\top}}$
denote the conjugate transpose of a vector or matrix. Suppose that $C\subset
D$ is a compact subset containing only simple eigenvalues $\lambda_{n}$,
$n=1,\dots,n(\gamma)$ with eigenvectors $v_{n}$, $w_{n}$ satisfying

$\displaystyle M(\lambda_{n})v_{n}=0,\quad w_{n}^{*}M(\lambda_{n})=0,\quad
w_{n}^{*}M^{\prime}(\lambda_{n})v_{n}=1.$

Then there exists some neighborhood $U$ of $C$ in $D$ and a holomorphic
function $R:U\rightarrow\mathbb{C}^{3n_{f}\times 3n_{f}}$ such that

$\displaystyle
M(z)^{-1}=\sum_{n=1}^{n(\gamma)}\frac{1}{z-\lambda_{n}}v_{n}w_{n}^{*}+R(z),\quad
z\in U\backslash\\{\lambda_{1},\dots,\lambda_{n(\gamma)}\\}.$ (22)

Let $\sigma(M)$ denote the spectrum of $M$ and let $\gamma\subset D$ be a
contour which satisfies $\gamma\cap\sigma(M)=\emptyset$ and containing
$n(\gamma)$ eigenvalues $\lambda_{1},\dots,\lambda_{n(\gamma)}$. Let
$f:D\rightarrow\mathbb{C}$ be holomorphic. Then, applying the residue theorem
to (22) we have

$\displaystyle\frac{1}{2\pi\mathrm{i}}\int_{\gamma}f(z)M(z)^{-1}\;\mathrm{d}z=\sum_{n=1}^{n(\gamma)}f(\lambda_{n})v_{n}w_{n}^{*}.$

For a specified contour $\gamma(t)=\mu+R\cos(t)+\text{i}R\sin(t)$,
$\mu\in\mathbb{R}$, $t\in[0,2\pi)$, $R\in\mathbb{R}_{>0}$ and the number
$n(\gamma)$ of eigenvalues enclosed in $\gamma$ (with multiplicity), we pick a
random matrix $\hat{V}\in\mathbb{C}^{3n_{f}\times\ell}$ with
$3n_{f}\gg\ell\geq n(\gamma)$. We compute the contour integrals

$\displaystyle\textstyle A_{0}$
$\displaystyle=\frac{1}{2\pi\text{i}}\int_{\gamma}M^{-1}(\kappa)\hat{V}\;\mathrm{d}s(\kappa)$
$\displaystyle=\frac{1}{2\pi\text{i}}\int_{0}^{2\pi}M^{-1}(\gamma(t))\hat{{V}}\gamma^{\prime}(t)\;\mathrm{d}s(t)$
$\displaystyle=\sum_{n=1}^{n(\gamma)}v_{n}w_{n}^{*}\hat{V}=VW^{*}\hat{V},$
$\displaystyle A_{1}$ $\displaystyle=\frac{1}{2\pi\text{i}}\int_{\gamma}\kappa
M^{-1}(\kappa)\hat{V}\;\mathrm{d}s(\kappa)$
$\displaystyle=\frac{1}{2\pi\text{i}}\int_{0}^{2\pi}\gamma(t)M^{-1}(\gamma(t))\hat{{V}}\gamma^{\prime}(t)\;\mathrm{d}s(t)$
$\displaystyle=\sum_{n=1}^{n(\gamma)}\lambda_{n}v_{n}w_{n}^{*}\hat{V}=V\Lambda
W^{*}\hat{V},$

with $\Lambda=\text{diag}(\lambda_{1},\dots,\lambda_{n(\gamma)})$. These
integrals are approximated using the trapezoidal rule:

$\displaystyle A_{0,K}$
$\displaystyle=\frac{1}{\text{i}K}\sum_{j=0}^{K-1}M^{-1}(\gamma(t_{j}))\hat{V}\gamma^{\prime}(t_{j}),\;A_{1,K}=\frac{1}{\text{i}K}\sum_{j=0}^{K-1}\gamma(t_{j})M^{-1}(\gamma(t_{j}))\hat{V}\gamma^{\prime}(t_{j}),$
(23)

with $t_{j}=\frac{2\pi j}{K},\ j=0,\dots,K$. Next, we compute a singular value
decomposition $A_{0,K}=V\Sigma W^{*}$ with
$V\in\mathbb{C}^{3n_{f}\times\ell}$, $\Sigma\in\mathbb{C}^{\ell\times\ell}$
and $W\in\mathbb{C}^{\ell\times\ell}$. For
$\Sigma=\text{diag}(\sigma_{1},\sigma_{2},\dots,\sigma_{\ell})$ and a given
tolerance $\epsilon_{\mathrm{sing}}$, we find $n(\gamma)$ such that

$\displaystyle\sigma_{1}\geq\dots\geq\sigma_{n(\gamma)}>\epsilon_{\mathrm{sing}}>\sigma_{n(\gamma)+1}\geq\dots\geq\sigma_{\ell}.$
(24)

The singular values $\sigma_{n(\gamma)+1}\geq\dots\geq\sigma_{\ell}$ which are
below the tolerance are likely to be different from zero only due to rounding
errors and are discarded. With the matrices $V_{0}:=(V_{ij})_{1\leq i\leq
m,1\leq j\leq n(\gamma)},\ \Sigma_{0}:=(\Sigma_{ij})_{1\leq i\leq
n(\gamma),1\leq j\leq n(\gamma)}$ as well as $W_{0}:=(W_{ij})_{1\leq
i\leq\ell,1\leq j\leq n(\gamma)}$ there is a regular matrix
$S\in\mathbb{C}^{n(\gamma)\times n(\gamma)}$ such that $V=V_{0}S$, and after
rewriting some matrices, we have that

$\displaystyle S\Lambda S^{-1}=V_{0}^{*}A_{1}W_{0}\Sigma_{0}^{-1}.$

The matrix on the right-hand side therefore has the same eigenvalues and
eigenvectors as $\Lambda$. We compute the $n(\gamma)$ eigenvalues and
eigenvectors $s_{i}$ of the matrix
$V_{0}^{*}A_{1,N}W_{0}\Sigma_{0}^{-1}\in\mathbb{C}^{n(\gamma)\times
n(\gamma)}$. The $i$-th eigenfunction $u_{i}$ is then approximated by
inserting $\phi=V_{0}s_{i}$ (see (21)) into

$\displaystyle
u_{i}(x)=\mathrm{DL}_{\kappa}[\psi](x)\approx\sum_{j=1}^{n_{f}}\sum_{k=1}^{3}\hat{a}_{j,k}(x)\cdot\psi(\tilde{v}_{j,k}),$

where

$\displaystyle\hat{a}_{j,k}(x)=\int_{0}^{1}\partial_{\nu(\tilde{m}_{j}(s))}\Phi_{\kappa}(x,\tilde{m}_{j}(s))\|\partial_{s}\tilde{m}_{j}(s)\|\widehat{L_{k}}(s)\;\mathrm{d}s(s),$

which again are approximated using a Gauß-Kronrod quadrature.

$0$$10$$20$$30$$40$$50$$60$$70$$80$$100$$150$$200$$250$$300$$350$$400$$450$$500$$\kappa$$n_{f}$$\epsilon^{(\lambda)}=2\cdotp
10^{-4}$$\epsilon^{(\eta)}=5\cdotp 10^{-6}$
$0$$10$$20$$30$$40$$50$$60$$70$$80$$100$$150$$200$$250$$300$$350$$400$$450$$500$$\kappa$$n_{f}$$\epsilon^{(\lambda)}=2\cdotp
10^{-4}$$\epsilon^{(\eta)}=5\cdotp 10^{-6}$
---
$0$$10$$20$$30$$40$$50$$60$$70$$80$$100$$150$$200$$250$$300$$350$$400$$450$$500$$\kappa$$n_{f}$$\epsilon^{(\lambda)}=1\cdotp
10^{-4}$$\epsilon^{(\eta)}=5\cdotp 10^{-6}$
$0$$10$$20$$30$$40$$50$$60$$70$$80$$100$$150$$200$$250$$300$$350$$400$$450$$500$$\kappa$$n_{f}$$\epsilon^{(\lambda)}=1\cdotp
10^{-4}$$\epsilon^{(\eta)}=5\cdotp 10^{-6}$
$0$$10$$20$$30$$40$$50$$60$$70$$80$$100$$150$$200$$250$$300$$350$$400$$450$$500$$\kappa$$n_{f}$$\epsilon^{(\lambda)}=2\cdotp
10^{-4}$$\epsilon^{(\eta)}=2\cdotp 10^{-6}$
$0$$10$$20$$30$$40$$50$$60$$70$$80$$100$$150$$200$$250$$300$$350$$400$$450$$500$$\kappa$$n_{f}$$\epsilon^{(\lambda)}=2\cdotp
10^{-4}$$\epsilon^{(\eta)}=2\cdotp 10^{-6}$
Figure 1: The number of boundary elements to approximate wavenumbers and their
respective eigenfunctions to stay below certain error tolerances
$\epsilon^{(\lambda)}$ and $\epsilon^{(\eta)}$. For each base function, we
computed a reference solution with $n_{f}^{\mathrm{ref}}=600$ boundary
elements and compared the same function with less boundary elements to it,
with $n_{f}\in\\{50,75,100,150,200,300,400,500\\}$ boundary elements used. The
plots on the left show raw data, with the least number of required boundary
elements. The plots on the right show linearly interpolated numbers of needed
boundary elements based off the data from the plots on the left. We can
observe the linear relationship between wavenumber and needed boundary
elements $n_{f}$, and that the needed $n_{f}$ increases as the error
tolerances decrease (in red are the values which exceed the error tolerance
for every tested $n_{f}$). The blue reference slopes are functions
$f(x)=mx+b$, for $m\in\\{7,8.5,10\\}$, $b\in\\{80,65,70\\}$, from top to
bottom.

###  3.4 Preprocessing: Building an orthonormal basis (ONB)

The analytical expressions of the Dirichlet eigenfunctions are known for some
shapes. In the following, we will consider a shape which we call the Peanut
shape (see Figure 3), whose Dirichlet eigenvalues and eigenfunctions are not
known analytically and whose parametrization for $t\in[0,2\pi)$ is given by

$\begin{split}\textstyle\left(\begin{array}[]{c}0.06\cdot((\cos(t)+2)\cdot(\cos(t+0.6)+2)\cdot(0.1\cdot\cos(3t)+2))-0.1\\\
0.06\cdot(\sin(t)+2)\cdot(\sin(t-0.5)+2)\cdot(0.4\cdot\cos(2t)+2)\cdot(0.1\cdot\sin(4t)+1))-0.06\end{array}\right).\end{split}$

For the boundary element collocation method, we choose
$\alpha=(1-\sqrt{3/5})/2$ (see (19)). Note that while any
$\alpha\in(0,\frac{1}{2})$ ensures at least cubic convergence, empirically it
was seen [33, Section 4.6] that the choice $\alpha=(1-\sqrt{3/5})/2$ yielded
superconvergence for domains with a $\mathcal{C}^{2}$ boundary. Whenever the
Beyn integral algorithm is used to compute its Dirichlet eigenvalues and
eigenfunctions, we use $\epsilon_{\mathrm{sing}}=10^{-4}$ as singular value
tolerance (see (24)) and $K=24$ for the periodic trapezoidal rule (see (23)).

####  3.4.1 Convergence and accuracy of base functions.

In previous work, the Beyn integral algorithm was used to compute Neumann
eigenvalues, and a convergence with rates up to three for a fixed wavenumber
could be observed [33]. However, the eigenvalues and eigenfunctions become
less accurate for higher wavenumbers, so that more boundary elements are
needed for higher wavenumbers. Previous publications from high-frequency
scattering theory using boundary element collocation methods to solve the
Helmholtz equation have pointed out that this effort rises linearly with the
wavenumber [11, 12], which we could confirm for the case of the Dirichlet
eigenfunctions and eigenvalues (see Figure 1). For the Peanut shape, we have
conducted a convergence analysis for the first 200 eigenfunctions. For each
wavenumber, we computed a reference function for $600$ boundary elements which
we compare to functions with lower accuracy (computed with 500, 400, 300, 200,
150, 100, 75 and 50 boundary elements, respectively). For speeding up the
computations, these functions were stored only on an equidistant $81\times 81$
grid inside $[0,1]^{2}$, while the final functions making up the ONB are
stored on a $301\times 301$ grid to obtain results with better resolution (for
the error analysis in Section 4, we will denote the grid resolution as $R$ and
will use $R=301$). The pointwise storage on the equidistant grid means that
the inner products, i.e. the integrals in the scheme (10) can be approximated
numerically. As it is easy to implement and has a convergence order of four
[21, p. 21 f.] (and hence, a convergence order of two for two-dimensional
functions), we use the composite two-dimensional Simpson rule and set the
numerically integrated function to be zero outside the boundary, so that the
Simpson rule can be applied on the unit square. The application of the Simpson
rule is possible since eigenfunctions of the Laplacian are always smooth [23,
§2.(i)].
Suppose that we use $n_{f}^{(i)}$ boundary elements for eigenfunction number
$i$, $i=1,\dots,N$. We fix error tolerances

$\displaystyle\epsilon^{(\lambda)}:=\max_{i=1,\dots,N}|\widetilde{\lambda_{i}}^{(n_{f}^{(i)})}-\tilde{\lambda}_{i}^{(600)}|,\
\epsilon^{(\eta)}$
$\displaystyle:=\max_{i=1,\dots,N}\|\tilde{e}_{i}^{(n_{f}^{(i)})}-\tilde{e}_{i}^{(600)}\|_{H},$

where the tilde denotes approximation and the superscript denotes that the
approximation is done with $n_{f}^{(i)}$ boundary elements. Depending on which
$\epsilon^{(\lambda)}$ and $\epsilon^{(\eta)}$ we choose, different wavenumber
dependent amounts of boundary elements are needed (see Figure 1). We will
later observe (see Example 5.1) that for the error bound derived in Theorem
4.4 it is crucial to reduce $\epsilon^{(\lambda)}$ and $\epsilon^{(\eta)}$ as
much as possible.

####  3.4.2 Computation of the orthonormal basis.

After precomputing the necessary numbers of boundary elements for each
function (see Section 3.4.1) and the spectrum of $-\Delta$ on $D$ is known up
to a few digits, the eigenvalues and eigenfunctions can be computed by
applying the Beyn algorithm to each eigenvalue by enclosing it by a contour
which does not enclose any other eigenvalue (see Figure 2). The Peanut shape
analyzed here is asymmetric and does not appear to have any multiple Dirichlet
eigenvalues. If multiple eigenvalues occur, it can happen that the different
eigenfunctions belonging to one eigenvalue are not returned by the algorithm
as orthogonal functions, but as linear combinations of them. This would
necessitate an additional orthogonalization routine applied to those
eigenfunctions, which would pose another error source.
Based on the convergence results for the Peanut shape, our base functions for
the Peanut have been computed such that the error to a reference function with
$n_{f}^{\mathrm{ref}}=600$ boundary elements stays below the tolerances
$\epsilon^{(\lambda)}=2\cdot 10^{-4}$ and $\epsilon^{(\eta)}=5\cdot 10^{-6}$.

$\mathbb{R}$$\mathbb{C}$$0$$\lambda_{1}$$\lambda_{2}$$\lambda_{3}$$\dots$$\times$$\times$$\times$$\hat{\gamma}_{1}$$\hat{\gamma}_{2}$$\hat{\gamma}_{3}$$\hat{\gamma}_{4}$$\mathbb{R}$$\mathbb{C}$$0$$\lambda_{1}$$\lambda_{2}$$\lambda_{3}$$\dots$$\times$$\times$$\times$$\gamma_{1}$$\gamma_{2}$$\gamma_{3}$
Figure 2: In the first step, the eigenvalues are still unknown, so the real
axis up to a certain point is covered with contours $\hat{\gamma}_{i}$ which
cover the real axis up to the point where the desired number $N$ of
eigenvalues has been found. For the subsequent, more precise computation,
contours $\gamma_{i}$, $i=1,\dots,N$ around the found eigenvalues are chosen.

In addition, we have computed some errors to even more accurate wavenumbers
and eigenfunctions computed with $n_{f}^{\mathrm{ref}}=1,200$ boundary
elements. A comparison of the reference errors is shown in Figure 3 for two
wavenumbers. It can be seen that there is a slight discrepancy for the
wavenumber error at $n_{f}=500$, but the reference error for the
eigenfunctions is almost the same for both reference functions, so it can be
assumed that the reference functions are sufficiently accurate.
On the computational cost, we note that for an orthonormal basis of a certain
number $N$ of eigenvalues and eigenfunctions the computational time with our
implementation behaved asymptotically as $\mathcal{O}(N^{2})$. In more detail,
it is shown in Figure 4 that, with growing number $n_{f}$ of used boundary
elements, the computation times for an eigenvalue rises as
$\mathcal{O}(n_{f}^{2})$, and for an eigenfunction as $\mathcal{O}(n_{f})$
(although it is more costly for low $n_{f}$). The computation time is very
similar for different eigenvalues and the corresponding eigenfunctions, and
neither increases nor decreases as the eigenvalue increases. Figure 5 shows an
estimate of the total computational cost to compute an orthonormal basis of
$N$ eigenvalues and eigenfunctions, i.e. the sum

$\displaystyle\sum_{i=1}^{N}\mathrm{ev}(n_{f}(\lambda_{i}))+\mathrm{ef}(n_{f}(\lambda_{i})),$
(25)

where $n_{f}(\lambda)$ is the needed number of boundary elements for a given
set of error tolerances, and $\mathrm{ev}(n_{f})$ and $\mathrm{ef}(n_{f})$ are
the computation times for eigenvalue and eigenfunction using $n_{f}$ boundary
elements. Based on empirical tests such as shown in Figure 4, we assumed that
$\mathrm{ev}(n_{f})=0.15\cdot n_{f}^{2}$ and $\mathrm{ef}(n_{f})=50\cdot
n_{f}$. Since $\lambda_{i}\sim i$, $n_{f}(i)\sim\sqrt{i}$ (see 1),
$\mathrm{ev}(i)\sim i^{2}$ and $\mathrm{ef}(i)\sim i$, the sum (25) contains a
term behaving as $i$ and one behaving as $\sqrt{i}$, so that the whole sum is
made up of two terms (one each for eigenvalues and eigenfunctions) behaving as
$\mathcal{O}(N^{2})+\mathcal{O}(N^{3/2})$.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$x$$y$
$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$x$$y$
---
$0$$200$$400$$600$$800$$1{,}000$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$n_{f}$$|\tilde{\lambda}_{1}^{(n_{f})}-\tilde{\lambda}_{1}^{(n_{f}^{\mathrm{ref}})}|$
$0$$200$$400$$600$$800$$1{,}000$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$n_{f}$$|\tilde{\lambda}_{56}^{(n_{f})}-\tilde{\lambda}_{56}^{(n_{f}^{\mathrm{ref}})}|$
$0$$200$$400$$600$$800$$1{,}000$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$n_{f}$$\|\tilde{e}_{1}^{(n_{f})}-\tilde{e}_{1}^{(n_{f}^{\mathrm{ref}})}\|$
$0$$200$$400$$600$$800$$1{,}000$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$n_{f}$$\|\tilde{e}_{56}^{(n_{f})}-\tilde{e}_{56}^{(n_{f}^{\mathrm{ref}})}\|$$\;\;n_{f}^{\mathrm{ref}}=1200$$n_{f}^{\mathrm{ref}}=600$
Figure 3: Convergence results for two eigenfunctions of the Peanut shape (on
the left, the first one, on the right, the 56th one), in red, the error is
measured with respect to reference eigenvalues and eigenfunctions computed for
$n_{f}^{\mathrm{ref}}=600$ boundary elements, in green, the reference data is
computed for $n_{f}^{\mathrm{ref}}=1200$ boundary elements.

##  4 Error analysis.

$100$$200$$300$$400$$500$$600$$10^{-1}$$10^{-2}$$10^{-3}$$10^{-4}$$10^{-5}$$n_{f}$Absolute
errors$10$$100$$200$$400$$600$$2$$4$$6$$\cdot 10^{4}$$n_{f}$Computing time
(s)EigenvalueEigenfunction
---
Figure 4: This plot shows, for different numbers of boundary elements $n_{f}$,
absolute errors for the 300th eigenvalue and its corresponding eigenfunction
in green (left axis) as well as CPU time for their computation in seconds
(right axis). It can be seen that the CPU time for the eigenvalues behaves
like $\mathcal{O}(n_{f}^{2})$, while the time for the eigenfunctions is
$\mathcal{O}(n_{f})$. The computing times are very similar for other
eigenvalues and eigenfunctions. $0$$100$$200$$300$$400$$500$$2$$4$$6$$\cdot
10^{7}$$N$Total CPU time (s)Scenario 3Scenario 2Scenario 1
---
Figure 5: The total computation time on our machine for an orthonormal basis
consisting of $N$ eigenvalues and eigenfunctions. Scenarios 1 to 3 correspond,
from top to bottom, to the three pairs of error tolerances shown in Figure 1.

We present an error analysis for the strong error in the Galerkin exponential
Euler scheme from [29] with approximated eigenfunctions and eigenvalues. We
denote $\langle\cdot,\cdot\rangle:=\langle\cdot,\cdot\rangle_{H}$ and
$\|\cdot\|:=\|\cdot\|_{H}$, and the exact solution at time $t$ by $U(t)$. Its
approximation by the spectral Galerkin exponential Euler scheme is denoted by
$V_{k}^{(1)}$. The approximation of $V_{k}^{(1)}$ by computing the inner
products as numerical integrals over $D$ is denoted by $V_{k}^{(2)}$. The
approximation of $V_{k}^{(2)}$ using approximate eigenvalues is denoted by
$V_{k}^{(3)}$. Finally, the approximation of $V_{k}^{(3)}$ using approximate
eigenfunctions is denoted by $V_{k}^{(4)}$. We denote

$\displaystyle V_{k,j}^{(i)}:=\langle V_{k}^{(i)},e_{j}\rangle,\quad
i=1,2,3,4,\quad\epsilon_{k,j}^{(i)}:=\mathbb{E}|V_{k,j}^{(i-1)}-V_{k,j}^{(i)}|,\quad
i=2,3,4,$

with $k=0,\dots,M$ and $j=1,\dots,N$. First, error estimates for
$\|V_{k}^{(i)}-V_{k}^{(i+1)}\|$ are given for $i=1,2,3$ in Lemmas 4.1, 4.2 and
4.3 which are then combined to give a general result on the error
$\|U(t_{k})-V_{k}^{(4)}\|$ in Theorem 4.4. Note that $V_{k}^{(4)}$ is the
output of the computer programs, so the error estimate in Theorem 4.4 contains
Fourier-Galerkin coefficients of $V_{k}^{(4)}$ which can then either be
inserted for a run-specific bound, or be bound by an empirically found
constant which is generally found to be an upper bound (which we have done in
Example 5.1).

######  Lemma 4.1.

We have

$\displaystyle\mathbb{E}\|V_{k}^{(1)}-V_{k}^{(2)}\|\leq\sum_{j=1}^{N}\epsilon_{k,j}^{(2)}$
(26)

and, for a constant $C>0$,

$\displaystyle\epsilon_{k+1,j}^{(2)}\leq\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\cdot\left[L\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(2)}+C\cdot
R^{-2}\right].$ (27)

Proof. We compute, distinguishing $f_{N}^{j}(V_{k}^{(2)})$ computed as an
integral and $\widetilde{f_{N}^{j}(V_{k}^{(2)})}$ computed as a numerical
integral,

$\displaystyle\begin{split}\mathbb{E}\|V_{k}^{(1)}-V_{k}^{(2)}\|&=\mathbb{E}\left\|\sum_{j=1}^{N}V_{k,j}^{(1)}\cdot
e_{j}-\sum_{j=1}^{N}V_{k,j}^{(2)}\cdot
e_{j}\right\|=\mathbb{E}\left\|\sum_{j=1}^{N}(V_{k,j}^{(1)}-V_{k,j}^{(2)})\cdot
e_{j}\right\|\\\
&\leq\sum_{j=1}^{N}\mathbb{E}|V_{k,j}^{(1)}-V_{k,j}^{(2)}|=\sum_{j=1}^{N}\epsilon_{k,j}^{(2)}\end{split}$
(28)

and

$\displaystyle\epsilon_{k+1,j}^{(2)}$
$\displaystyle=\mathbb{E}\left|\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(1)}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(1)})+\left(\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})\right)^{1/2}R_{k}^{j}\right.$
$\displaystyle\quad-\left.\left[\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\widetilde{f_{N}^{j}(V_{k}^{(2)})}+\left(\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})\right)^{1/2}R_{k}^{j}\right]\right|$
$\displaystyle\leq\
\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\mathbb{E}|f_{N}^{j}(V_{k}^{(1)})-\widetilde{f_{N}^{j}(V_{k}^{(2)})}|.$
$\displaystyle\mathbb{E}|f_{N}^{j}(V_{k}^{(1)})-\widetilde{f_{N}^{j}(V_{k}^{(2)})}|\leq\mathbb{E}|f_{N}^{j}(V_{k}^{(1)})-f_{N}^{j}(V_{k}^{(2)})|+\mathbb{E}|f_{N}^{j}(V_{k}^{(2)})-\widetilde{f_{N}^{j}(V_{k}^{(2)})}|$
$\displaystyle\leq\ $
$\displaystyle\mathbb{E}\int_{D}\left[f_{N}(V_{k}^{(1)}(x))-f_{N}(V_{k}^{(1)}(x)+\eta_{k}(x))\right]e_{j}(x)\;\mathrm{d}x+C\cdot
R^{-2}$ $\displaystyle\leq\ $
$\displaystyle\mathbb{E}\left[\int_{D}C[L\cdot\eta_{k}(x)]^{2}\;\mathrm{d}x\right]^{1/2}+C\cdot
R^{-2}=L\cdot\|\eta_{k}\|+C\cdot R^{-2}$ $\displaystyle\leq\ $ $\displaystyle
L\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(2)}+C\cdot R^{-2},$

where the Hölder inequality and the Lipschitz continuity of $f_{N}$ were used.
In total, we have

$\epsilon_{k+1,j}^{(2)}\leq\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\cdot\left[L\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(2)}+C\cdot
R^{-2}\right].\hfill\qed$

######  Lemma 4.2.

We have

$\displaystyle\mathbb{E}\|V_{k}^{(2)}-V_{k}^{(3)}\|\leq\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}$
(29)

and

$\displaystyle\epsilon_{k+1,j}^{(3)}$
$\displaystyle\leq\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(3)}+\mathrm{e}^{-\lambda_{j}h}|1-\mathrm{e}^{-\epsilon_{j}^{(\lambda)}h}|\cdot(\epsilon_{k,j}^{(4)}+\mathbb{E}|V_{k,j}^{(4)}|)$
(30)
$\displaystyle\quad+L\cdot\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}$
$\displaystyle\quad+\frac{1}{\lambda_{j}\tilde{\lambda}_{j}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{j}h}\lambda_{j}-\mathrm{e}^{-\lambda_{j}h}\tilde{\lambda}_{j}|\right]\cdot\left[L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right]$
$\displaystyle\quad+\left|\sqrt{\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})}-\sqrt{\frac{q_{j}}{2(\lambda_{j}+\epsilon_{j}^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{j}+\epsilon_{j}^{(\lambda)})h})}\right|$

Proof. Suppose that we have approximate eigenvalues
$\tilde{\lambda}_{j}=\lambda_{j}+\epsilon_{j}^{(\lambda)}$, $j=1,\dots,N$.
Then by the same steps as in (28), the resulting error of the whole function
is

$\displaystyle\mathbb{E}\|V_{k}^{(2)}-V_{k}^{(3)}\|\leq\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}.$

The coefficient errors after the next time step are given by

$\displaystyle\epsilon_{k+1,j}^{(3)}$
$\displaystyle=\mathbb{E}\left|\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(2)})+\left(\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-\lambda_{j}h})\right)^{1/2}R_{k}^{j}\right.$
$\displaystyle\quad-\left.\left[\mathrm{e}^{-\tilde{\lambda}_{j}h}V_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}f_{N}^{j}(V_{k}^{(3)})+\left(\frac{q_{j}}{2\tilde{\lambda}_{j}}(1-\mathrm{e}^{-2\tilde{\lambda}_{j}h})\right)^{1/2}R_{k}^{j}\right]\right|$
$\displaystyle\leq\mathbb{E}|\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(2)}-\mathrm{e}^{-\tilde{\lambda}_{j}h}V_{k,j}^{(3)}|$
$\displaystyle\quad+\mathbb{E}\left|\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(2)})-\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}f_{N}^{j}(V_{k}^{(3)})\right|$
$\displaystyle\quad+\mathbb{E}\left|\left[\left(\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})\right)^{1/2}-\left(\frac{q_{j}}{2\tilde{\lambda}_{j}}(1-\mathrm{e}^{-2\tilde{\lambda}_{j}h})\right)^{1/2}\right]R_{k}^{j}\right|.$

Going through this term by term, we have that

$\displaystyle\quad\mathbb{E}|\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(2)}-\mathrm{e}^{-\tilde{\lambda}_{j}h}V_{k,j}^{(3)}|$
$\displaystyle\leq\mathbb{E}|\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(2)}-\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(3)}|+\mathbb{E}|\mathrm{e}^{-\lambda_{j}h}V_{k,j}^{(3)}-\mathrm{e}^{-\tilde{\lambda}_{j}h}V_{k,j}^{(3)}|$
$\displaystyle=\mathrm{e}^{-\lambda_{j}h}\mathbb{E}|V_{k,j}^{(2)}-V_{k,j}^{(3)}|+|\mathrm{e}^{-\lambda_{j}h}-\mathrm{e}^{-(\lambda_{j}+\epsilon_{j}^{(\lambda)})h}|\cdot\mathbb{E}|V_{k,j}^{(3)}|$
$\displaystyle\leq\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(3)}+\mathrm{e}^{-\lambda_{j}h}|1-\mathrm{e}^{-\epsilon_{j}^{(\lambda)}h}|\cdot\mathbb{E}|V_{k,j}^{(3)}|$
$\displaystyle\leq\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(3)}+\mathrm{e}^{-\lambda_{j}h}|1-\mathrm{e}^{-\epsilon_{j}^{(\lambda)}h}|\cdot(\mathbb{E}|V_{k,j}^{(3)}-V_{k,j}^{(4)}|+\mathbb{E}|V_{k,j}^{(4)}|)$
$\displaystyle\leq\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(3)}+\mathrm{e}^{-\lambda_{j}h}|1-\mathrm{e}^{-\epsilon_{j}^{(\lambda)}h}|\cdot(\epsilon_{k,j}^{(4)}+\mathbb{E}|V_{k,j}^{(4)}|),$
(31)

where we have made the last estimate since unlike $\mathbb{E}[V_{k}^{(3)}]$,
$\mathbb{E}[V_{k}^{(4)}]$ can be approximated from the numerical
implementation. Next, we have

$\displaystyle\quad\mathbb{E}\left|\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(2)})-\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}f_{N}^{j}(V_{k}^{(3)})\right|$
$\displaystyle\leq\mathbb{E}\left|\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(2)})-\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(3)})\right|$
$\displaystyle\quad+\mathbb{E}\left|\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}f_{N}^{j}(V_{k}^{(3)})-\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}f_{N}^{j}(V_{k}^{(3)})\right|$
$\displaystyle=\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\mathbb{E}|f_{N}^{j}(V_{k}^{(2)})-f_{N}^{j}(V_{k}^{(3)})|+\left|\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}-\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\right|\cdot\mathbb{E}|f_{N}^{j}(V_{k}^{(3)})|$
$\displaystyle\leq
L\cdot\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}$
$\displaystyle\quad+\frac{1}{\lambda_{j}\tilde{\lambda}_{j}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{j}h}\lambda_{j}-\mathrm{e}^{-\lambda_{j}h}\tilde{\lambda}_{j}|\right]\cdot\left[\mathbb{E}|f_{N}^{j}(V_{k}^{(3)})-f_{N}^{j}(V_{k}^{(4)})|+\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right]$
(32) $\displaystyle\leq
L\cdot\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}$
(33)
$\displaystyle\quad+\frac{1}{\lambda_{j}\tilde{\lambda}_{j}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{j}h}\lambda_{j}-\mathrm{e}^{-\lambda_{j}h}\tilde{\lambda}_{j}|\right]\cdot\left[L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right].$

Finally, we have (using that $\mathbb{E}|X|\leq[\mathbb{E}|X|^{2}]^{1/2}$)

$\displaystyle\mathbb{E}\left|\left[\left(\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})\right)^{1/2}-\left(\frac{q_{j}}{2\tilde{\lambda}_{j}}(1-\mathrm{e}^{-2\tilde{\lambda}_{j}h})\right)^{1/2}\right]R_{k}^{j}\right|$
$\displaystyle=\left|\sqrt{\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})}-\sqrt{\frac{q_{j}}{2(\lambda_{j}+\epsilon_{j}^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{j}+\epsilon_{j}^{(\lambda)})h})}\right|\mathbb{E}|(R_{k}^{j})^{2}|^{1/2},$
(34)

and $\mathbb{E}|(R_{k}^{j})^{2}|=\text{Var}(R_{k}^{j})=1$. Summing (31), (33)
and (34), we obtain (30).∎

######  Lemma 4.3.

We have

$\displaystyle\mathbb{E}\|V_{k}^{(3)}-V_{k}^{(4)}\|\leq\sum_{j=1}^{N}\left[\epsilon_{k,j}^{(4)}+\epsilon^{(\eta)}\mathbb{E}|V_{k,j}^{(4)}|\right]$
(35)

and

$\displaystyle\begin{split}\epsilon_{k+1,j}^{(4)}&\leq\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{k,j}^{(4)}+\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{j}^{(\eta)}\mathbb{E}\|V_{k}^{(4)}\|\\\
&\quad+\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\left[L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\epsilon_{j}^{(\eta)}\cdot\mathbb{E}\|f_{N}(V_{k}^{(4)})\|\right].\end{split}$
(36)

Proof. Suppose that we have approximate eigenfunctions

$\displaystyle\tilde{e}_{j}=e_{j}+\eta_{j},\quad\|\eta_{j}\|\leq\epsilon_{j}^{(\eta)}.$

Then, the Fourier-Galerkin expansion of $V_{k}^{(4)}$ also involves the
approximate eigenfunctions, i.e.
$V_{k}^{(4)}=\sum_{j=1}^{N}(e_{j}+\eta_{j})V_{k,j}^{(4)}$, and the error for
the function is

$\displaystyle\mathbb{E}\|V_{k}^{(3)}-V_{k}^{(4)}\|$
$\displaystyle=\mathbb{E}\left\|\sum_{j=1}^{N}e_{j}V_{k,j}^{(3)}-\sum_{j=1}^{N}(e_{j}+\eta_{j})V_{k,j}^{(4)}\right\|$
$\displaystyle=\mathbb{E}\left\|\sum_{j=1}^{N}e_{j}(V_{k,j}^{(3)}-V_{k,j}^{(4)})-\sum_{j=1}^{N}\eta_{j}V_{k,j}^{(4)}\right\|$
$\displaystyle\leq\;\;\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\sum_{j=1}^{N}\epsilon^{(\eta)}\cdot\mathbb{E}|V_{k,j}^{(4)}|$
$\displaystyle=\;\;\sum_{j=1}^{N}\left[\epsilon_{k,j}^{(4)}+\epsilon^{(\eta)}\mathbb{E}|V_{k,j}^{(4)}|\right].$

The error for the next coefficient errors is (with the wide tilde denoting
inner products with $e_{j}+\eta_{j}$ instead of with $e_{j}$)

$\displaystyle\epsilon_{k+1,j}^{(4)}$
$\displaystyle=\mathbb{E}\left|\mathrm{e}^{-\tilde{\lambda}_{j}h}V_{k,j}^{(3)}+\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}f_{N}^{j}(V_{k}^{(3)})+\left(\frac{q_{j}}{2\tilde{\lambda}_{j}}(1-\mathrm{e}^{-2\tilde{\lambda}_{j}h})\right)^{1/2}R_{k}^{j}\right.$
$\displaystyle\quad-\left.\left[\mathrm{e}^{-\tilde{\lambda}_{j}h}\widetilde{V_{k,j}^{(4)}}+\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\widetilde{f_{N}^{j}(V_{k}^{(4)})}+\left(\frac{q_{j}}{2\tilde{\lambda}_{j}}(1-\mathrm{e}^{-2\tilde{\lambda}_{j}h})\right)^{1/2}R_{k}^{j}\right]\right|$
$\displaystyle=\mathbb{E}\left|\mathrm{e}^{-\tilde{\lambda}_{j}h}(V_{k,j}^{(3)}-\widetilde{V_{k,j}^{(4)}})+\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}[f_{N}^{j}(V_{k}^{(3)})-\widetilde{f_{N}^{j}(V_{k}^{(4)})}]\right|.$

Dealing with these terms separately, we have

$\displaystyle\mathrm{e}^{-\tilde{\lambda}_{j}h}\mathbb{E}|V_{k,j}^{(3)}-V_{k,j}^{(4)}|+\mathrm{e}^{-\tilde{\lambda}_{j}h}\mathbb{E}|V_{k,j}^{(4)}-\widetilde{V_{k,j}^{(4)}}|$
$\displaystyle\leq\ $
$\displaystyle\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{k,j}^{(4)}+\mathrm{e}^{-\tilde{\lambda}_{j}h}\mathbb{E}|\langle
V_{k}^{(4)},e_{j}\rangle-\langle V_{k}^{(4)},e_{j}+\eta_{j}\rangle|$
$\displaystyle=\ $
$\displaystyle\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{k,j}^{(4)}+\mathrm{e}^{-\tilde{\lambda}_{j}h}\mathbb{E}|\langle
V_{k}^{(4)},-\eta_{j}\rangle|$ $\displaystyle\leq\ $
$\displaystyle\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{k,j}^{(4)}+\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{j}^{(\eta)}\mathbb{E}\|V_{k}^{(4)}\|$
(37)

and

$\displaystyle\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\left[\mathbb{E}|f_{N}^{j}(V_{k}^{(3)})-f_{N}^{j}(V_{k}^{(4)})|+\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})-\widetilde{f_{N}^{j}(V_{k}^{(4)})}|\right]$
$\displaystyle\leq\ $
$\displaystyle\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\left[\mathbb{E}|f_{N}^{j}(V_{k}^{(3)})-f_{N}^{j}(V_{k}^{(4)})|+\mathbb{E}|\langle
f_{N}(V_{k}^{(4)}),e_{j}\rangle-\langle
f_{N}(V_{k}^{(4)}),e_{j}+\eta_{j}\rangle|\right]$ $\displaystyle\leq\ $
$\displaystyle\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\left[\mathbb{E}|f_{N}^{j}(V_{k}^{(3)})-f_{N}^{j}(V_{k}^{(4)})|+\mathbb{E}|\langle
f_{N}(V_{k}^{(4)}),-\eta_{j}\rangle|\right]$ $\displaystyle\leq\ $
$\displaystyle\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\left[L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\epsilon_{j}^{(\eta)}\cdot\mathbb{E}\|f_{N}(V_{k}^{(4)})\|\right].$
(38)

Summing (37) and (38), we obtain (36).∎

######  Theorem 4.4.

Suppose that we have fixed parameters $T$, $N$ eigenfunctions, $M$ time steps,
a grid resolution $R$, a Lipschitz constant $L=\|f^{\prime}\|_{\infty}$,
Dirichlet eigenvalues
$0<\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\dots\leq\lambda_{n}$, eigenvalue
error tolerance $\epsilon^{(\lambda)}$, and eigenfunction error tolerance
$\epsilon^{(\eta)}$. We denote the Fourier coefficients $V_{k,j}^{(4)}$ of the
approximated solution $V^{(4)}_{k}$ and the Fourier coefficients
$f_{N}^{j}(V_{k}^{(4)})$ of $f_{N}(U_{t_{k}})$ in the matrices
$\bar{V}=(V_{k,j}^{(4)})$, $f_{N}(\bar{V})=(f_{N}^{j}(V_{k}^{(4)}))$,
$k=1,\dots,M$, $j=1,\dots,N$. Then, for any $\epsilon\in(0,1)$, we have a
bound for the strong error of the complete scheme at time $T$

$\displaystyle\begin{split}\mathbb{E}\|U(T)-V_{M}^{(4)}\|\leq\mathcal{E}^{s}(T,N,M,R,L,\lambda_{1},\lambda_{N},\epsilon^{(\lambda)},\epsilon^{(\eta)},\bar{V},f_{N}(\bar{V}))=\\\
\mathcal{E}_{M}^{(1)}(T,N,M,\lambda_{N})+\mathcal{E}_{M}^{(234)}(T,N,M,R,L,\lambda_{1},\lambda_{N},\epsilon^{(\lambda)},\epsilon^{(\eta)},\bar{V},f_{N}(\bar{V})),\end{split}$
(39)

where, for a constant $C_{T}>0$,

$\displaystyle\mathcal{E}_{M}^{(1)}(T,N,M,\lambda_{N})\leq
C_{T}\left(\lambda_{N}^{\epsilon-1}+\frac{\log(M)}{M}\right)$ (40)

and

$\displaystyle\mathcal{E}_{M}^{(234)}(T,N,M,R,L,\lambda_{1},\epsilon^{(\lambda)},\epsilon^{(\eta)},\bar{V},f_{N}(\bar{V}))$
$\displaystyle\quad\leq
a_{2}^{M}\cdot\epsilon_{0}^{(2)}+b_{2}\cdot\frac{1-a_{2}^{M}}{1-a_{2}}+a_{2}^{M}\cdot\epsilon_{0}^{(3)}+b_{3}\cdot\frac{1-a_{2}^{M}}{1-a_{2}}+\tilde{a}_{2}^{M}\cdot\epsilon_{0}^{(4)}+b_{4}\cdot\frac{1-\tilde{a}_{2}^{M}}{1-\tilde{a}_{2}}+b_{5}$
(41)
$\displaystyle\quad\rightarrow\frac{b_{2}}{1-a_{2}}+\frac{b_{3}}{1-a_{2}}+\frac{b_{4}}{1-\tilde{a}_{2}}+b_{5}\
\text{as}\ M\rightarrow\infty,\ \text{if}\ a_{2}<1\ (\text{since
}\tilde{a}_{2}\leq a_{2}).$

where, for some constant $C>0$,

$\displaystyle a_{2}$ $\displaystyle=\mathrm{e}^{-\lambda_{1}h}+L\cdot
N\cdot\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}},\quad\tilde{a}_{2}=\mathrm{e}^{-\tilde{\lambda}_{1}h}+L\cdot
N\cdot\frac{1-\mathrm{e}^{-\tilde{\lambda}_{1}h}}{\tilde{\lambda}_{1}},$
$\displaystyle b_{2}$
$\displaystyle=\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\cdot N\cdot
C\cdot R^{-2},$ $\displaystyle b_{3}$
$\displaystyle=\mathrm{e}^{-\lambda_{1}h}|1-\mathrm{e}^{-\epsilon^{(\lambda)}h}|\cdot\sup_{k=0,\dots,M}\left[\sum_{j=1}^{N}(\epsilon_{k,j}^{(4)}+\mathbb{E}|V_{k,j}^{(4)}|)\right]$
$\displaystyle\quad+\frac{1}{\lambda_{1}\tilde{\lambda}_{1}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{1}h}\lambda_{1}-\mathrm{e}^{-\lambda_{1}h}\tilde{\lambda}_{1}|\right]$
$\displaystyle\quad\cdot\sup_{k=0,\dots,M}\left[N\cdot
L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\sum_{j=1}^{N}\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right]$
$\displaystyle\quad+N\cdot\left|\sqrt{\frac{1}{2\lambda_{1}}(1-\mathrm{e}^{-2\lambda_{1}h})}-\sqrt{\frac{1}{2(\lambda_{1}+\epsilon^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{1}+\epsilon^{(\lambda)})h})}\right|$
$\displaystyle b_{4}$
$\displaystyle=N\cdot\sup_{j=1,\dots,N}\Biggl{\\{}(1+\mathbb{E}|V^{(4)}_{M,j}|)$
$\displaystyle\quad\cdot\left[\mathrm{e}^{-\tilde{\lambda}_{1}h}\cdot\epsilon^{(\eta)}\cdot\sup_{k=0,\dots,M}\mathbb{E}\|V_{k}^{(4)}\|+\epsilon^{(\eta)}\cdot\frac{1-\mathrm{e}^{-\tilde{\lambda}_{1}h}}{\tilde{\lambda}_{1}}\cdot\sup_{k=0,\dots,M}\mathbb{E}\|f_{N}(V_{k}^{(4)})\|\right]\Biggl{\\}},$
$\displaystyle b_{5}$
$\displaystyle=\epsilon^{(\eta)}\cdot\sum_{j=1}^{N}\mathbb{E}|V_{M,j}^{(4)}|.$

######  Remark 4.5.

For the case of the simpler error bound with $M\rightarrow\infty$, the
condition $a_{2}<1$ is equivalent to $L\cdot N<\lambda_{1}$. This turns out to
be a rather strict condition (see Section 5). The error bound (39) requires
very strict conditions on the parameters in order to be usable, we address
this in Example 5.1. This strictness stems from some crude estimates in the
proofs to simplify expressions, for example the step of replacing all
$\lambda_{n}$ by $\lambda_{1}$ in the proof of Theorem 4.4.

Proof. The error of the scheme can be split up as

$\displaystyle\mathbb{E}(\|U(t_{k})-V_{k}^{(4)}\|)$
$\displaystyle\leq\mathbb{E}(\|U(t_{k})-V_{k}^{(1)}\|)+\mathbb{E}(\|V_{k}^{(1)}-V_{k}^{(2)}\|)+\mathbb{E}(\|V_{k}^{(2)}-V_{k}^{(3)}\|)+\mathbb{E}(\|V_{k}^{(3)}-V_{k}^{(4)}\|)$
$\displaystyle=:\mathbb{E}(\|U(t_{k})-V_{k}^{(1)}\|)+\epsilon_{k}^{(2)}+\epsilon_{k}^{(3)}+\epsilon_{k}^{(4)}\leq\mathcal{E}_{k}^{(1)}+\mathcal{E}_{k}^{(2)}+\mathcal{E}_{k}^{(3)}+\mathcal{E}_{k}^{(4)},\quad
k=1,\dots,M,$

where the error bound $\mathcal{E}_{M}^{(1)}$ in (40) is from [29]. The other
terms relate to the second error term above as

$\displaystyle\mathcal{E}_{k}^{(234)}=\mathcal{E}_{k}^{(2)}+\mathcal{E}_{k}^{(3)}+\mathcal{E}_{k}^{(4)},\quad
k=1,\dots,M.$

The error bounds $\mathcal{E}_{k}^{(i)}$, $i=2,3,4$, are established in Lemmas
4.1, 4.2 and 4.3 as

$\displaystyle\epsilon_{k}^{(2)}\leq\sum_{j=1}^{N}\epsilon_{k,j}^{(2)}=:\mathcal{E}_{k}^{(2)},\quad\epsilon_{k}^{(3)}\leq\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}=:\mathcal{E}_{k}^{(3)},\quad\epsilon_{k}^{(4)}\leq\sum_{j=1}^{N}(1+|V_{k,j}^{(4)}|)\epsilon_{k,j}^{(4)}=:\mathcal{E}_{k}^{(4)}.$

Furthermore, the inequalities (27), (30) and (36) establish relations between
$\epsilon_{k,j}^{(i)}$ and $\epsilon_{k+1,j}^{(i)}$ for $k=0,\dots,M-1$,
$j=1,\dots,N$ and $i=2,3,4$. First, we consider in each time step the sum of
all errors (summing over index $j$) and their bounds given by

$\displaystyle\sum_{j=1}^{N}\epsilon_{k+1,j}^{(2)}$ $\displaystyle\leq\
\sum_{j=1}^{N}\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(2)}+\sum_{j=1}^{N}\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\left[L\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(2)}+D\cdot
R^{-2}\right],$ $\displaystyle\sum_{j=1}^{N}\epsilon_{k+1,j}^{(3)}$
$\displaystyle\leq\
\sum_{j=1}^{N}\mathrm{e}^{-\lambda_{j}h}\epsilon_{k,j}^{(3)}+\sum_{j=1}^{N}\mathrm{e}^{-\lambda_{j}h}|1-\mathrm{e}^{-\epsilon^{(\lambda)}h}|\cdot(\epsilon_{k,j}^{(4)}+\mathbb{E}|V_{k,j}^{(4)}|)$
$\displaystyle\quad+\sum_{j=1}^{N}L\cdot\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(3)}$
$\displaystyle\quad+\sum_{j=1}^{N}\frac{1}{\lambda_{j}\tilde{\lambda}_{j}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{j}h}\lambda_{j}-\mathrm{e}^{-\lambda_{j}h}\tilde{\lambda}_{j}|\right]\cdot\left[L\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(4)}+\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right]$
$\displaystyle\quad+\sum_{j=1}^{N}\left|\sqrt{\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})}-\sqrt{\frac{q_{j}}{2(\lambda_{j}+\epsilon^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{j}+\epsilon^{(\lambda)})h})}\right|,$
$\displaystyle\sum_{j=1}^{N}\epsilon_{k+1,j}^{(4)}\leq\ $
$\displaystyle\sum_{j=1}^{N}\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon_{k,j}^{(4)}+\sum_{j=1}^{N}\mathrm{e}^{-\tilde{\lambda}_{j}h}\epsilon^{(\eta)}\mathbb{E}\|V_{k}^{(4)}\|$
$\displaystyle\quad+\sum_{j=1}^{N}\frac{1-\mathrm{e}^{-\tilde{\lambda}_{j}h}}{\tilde{\lambda}_{j}}\left[L\cdot\sum_{i=1}^{N}\epsilon_{k,i}^{(4)}+\epsilon^{(\eta)}\cdot\mathbb{E}\|f_{N}(V_{k}^{(4)})\|\right].$

Since the terms in the above expressions involving $\lambda_{j}$ like
$\mathrm{e}^{-\lambda_{j}h}$,
$\frac{1-\mathrm{e}^{-\lambda_{j}h}}{\lambda_{j}}$,$\frac{1}{\lambda_{j}\tilde{\lambda}_{j}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{j}h}\lambda_{j}-\mathrm{e}^{-\lambda_{j}h}\tilde{\lambda}_{j}|\right]$
and

$\textstyle\left|\sqrt{\frac{q_{j}}{2\lambda_{j}}(1-\mathrm{e}^{-2\lambda_{j}h})}-\sqrt{\frac{q_{j}}{2(\lambda_{j}+\epsilon^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{j}+\epsilon^{(\lambda)})h})}\right|$

are monotone decreasing in $\lambda_{j}$, we can bound all the above terms
from above by replacing all $\lambda_{j}$ by $\lambda_{1}$. We obtain

$\displaystyle\sum_{j=1}^{N}\epsilon_{k+1,j}^{(2)}$
$\displaystyle\leq\mathrm{e}^{-\lambda_{1}h}\sum_{j=1}^{N}\epsilon_{k,j}^{(2)}+N\cdot
L\cdot\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\sum_{j=1}^{N}\epsilon_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\cdot
N\cdot D\cdot R^{-2}$ $\displaystyle=\left[\mathrm{e}^{-\lambda_{1}h}+N\cdot
L\cdot\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\right]\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(2)}+\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\cdot
N\cdot D\cdot R^{-2},$ (42)
$\displaystyle\sum_{j=1}^{N}\epsilon_{k+1,j}^{(3)}$
$\displaystyle\leq\mathrm{e}^{-\lambda_{1}h}\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}+\mathrm{e}^{-\lambda_{1}h}|1-\mathrm{e}^{-\epsilon^{(\lambda)}h}|\sum_{j=1}^{N}(\epsilon_{k,j}^{(4)}+\mathbb{E}|V_{k,j}^{(4)}|)$
$\displaystyle\quad+N\cdot
L\cdot\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}$
$\displaystyle\quad+\frac{1}{\lambda_{1}\tilde{\lambda}_{1}}$ (43)
$\displaystyle\quad\cdot\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{1}h}\lambda_{1}-\mathrm{e}^{-\lambda_{1}h}\tilde{\lambda}_{1}|\right]\cdot\left[N\cdot
L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\sum_{j=1}^{N}\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right]$
$\displaystyle\quad+N\cdot\left|\sqrt{\frac{1}{2\lambda_{1}}(1-\mathrm{e}^{-2\lambda_{1}h})}-\sqrt{\frac{1}{2(\lambda_{1}+\epsilon^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{1}+\epsilon^{(\lambda)}h)})}\right|$
$\displaystyle\leq\left[\mathrm{e}^{-\lambda_{1}h}+N\cdot
L\cdot\frac{1-\mathrm{e}^{-\lambda_{1}h}}{\lambda_{1}}\right]\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(3)}$
(44)
$\displaystyle\quad+\mathrm{e}^{-\lambda_{1}h}|1-\mathrm{e}^{-\epsilon^{(\lambda)}h}|\cdot\sup_{k=0,\dots,M}\left[\sum_{j=1}^{N}(\epsilon_{k,j}^{(4)}+\mathbb{E}|V_{k,j}^{(4)}|)\right]$
$\displaystyle\quad+\frac{1}{\lambda_{1}\tilde{\lambda}_{1}}\left[\epsilon^{(\lambda)}+|\mathrm{e}^{-\tilde{\lambda}_{1}h}\lambda_{1}-\mathrm{e}^{-\lambda_{1}h}\tilde{\lambda}_{1}|\right]$
$\displaystyle\quad\cdot\sup_{k=0,\dots,M}\left[N\cdot
L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\sum_{j=1}^{N}\mathbb{E}|f_{N}^{j}(V_{k}^{(4)})|\right]$
$\displaystyle\quad+N\cdot\left|\sqrt{\frac{1}{2\lambda_{1}}(1-\mathrm{e}^{-2\lambda_{1}h})}-\sqrt{\frac{1}{2(\lambda_{1}+\epsilon^{(\lambda)})}(1-\mathrm{e}^{-2(\lambda_{1}+\epsilon^{(\lambda)}h)})}\right|,$
$\displaystyle\sum_{j=1}^{N}\epsilon_{k+1,j}^{(4)}$
$\displaystyle\leq\mathrm{e}^{-\tilde{\lambda}_{1}h}\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\mathrm{e}^{-\tilde{\lambda}_{1}h}\cdot
N\cdot\epsilon^{(\eta)}\cdot\mathbb{E}\|V_{k}^{(4)}\|$
$\displaystyle\quad+N\cdot\frac{1-\mathrm{e}^{-\tilde{\lambda}_{1}h}}{\tilde{\lambda}_{1}}\left[L\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}+\epsilon^{(\eta)}\cdot\mathbb{E}\|f_{N}(V_{k}^{(4)})\|\right]$
$\displaystyle\leq\left[\mathrm{e}^{-\tilde{\lambda}_{1}h}+N\cdot
L\cdot\frac{1-\mathrm{e}^{-\tilde{\lambda}_{1}h}}{\tilde{\lambda}_{1}}\right]\cdot\sum_{j=1}^{N}\epsilon_{k,j}^{(4)}$
(45) $\displaystyle\quad+\mathrm{e}^{-\tilde{\lambda}_{1}h}\cdot
N\cdot\epsilon^{(\eta)}\cdot\sup_{k=0,\dots,M}\mathbb{E}\|V_{k}^{(4)}\|$ (46)
$\displaystyle\quad+N\cdot\epsilon^{(\eta)}\cdot\frac{1-\mathrm{e}^{-\tilde{\lambda}_{1}h}}{\tilde{\lambda}_{1}}\cdot\sup_{k=0,\dots,M}\mathbb{E}\|f_{N}(V_{k}^{(4)})\|.$

We can now see that the error sum at $t_{k+1}$ can be bounded by an expression
which is the image of the affine map

$\displaystyle\psi_{a,b}:\mathbb{R}\rightarrow\mathbb{R},\quad\psi_{a,b}(x)=ax+b$

for $a>0,\ b\in\mathbb{R}$. For a given start value, the behavior of the
iteration $\psi_{a,b}^{M}(x_{0})$, $M>0$, depends on whether $a<1$, $a=1$ or
$a>1$. We have the following cases:

$\displaystyle\psi_{a,b}^{M}(x_{0})=\ $ $\displaystyle
a^{M}x_{0}+b\sum_{i=1}^{M-1}a^{i}=a^{M}x_{0}+b\cdot\frac{1-a^{M}}{1-a}\stackrel{{\scriptstyle
M\rightarrow\infty}}{{\longrightarrow}}\frac{b}{1-a}\quad$ $a<1$, (47)
$\displaystyle\psi_{a,b}^{M}(x_{0})=\ $ $\displaystyle x_{0}+Mb,\quad$ $a=1$,
$\displaystyle\psi_{a,b}^{M}(x_{0})=\ $ $\displaystyle
a^{M}\left(x_{0}+\frac{b}{a-1}\right)-\frac{b}{a-1},\quad$ $a>1$. (48)

For our three error sums, the constant $a$ is given by $a_{2}$ and
$\tilde{a}_{2}$, and the constant $b$ is given by $b_{2}$, $b_{3}$ and
$b_{4}$, respectively. Applying the behavior of the iterated affine map
$\psi_{a,b}$ as shown in (47) to (48) to the maps (42), (44) and (45) with the
constants $a_{2}$, $\tilde{a}_{2}$, $b_{2}$, $b_{3}$ and $b_{4}$, as well as
the bounds (26), (29) and (35) shown in the three lemmas, we obtain the terms
shown in (41). ∎

##  5 Numerical experiments.

For the case where the function $f$ associated to the operator $F$ in equation
1 is linear and given by $f(x)=x$, we can compute an exact solution by taking
the function inside the linear operator $A$, so that we consider the
eigenvalues of the operator $A+I$ which are given by $-\lambda_{j}+1$,
$j\in\mathbb{N}$. For the linear function $f_{1}(x)=x$, we present convergence
plots in Figure 6 and compare the approximated solution to the exact solution.
For the nonlinear function $f_{2}(x)=\frac{1}{1+x^{2}}$, we additionally
compare the approximate solution to a reference solution with $N=400$ and
$M=100$ and show the convergence with respect to $N$ in the second plot. Since
the constant for the noise term scales nonlinearly with the step size (see the
final term in (10)), we cannot compute a reference solution to compare the
other approximations to, as we cannot sensibly merge the finer random
increments into more rough ones. For the convergence plots with $f=f_{1}$, the
error was averaged over 100 independent realizations, while for $f=f_{2}$, 10
independent realizations were computed. We can see in all cases that the
convergence with respect to $M$ and $N$ is of order one. We note at this point
that we were not able to verify that the nonlinear function $f_{2}$ indeed
satisfies the conditions (3) to (5), but convergence to a reference solution
is observed nevertheless. We also remark that if strictly nonlinear $f$ is
used, the computational effort is substantially higher since in every time
step $N$ two-dimensional numerical integrals have to be computed.
We now shed some light on how the error bound from Theorem 4.4 behaves here:
We first note that for the Peanut shape it is $\lambda_{1}\approx 6.5155$, so
in order for the condition $a_{2}<1$ to be fulfilled we would need $L\cdot
N<6.5155.$ This is a very strict condition as either $N$ or $L$ needs to be
quite small; for instance, if $N=100$, only nonlinearities with $L<0.065155$
are covered by the case $M\rightarrow\infty$ in Theorem 4.4. But as this is
not a necessary condition to use Theorem 4.4, we can still give bounds also if
$a_{2}>1$.

######  Example 5.1.

We give a concrete example for an error bound resulting from Theorem 4.4.
Suppose $T=0.1$, $N=100$, $M=50$, $R=301$, $L=0.01$, $\lambda_{1}=6.5155$,
$|V_{k,j}^{(4)}|\leq 1$ for all $k=0,\dots,M$, $j=1,\dots,N$ (we have observed
in practice that this has always been the case). We also assume $C=1$ (as
Simpson’s rule is very accurate for smooth functions) and
$\epsilon_{k,j}^{(4)}\leq 1$ for all $k,j$ (which is a very loose assumption,
but does not matter much as the corresponding terms in $b_{3}$ are small
anyway). We assume error tolerances $\epsilon^{(\lambda)}=2\cdot 10^{-4}$ and
$\epsilon^{(\eta)}=5\cdot 10^{-6}$. We also assume
$\epsilon_{0}^{(2)}=\epsilon_{0}^{(3)}=\epsilon_{0}^{(4)}=10^{-4}$ (as the
initial condition is smooth, its approximation can be assumed to be accurate).
With the above assumptions on the parameters, we obtain $a_{2}\approx 0.98904$
(with $\tilde{a}_{2}$, $b_{2}\approx 2.19316\cdot 10^{-6}$, $b_{3}\approx
1.96257\cdot 10^{-4}$, $b_{4}\approx 1.83156\cdot 10^{-4}$ and $b_{5}=5\cdot
10^{-5}$ and Theorem 4.4 yields a total error bound of
$\mathcal{E}_{M}^{(234)}\leq 0.2234$. While this is clearly not a satisfactory
bound especially compared to $\mathcal{E}_{M}^{(1)}$, we note that this bound
depends critically on the error tolerances $\epsilon^{(\lambda)}$ and
$\epsilon^{(\eta)}$. For example, for $\epsilon^{(\lambda)}=2\cdot 10^{-5}$
and $\epsilon^{(\eta)}=5\cdot 10^{-7}$, leaving all other parameters
unchanged, we obtain $\mathcal{E}_{M}^{(234)}\leq 0.02257$. This shows that
the error bound is usable in practice, but only with an as yet infeasible
effort of computation (see Figure 5).

In the plots in Figure 7, we give a visual impression of a solution of an SPDE
(1) on the Peanut shape by showing realizations of the approximated solution
again for noise with $q_{n}=n^{-2}$, $A=\Delta$, $N=400$, $M=100$ on the
asymmetric Peanut shape introduced in Section 3. For each simulation, the
initial condition was picked to be the bump function

$\displaystyle
B_{x_{1},x_{2},y_{1},y_{2}}(x,y)=\begin{cases}\exp\left(-\frac{1}{1-r^{2}}\right),\quad&r^{2}=\left(x-\frac{x_{1}+x_{2}}{2}\right)^{2}+\left(y-\frac{y_{1}+y_{2}}{2}\right)^{2}<1,\\\
0,\ &r^{2}\geq 1\end{cases}$

supported on an ellipse within the rectangle
$[x_{1},x_{2}]\times[y_{1},y_{2}]$, where we picked $[x_{1},x_{2}]=[0.4,0.6]$,
$[y_{1},y_{2}]=[0.3,0.5]$. The random input, i.e. the realizations of the
Wiener increments was also the same for each simulation. Only the nonlinear
function was changed for the different simulations: It was chosen to be a
function

$\displaystyle b_{p}(x)=\exp(-10\cdot(x-p)^{2}),$ (49)

with one maximum at $p\in\\{-0.4,-0.2,0.2,0.4\\}$.

$10^{1}$$10^{2}$$10^{-4}$$10^{-5}$$M$$\|\tilde{U}_{T}^{M_{\mathrm{ref}},N}-\tilde{U}_{T}^{M,N}\|$
$10^{1}$$10^{2}$$1\cdot 10^{-2}$$1\cdot 10^{-3}$$2\cdot
10^{-4}$$N$$\|\tilde{U}_{T}^{M,N_{\mathrm{ref}}}-\tilde{U}_{T}^{M,N}\|$$f=f_{1}$$f=f_{2}$
---
Figure 6: Convergence tests for the Peanut shape. The first plot shows the
convergence with respect to $M$ for the linear function $f_{1}(x)=x$, the
second plot shows the convergence with respect to $N$ for $f_{1}(x)=x$ as well
as $f_{2}(x)=\frac{1}{1+x^{2}}$. The dashed lines in blue are reference lines
for an order of convergence equal to one.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$x$$y$
$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$x$$y$
---
$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$x$$y$
$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$$x$$y$

Figure 7: The approximated solution of the SPDE (1) on the Peanut shape for
$N=400$, $T=0.1$, $M=100$ and the nonlinear functions (49) for
$p\in\\{-0.4,-0.2,0.2,0.4\\}$ from top-left to bottom-right. It appears that
the process is ‘dimmed down’ from using nonlinear functions with a peak of
higher absolute value.

##  6 Outlook.

The advantages of the presented method are its flexibility in choosing the
two-dimensional domain and the possibility to freely choose the points inside
the domain at which the eigenfunctions are computed. This renders the solution
of SPDEs with high resolution on a big class of domains possible. Due to the
general setting of the spectral method, the scheme could be extended to a
larger class of differential operators such as general second-order elliptic
differential operators, whose Dirichlet eigenfunctions also form a complete
orthogonal basis in $L^{2}(D)$ [10]. Difficulties might arise while using a
boundary integral equation method, as explicit knowledge of the fundamental
solution of the differential operator is required. The scheme could also be
applied to three-dimensional domains with a $\mathcal{C}^{2}$ boundary: In
order to do this, the fundamental solution needs to be changed to the one for
the three-dimensional Helmholtz operator, the curved elements are now two-
dimensional (see [34, Section 5]) and the linear system resulting from
applying Beyn’s contour integral algorithm is larger. However, for a small
error, a very large number of eigenfunctions would have to be computed due to
Weyl’s law (7) and (40).
The requirements on the nonlinear operator $F$ could be relaxed by
implementing a different scheme like the one from [27], where the requirements
on $F$ are less strict. Future work might also include deriving error bounds
for higher moments such as $\mathbb{E}[\|U(t_{k})-V_{k}^{(4)}\|^{2}]^{1/2}$.

## Acknowledgments

We thank Arnulf Jentzen for valuable comments and a fruitful conversation.

##  Appendix A Joint eigenbasis of $Q$ and $\Delta$ in two dimensions.

Let $D\subset\mathbb{R}^{d}$ be a bounded domain, $d\geq 1$, with simple
Dirichlet eigenvalues. It is well-known that a $Q$-Wiener process always has
the series expansion given by (6) in terms of the eigenfunctions and
eigenvalues of $Q$ [16, Proposition 4.3]. In our setting (see end of Section
3), we have
$\sum_{n=1}^{\infty}\|Qe^{Q}_{n}\|^{2}=\sum_{n=1}^{\infty}q_{n}^{2}=\sum_{n=1}^{\infty}n^{-(4+2\epsilon_{q})}<\infty$,
so $Q$ is a Hilbert-Schmidt operator. Furthermore, as is recalled in [7,
Section 3], the operator $Q$ is Hilbert-Schmidt if and only if there is a
kernel $q\in L^{2}(D\times D)$ such that

$\displaystyle(Qf)(x)=\int_{D}q(x,y)f(y)\;\mathrm{d}y,\quad f\in L^{2}(D).$
(50)

It is shown in [7, Theorem 3.3] that if $Q$ is trace class (in our setting, it
is, since $\mathrm{tr}(Q)=\sum_{n=1}^{\infty}q_{n}<\infty$), the kernel $q$
can be understood as the correlation functional of a generalized Gaussian
stochastic process (GSP) $\xi=\partial_{t}W$ satisfying

$\displaystyle\mathbb{E}[\xi(t,x)]=0,\quad\mathbb{E}[\xi(t,x)\xi(s,y)]=\delta(t-s)q(x,y),$
(51)

with $0\leq s\leq t\leq T$ and $x,y\in D$. (We note here that this
correspondence to a GSP relies on a geometric assumption [7, Assumption 3.1]
on $D$ which is not fulfilled by the Peanut shape considered in this paper,
but it is also remarked on [7, p. 259] that the results are true for ‘much
more general domains’.)
A necessary and sufficient condition for $Q$ to have a joint eigenbasis with
$\Delta$ is given in terms of $q$ [7, Theorem 4.2]: If
$q\in\mathcal{C}^{2}(D\times D)$, then $Q$ and $\Delta$ have a joint
eigenbasis if and only if

$\displaystyle\Delta_{y}q(x,y)=\Delta_{x}q(x,y)\quad\text{and }\quad
q(\cdot,\zeta),\ q(\zeta,\cdot)\in\mathcal{C}_{\tau}^{2}(D)$ (52)

for all $x,y,\zeta\in D$, where
$\mathcal{C}^{2}_{\tau}(D)=\\{u\in\mathcal{C}^{2}(D)\ |\ u=0\ \text{on
}\partial D\\}$, and with $\Delta_{x}$ denoting the Laplacian in $x$ with
Dirichlet boundary conditions and likewise for $y$.
Solving the first equation in (52) by separation of variables, we write
$q(x,y)=q_{1}(x)q_{2}(y)$ and obtain

$\displaystyle
q_{1}(x)\Delta_{y}q_{2}(y)=\Delta_{x}q_{1}(x)q_{2}(y)=\lambda\Leftrightarrow\frac{\Delta_{y}q_{2}(y)}{q_{2}(y)}=\frac{\Delta_{x}q_{1}(x)}{q_{1}(x)}=\lambda$
(53)

for a constant $\lambda$, for all $x,y\in D$ where $q_{2}(y)\neq 0$,
$q_{1}(x)\neq 0$. We retrieve from (53) the Helmholtz equations

$\displaystyle\Delta q_{1}(x)=\lambda q_{1}(x),\quad\Delta q_{2}(y)=\lambda
q_{2}(y)$ (54)

with homogeneous Dirichlet boundary conditions for $q_{1}$ and $q_{2}$, with
the solutions $\\{e_{n}\\}_{n\in\mathbb{N}}$,
$\\{\lambda_{n}\\}_{n\in\mathbb{N}}$ the Dirichlet eigenfunctions and
eigenvalues on $D$, resulting in

$\displaystyle
q(x,y)=\left[\sum_{i=1}^{\infty}a_{i}e_{i}(x)\right]\cdot\left[\sum_{j=1}^{\infty}b_{j}e_{j}(y)\right]=\sum_{i=1}^{\infty}\sum_{j=1}^{i}a_{i}b_{i+1-j}e_{i}(x)e_{i+1-j}(y)$
(55)

for real coefficients $a_{n},b_{n}$, $n\in\mathbb{N}$. The first condition in
(52) now implies

$\displaystyle\sum_{i=1}^{\infty}\sum_{j=1}^{i}\lambda_{i}a_{i}b_{i+1-j}e_{i}(x)e_{i+1-j}(y)=\sum_{i=1}^{\infty}\sum_{j=1}^{i}\lambda_{i+1-j}a_{i}b_{i+1-j}e_{i}(x)e_{i+1-j}(y)$
(56)

and since the eigenvalues are simple, due to orthogonality only the terms
remain where $i=i+1-j$. For real coefficients $c_{n}$, $n\in\mathbb{N}$, we
obtain

$\displaystyle q(x,y)=\sum_{n=1}^{\infty}c_{n}e_{n}(x)e_{n}(y),$ (57)

and in our setting in two dimensions we choose $c_{n}=n^{-(2+\epsilon_{q})}$
(see Section 2). Indeed, as $Q$ is symmetric and positive semidefinite [16, p.
47 f.], and therefore so is $q$ [7, p. 260], due to Mercer’s theorem [15, p.
138], $q$ always has a series representation (57) in terms of the
eigenfunctions of $Q$, and we have shown that the Dirichlet eigenfunctions are
the suitable choice for the Dirichlet boundary conditions.
On the second part in (52), we note that

$\displaystyle\Delta_{x}q(x,y)=\sum_{n=1}^{\infty}n^{-(2+q_{\epsilon})}\lambda_{n}e_{n}(x)e_{n}(y),$
(58)

is absolutely convergent due to
$n^{-(2+q_{\epsilon})}\lambda_{n}=\mathcal{O}(n^{-(1+q_{\epsilon})})$ being
summable. We also note that the family $\\{e_{n}\\}_{n\in\mathbb{N}}$ is
uniformly bounded in the $\|\cdot\|_{\infty}$ norm since due to the Hölder
inequality

$\displaystyle\|1\cdot
e_{n}^{2}\|_{1}\leq\|1\|_{1}\|e_{n}^{2}\|_{\infty}\Leftrightarrow\|e_{n}\|_{2}^{2}\leq\mu(D)\cdot\|e_{n}\|_{\infty}^{2},$
(59)

so $\|e_{n}\|_{\infty}\leq\mu(D)^{-1/2}$, as $\|e_{n}\|_{2}=1$ for each
$n\in\mathbb{N}$, where $\mu(D)$ is the area of $D$. It now follows from the
Weierstraß convergence criterion [20, p. 259] that (58) converges absolutely
and uniformly, and therefore as a uniform limit of continuous functions is
continuous.

## References

  * [1] A. Alabert and I. Gyöngy.  On Numerical Approximation of Stochastic Burgers’ Equation, pages 1–15.  In: Y. Kabanov, R. Lipster and J. Stoyanov (editors) From Stochastic Calculus to Mathematical Finance. Springer, 2006.
  * [2] E. J. Allen, S. J. Novosel, and Z. Zhang.  Finite element and difference approximation of some linear stochastic partial differential equations.  Stochastics and Stochastic Reports, 64(1–2):117–142, 1998.
  * [3] K. E. Atkinson.  The Numerical Solution of Integral Equations of the Second Kind. Cambridge University Press, 1997.
  * [4] D. J. Benson.  Music: A Mathematical Offering.  Cambridge University Press, 2006.
  * [5] W.-J. Beyn.  An integral method for solving nonlinear eigenvalue problems.  Linear Algebra and its Applications, 436(10):3839–3863, 2012.
  * [6] D. Blömker and A. Jentzen.  Galerkin approximations for the stochastic Burgers equation.  SIAM Journal on Numerical Analysis, 51(1):694–715, 2009.
  * [7] D. Blömker.  Nonhomogeneous noise and Q-Wiener processes on bounded domains.  Stochastic Analysis and Applications, 23(2):255–273, 2005.
  * [8] S. Bonaccorsi and G. Guatteri.  Classical solutions for SPDEs with Dirichlet boundary conditions  In: Seminar on Stochastic Analysis, Random Fields and Applications III: Centro Stefano Franscini, Ascona, September 1999, pages 33–44  Springer, 2002.
  * [9] C.-E. Bréhier and G. Vilmart.  High order integrator for sampling the invariant distribution of a class of parabolic stochastic PDEs with additive space-time noise.  SIAM Journal on Scientific Computing, 38(4):A2283–A2306, 2016.
  * [10] F. E. Browder.  On the eigenfunctions and eigenvalues of the general linear elliptic differential operator.  Proceedings of the National Academy of Sciences, 39(5):433–439, 1953.
  * [11] O. P. Bruno and F. Reitich.  High-order methods for high-frequency scattering applications.  In H. Ammari, editor, Modeling and Computations in Electromagnetics, pages 129–163. Springer, 2008.
  * [12] S. N. Chandler-Wilde, D. P. Hewett, S. Langdon, and A. Twigger.  A high frequency boundary element method for scattering by a class of nonconvex obstacles.  Numerische Mathematik, 129(4):647–689, 2015.
  * [13] P.-L. Chow.  Stochastic Partial Differential Equations.  CRC Press, 2007.
  * [14] D. Colton and R. Kress.  Inverse Acoustic and Electromagnetic Scattering Theory.  Springer, 2013.
  * [15] R. Courant and D. Hilbert.  Methods of Mathematical Physics, Vol. I.  Interscience Publishers, Inc., 1953.
  * [16] G. Da Prato and J. Zabczyk.  Stochastic Equations in Infinite Dimensions.  Cambridge University Press, 2nd edition, 2014.
  * [17] A. M. Davie and J. G. Gaines.  Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations.  Mathematics of Computation, 70(233):121–134, 2000.
  * [18] G. Di Nunno, S. Ortiz-Latorre and A. Petersson.  SPDE bridges with observation noise and their spatial approximation  Stochastic Processes and their Applications, 158:170–207, 2023.
  * [19] M. Dozzi and J. A. López-Mimbela.  Finite-time blowup and existence of global positive solutions of a semi-linear SPDE  Stochastic Processes and their Applications, 120(6):767–776, 2010.
  * [20] O. Forster.  Analysis 1. Springer, 11th edition, 2013.
  * [21] G. Evans.  Practical Numerical Integration. Wiley, 1993
  * [22] M. Gradinaru, I. Nourdin, and S. Tindel.  Itô’s and Tanaka’s-type formulae for the stochastic heat equation: the linear case.  Journal of Functional Analysis, 228(1):114–143, 2005.
  * [23] D. S. Grebenkov and B.-T. Nguyen.  Geometrical structure of Laplacian eigenfunctions.  SIAM Review, 55(4):601–667, 2013.
  * [24] W. Grecksch and P. E. Kloeden.  Time-discretised Galerkin approximations of parabolic stochastic PDEs.  Bulletin of the Australian Mathematical Society, 54(1):79–85, 1996\.
  * [25] E. Hausenblas.  Approximation for semilinear stochastic evolution equations.  Potential Analysis, 18(2):141–186, 2003.
  * [26] A. Jentzen.  Pathwise numerical approximation of SPDEs with additive noise under non-global Lipschitz coefficients.  Potential Analysis, 31:375–404, 2009.
  * [27] A. Jentzen.  Higher order pathwise numerical approximations of SPDEs with additive noise.  SIAM Journal on Numerical Analysis, 49(2):642–667, 2011.
  * [28] A. Jentzen, P. E. Kloeden, and G. Winkel.  Efficient simulation of nonlinear parabolic SPDEs with additive noise.  The Annals of Applied Probability, 21(3):908–950, 2011.
  * [29] A. Jentzen and P. E. Kloeden.  Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise.  Proceedings of the Royal Society A, 465:649–667, 2009.
  * [30] M. A. Katsoulakis, G. T. Kossioris, and O. Lakkis.  Noise regularization and computations for the 1-dimensional stochastic Allen-Cahn problem.  Interfaces and Free Boundaries, 9(1):1–30, 2007.
  * [31] A. Kleefeld.  A numerical method to compute interior transmission eigenvalues.  Inverse Problems, 29(10):104012, 2013.
  * [32] A. Kleefeld.  Shape optimization for interior Neumann and transmission eigenvalues.  In C. Constanda and P. Harris, editors, Integral Methods in Science and Engineering. Birkhäuser, 2019.
  * [33] A. Kleefeld.  The hot spots conjecture can be false: some numerical examples.  Advances in Computational Mathematics, 47(6):85, 2021.
  * [34] A. Kleefeld and T.-C. Lin.  The nonlinear Landweber method applied to an inverse scattering problem for sound-soft obstacles in 3D  Computer Physics Communications, 182(12):2550–2560, 2011.
  * [35] W. Liu and M. Röckner.  Stochastic Partial Differential Equations: An Introduction.  Springer, 2015.
  * [36] G. J. Lord and A. Tambue.  Stochastic exponential integrators for the finite element discretization of SPDEs for multiplicative and additive noise.  IMA Journal of Numerical Analysis, 33(2):515–543, 2013.
  * [37] S. V. Lototsky and B. L. Rozovsky.  Stochastic Partial Differential Equations.  Springer, 2017.
  * [38] P. Pušnik.  Strong convergence rates for full-discrete numerical approximations of stochastic partial differential equations with non-globally Lipschitz continuous nonlinearities.  PhD thesis, ETH Zurich, 2020.
  * [39] T. Shardlow.  Numerical methods for stochastic parabolic PDEs.  Numerical Functional Analysis and Optimization, 20(1–2):121–145, 1999.
  * [40] J. Sun and A. Zhou.  Finite Element Methods for Eigenvalue Problems.  CRC Press, 2016.
  * [41] J. B. Walsh.  An introduction to stochastic partial differential equations.  In: Hennequin, P. L. (eds) École d’Été de Probabilités de Saint Flour XIV - 1984. Lecture Notes in Mathematics, vol 1180. Springer, 2006.
  * [42] J. B. Walsh.  Finite element methods for parabolic stochastic PDEs.  Potential Analysis, 23:1–43, 2005.