A splitting/polynomial chaos expansion approach for stochastic evolution equations

TL;DR

This paper introduces a novel approach combining splitting methods with polynomial chaos expansion to efficiently solve stochastic evolution equations, supported by convergence analysis and numerical validation.

Contribution

It develops a new method integrating deterministic splitting with polynomial chaos for stochastic PDEs, enhancing solution efficiency and accuracy.

Findings

Method effectively reduces stochastic PDEs to deterministic systems.

Convergence analysis confirms the method's reliability.

Numerical experiments validate the approach's accuracy.

Abstract

In this paper, we combine deterministic splitting methods with a polynomial chaos expansion method for solving stochastic parabolic evolution problems. The stochastic differential equation is reduced to a system of deterministic equations that we solve explicitly by splitting methods. The method can be applied to a wide class of problems where the related stochastic processes are given uniquely in terms of stochastic polynomials. A comprehensive convergence analysis is provided and numerical experiments validate our approach.

Figures

Tables

Table 1. Table 4.1. Average computational time (in seconds) for the calculation of one solution u m subscript 𝑢 𝑚 u_{m} for different degrees of freedom N 𝑁 N , i.e., the number of computational points used in the discretization of D D {\rm D} and the operator ℒ ℒ \mathcal{L} . The employed methods are: Crank–Nicolson (CN), modified Lie splitting (MLSPL) and trapezoidal splitting (TSPL).

$N\times N$	CN $\left[s\right]$	MLSPL $\left[s\right]$	TSPL $\left[s\right]$
$4\times 4$	0.0133	0.0373	0.0530
$8\times 8$	0.0214	0.0241	0.0379
$16\times 16$	0.1005	0.1008	0.0948
$32\times 32$	0.4098	0.3652	0.4047
$64\times 64$	2.7620	1.8051	1.8237
$128\times 128$	41.1284	9.3921	13.5091

Equations

$$\begin{split}du(t)&=\big{(}(A+B)\,u(t)+f(t)\big{)}\,dt+\big{(}C\,u(t)+g(t)\big{)}\,dB(t)\\ u(0)&=u^{0},\end{split}$$

$$\begin{array}{rl}{u}_t(t,x,\omega )& =(A+B)u(t,x,\omega )+G(t,x,\omega )\\ u(0,x,\omega )& =u^0(x,\omega ),\end{array}$$

$$(F(\omega ),G(\omega ){)}_{(L{)}^2}=\mathbb{E}(F(\omega )G(\omega )),$$

$$\begin{array}{rl}{h}_0(x)& =1,\\ h_1(x)& =x,\\ h_{n+1}(x)& =x{h}_n(x)+n{h}_{n-1}(x)\text{for}n\ge 2,x\in \mathbb{R}.\end{array}$$

$$\begin{array}{r}{H}_{\alpha }(\mathbf{\xi }(\omega ))=H_{({\alpha }_1,{\alpha }_2,...)}(({\xi }_1(\omega ),{\xi }_2(\omega ),...))=\prod _{i\in \mathbb{N}}{h}_{{\alpha }_i}({\xi }_i(\omega )),\end{array}$$

$$\begin{array}{rl}{H}_{(0,0,\dots )}(\mathbf{\xi }(\omega ))& =\prod _{i\in \mathbb{N}}{h}_0({\xi }_i(\omega ))=1\text{and}\\ H_{{\epsilon }^{(k)}}(\mathbf{\xi }(\omega ))& =h_1({\xi }_k(\omega ))\prod _{i\ne k,i\in \mathbb{N}}{h}_0({\xi }_i(\omega ))={\xi }_k(\omega ),k\in \mathbb{N}.\end{array}$$

$$F(\mathbf{\xi }(\omega ))=\sum _{\alpha \in \mathcal{I}}{f}_{\alpha }{H}_{\alpha }(\mathbf{\xi }(\omega )),$$

$$\Vert F{\Vert }_{(L{)}^2}^2=\sum _{\alpha \in \mathcal{I}}{f}_{\alpha }^2\Vert H_{\alpha }{\Vert }_{(L{)}^2}^2<\infty .$$

$$\mathbb{E}(H_{\alpha }{H}_{\beta })=\alpha !{\delta }_{\alpha \beta },$$

$$\mathbb{E}(H_{(0,0,\dots )})=1\text{and}\mathbb{E}(H_{\alpha })=0\text{for}|\alpha |>0.$$

$${\int }_{\mathbb{R}}{h}_n(x)h_m(x)d\mu (x)=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}{h}_n(x)h_m(x)e^{-\frac{x^2}{2}}dx=n!{\delta }_{m,n}$$

$$\begin{array}{rl}{p}_0(x)& =1,\\ p_1(x)& =x,\\ (n+1)p_{n+1}(x)& =(2n+1)x{p}_n(x)-n{p}_{n-1}(x)\text{for}n\ge 1.\end{array}$$

$$p_n(x)=\frac{1}{2^{n}n!}\frac{d^n}{d{x}^n}(x^2-1{)}^n.$$

$${\int }_{-1}^1{p}_m(x)p_n(x)dx=\frac{2}{2n+1}{\delta }_{m,n},m,n\in {\mathbb{N}}_0.$$

$$L_{\alpha }(\mathbf{\xi }(\omega ))=\prod _{i\in \mathbb{N}}{p}_{{\alpha }_i}({\xi }_i(\omega )),\alpha =({\alpha }_1,{\alpha }_2,\dots )\in \mathcal{I},$$

$$\begin{array}{rl}{L}_{(0,0,\dots )}(\mathbf{\xi }(\omega ))& =1\text{and}\\ L_{{\epsilon }^{(k)}}(\mathbf{\xi }(\omega ))& ={\xi }_k(\omega )\text{for}k\in \mathbb{N}.\end{array}$$

$$\mathbb{E}(L_{(0,0,\dots )})=1\text{and}\mathbb{E}(L_{\alpha }(\mathbf{\xi }(\omega )))=0\text{for}|\alpha |>0,$$

$$\mathbb{E}(L_{\alpha }{L}_{\beta })=\mathbb{E}{L}_{\alpha }^2{\delta }_{\alpha ,\beta }=\frac{1}{{\prod }_{k\in \mathbb{N}}(2{\alpha }_k+1)}{\delta }_{\alpha \beta }$$

$$F(\mathbf{\xi }(\omega ))=\sum _{\alpha \in \mathcal{I}}{f}_{\alpha }{L}_{\alpha }(\mathbf{\xi }(\omega )),$$

$$f_{\alpha }=\frac{1}{\mathbb{E}(L_{\alpha }^2)}\mathbb{E}(F{L}_{\alpha }),\alpha \in \mathcal{I}$$

$$\Vert F{\Vert }_{(L{)}^2}^2=\sum _{\alpha \in \mathcal{I}}{f}_{\alpha }^2\mathbb{E}{L}_{\alpha }^2=\sum _{\alpha \in \mathcal{I}}\frac{f_{\alpha }^2}{{\prod }_{k\in \mathbb{N}}(2{\alpha }_k+1)}<\infty .$$

$$u(t,x,\omega )=\sum _{\alpha \in \mathcal{I}}{u}_{\alpha }(t,x)L_{\alpha }(\mathbf{\xi }(\omega )).$$

$$\sum _{\alpha \in \mathcal{I}}\Vert u_{\alpha }{\Vert }_{C([0,T],X)}^2\mathbb{E}{L}_{\alpha }^2=\sum _{\alpha \in \mathcal{I}}\underset{t\in [0,T]}{\sup }\Vert u_{\alpha }(t){\Vert }_X^2\mathbb{E}{L}_{\alpha }^2<\infty$$

$$\sum _{\alpha \in \mathcal{I}}\Vert u_{\alpha }{\Vert }_{C^1([0,T],X)}^2\mathbb{E}{L}_{\alpha }^2<\infty .$$

$$v(x,\omega )=\bar{v}(x)+\sum _{k\in \mathbb{N}}\sqrt{{\lambda }_k}{e}_k(x)Z_k(\omega ),$$

$${\int }_{\mathrm{D}}{C}_v(x_1,x_2)e_k(x_2)d{x}_2={\lambda }_k{e}_k(x_1),x_1\in \mathrm{D},k\in \mathbb{N},$$

$$\tilde{v}(x,\omega )=\bar{v}(x)+\sum _{k=1}^n\sqrt{{\lambda }_k}{e}_k(x)Z_k(\omega ),$$

$$v(x,\omega )=\bar{v}(x)+\sum _{k\in \mathbb{N}}{v}_{{\epsilon }^{(k)}}(x)L_{{\epsilon }^{(k)}}(\mathbf{Z}(\omega )),$$

$$\tilde{v}(x,\omega )=\bar{v}(x)+\sum _{k=1}^n{v}_{{\epsilon }^{(k)}}(x)L_{{\epsilon }^{(k)}}(\mathbf{Z}(\omega )).$$

$$\mathcal{L}u=\sum _{\alpha \in \mathcal{I}}\mathcal{L}(u_{\alpha })L_{\alpha }$$

Full text

1Andreas KOFLER, Department of Radiology, Charité-Universitätsmedizin Berlin, Germany, andreas.kofler@charite.de

2Tijana LEVAJKOVIĆ, Institute of Stochastics and Business Mathematics, Vienna University of Technology, Austria, tijana.levajkovic@tuwien.ac.at

3Hermann MENA, Department of Mathematics, Yachay Tech University, Urcuquí, Ecuador, mena@yachaytech.edu.ec & Department of Mathematics, University of Innsbruck, Austria, hermann.mena@uibk.ac.at

4Alexander OSTERMANN, Department of Mathematics, University of Innsbruck, Austria, alexander.ostermann@uibk.ac.at

1. Introduction

Splitting methods are numerical methods for solving differential equations, both ordinary and partial differential equations (PDEs), involving operators that are decomposable into a sum of (differential) operators. These methods are used to improve the speed of calculations for problems involving decomposable operators and to solve multidimensional PDEs by reducing them to a sum of one-dimensional problems [10]. Splitting methods have been successfully applied to many types of PDEs, e.g. [14], [16]. Exponential splitting methods are applied in cases when the explicit solution of a splitted equation can be computed. Such computations often rely on applying fast Fourier techniques, see for instant [38]. Resolvent splitting is used in cases when the splitted equation cannot be solved explicitly [17], [34]; here we consider this type of methods.

There are also many results in the literature about the approximation of solutions of SPDEs using splitting methods, see e.g. [2], [3], [4], [5], [9], [12], [15] and references therein. In [12] a splitting method for nonlinear stochastic equations of Schrödinger type is proposed. There the authors approximate the solution of the problem by a sequence of solutions of two types of equations: one without stochastic term and other containing only the stochastic term. They prove that an appropriate combination of the solutions of these equations converges strongly to the solution of the original problem. Exponential integrators for nonlinear Schrödinger equations with white noise dispersion were proposed in [5]. For a stochastic incompressible time-dependent Stokes equation different time-splitting methods were studied in [4]. In [2] the convergence of a Douglas–Rachford type splitting algorithm is presented for general SPDEs driven by linear multiplicative noise. In this work a splitting/polynomial chaos expansion is considered for stochastic evolution equations. Our approach has not been considered in the literature for solving these types of SPDEs so far.

We consider stochastic evolution equations of the form

[TABLE]

where $A$, $B$ and $C$ are differential operators acting on Hilbert space valued stochastic processes, $\{B_{t}\}_{t\geq 0}$ is a cylindrical Brownian motion on a given probability space $(\Omega,\mathcal{F},\mathbb{P})$ and $f$ and $g$ are deterministic functions. In [30] equation (1.1) involving Gaussian noise terms was solved in an appropriate weighted Wiener chaos space. The deterministic problem that corresponds to (1.1), i.e., the case where $C=0$ and $g=0$, for particular $Au=\partial_{x}(a\partial_{x}u)$, $Bu=\partial_{y}(b\partial_{y}u)$ and $f$ was studied in [10]. We consider equation (1.1) involving a non-Gaussian noise term. Namely, we consider inhomogeneous parabolic evolution equations involving the operators that can be split in $A+B$ and uniformly distributed random inputs. These equations, can be also written in the form

[TABLE]

where $G$ represents the noise term, see e.g. [20], [25], [28], [29], [30]. The existence of a random parameter $\omega$ is due to uncertainties coming from initial conditions and/or a random force term. Therefore, the solution is considered to be a stochastic process.

Stochastic processes with finite second moments on white noise spaces can be represented in series expansion form in terms of a family of orthogonal stochastic polynomials. The classes of orthogonal polynomials are chosen depending on the underlying probability measure [19], [20]. Namely, the Askey scheme of hypergeometric orthogonal polynomials and the Sheffer system [36], [37] can be used to define several discrete and continuous distribution types [39]. For example, in the case of the Gaussian measure, the orthogonal basis of the space of random variables with finite second moments is constructed by the use of the Hermite polynomials. We consider problems with non-Gaussian random inputs. The noise term is considered to be uniformly distributed. It is known that in order to obtain a square integrable solution of (1.1) with deterministic initial condition, it is enough to assume that the operator $A-\frac{1}{2}CC^{*}$ is elliptic and that the stochastic part (the noise term) is sufficiently regular, see e.g. [8]. In this work, the assumptions on the input data for problem (1.2) will be set such that the existence of a square integrable solution is always established. We do not consider solutions which are generalized stochastic processes as in [28], [30], since our focus is on numerical treatment.

Our approach is general enough to be applied to problems with additive noise, problems involving multiplicative noise and problems with convolution-type noise [28]. For instance, with this approach the heat equation with random potential, the heat equation in random (inhomogeneous and anisotropic) media and the Langevin equation can be solved. If (1.1) does have a sufficiently regular solution, this solution can be projected on an orthonormal basis in some Hilbert space, resulting in a system of equations for the corresponding Fourier coefficients. Thus, we use the so-called polynomial chaos method or the chaos expansion method and define the solution of (1.1) as a formal Fourier series with the coefficients computed by solving the corresponding system of deterministic PDEs [30]. With this method, the deterministic part of a solution is separated from its random part. Particularly, in the case of Gaussian noise, the orthonormal basis of stochastic polynomials involves the Hermite polynomials and in the case when the noise term is uniformly distributed, the orthonormal basis involves the Legendre polynomials [36]. By construction, the solution is strong in the probabilistic sense. It is uniquely determined by the coefficients, free terms, initial condition and the noise term. The coefficients in the Fourier series are uniquely determined by equation (1.1) and are computed by solving (numerically) the corresponding lower-triangular system of deterministic parabolic equations. The polynomial chaos method has been successfully applied for solving general classes of SPDEs. The list of references is long, here we mention just a few [20], [28], [32], [33]. In [25], [26], [27] this approach has been recently applied to the stochastic optimal regulator control problem [13].

Practical application of the Wiener polynomial chaos involves two truncations, truncation with respect to the number of the random variables and truncation with respect to the order of the orthogonal Askey polynomials used (in the particular case considered, the Legendre polynomials), see e.g. [21].

The paper is organized as follows. In Section 2 we introduce the notation and basic concepts used in the following sections. In Section 3 we present the splitting/polynomial chaos expansion approach and provide a complete convergence analysis. Finally, in Section 4 we validate our approach with a numerical experiment.

2. Preliminaries

In this section we briefly recall polynomial chaos representations of random variables and stochastic processes. Particular emphasis is given to Legendre polynomials and the corresponding Wiener–Legendre expansion, and to the Karhunen–Loève expansion.

2.1. Polynomial chaos representation

Let $\mathcal{I}=(\mathbb{N}_{0}^{\mathbb{N}})_{c}$ be the set of sequences of non-negative integers which have only finitely many nonzero components $\alpha=(\alpha_{1},\alpha_{2},\ldots,\alpha_{m},0,0,\ldots)$, $\alpha_{i}\in\mathbb{N}_{0}$, $i=1,2,...,m$, $m\in\mathbb{N}$. Particularly, $(0,0,\dots)$ is the zero vector. We denote by $\varepsilon^{(k)}=(0,\cdots,0,1,0,\cdots)$, $k\in\mathbb{N}$ the $k$th unit vector. The length of $\alpha\in\mathcal{I}$ is the sum of its components $|\alpha|=\sum_{k=1}^{\infty}\alpha_{k}$.

First, we briefly recall the main results from the Wiener–Itô chaos expansion. Let $(\Omega,\mathcal{F},\mu)$ be a probability space with the Gaussian probability measure $\mu$ and let $(L)^{2}=L^{2}(\Omega,\mathcal{F},\mu)$ denote the space of random variables with finite second moments on the probability space $(\Omega,\mathcal{F},\mu)$. The space $(L)^{2}$ is a Hilbert space. The scalar product of two random variables $F,G\in(L)^{2}$ is given by

[TABLE]

where $\mathbb{E}$ denotes the expectation with respect to the measure $\mu$.

Let $\{h_{n}\}_{n\in\mathbb{N}_{0}}$ be the Hermite polynomials given through the recursion

[TABLE]

Define the $\alpha$th Fourier–Hermite polynomial as the product

[TABLE]

represented in terms of the Hermite polynomials evaluated at appropriate components of the sequence $\boldsymbol{\xi}=(\xi_{1},\xi_{2},...)$ of independent Gaussian variables with zero mean and unit variance. Especially,

[TABLE]

Theorem 2.1** (Wiener–Itô chaos expansion theorem, [20]).**

Each square integrable random variable $F\in(L)^{2}$ can be uniquely represented in the form

[TABLE]

where $f_{\alpha}\in\mathbb{R}$ for $\alpha\in\mathcal{I}$. Moreover, it holds

[TABLE]

The family of stochastic polynomials $\{H_{\alpha}\}_{\alpha\in\mathcal{I}}$ forms an orthogonal basis of $(L)^{2}$ such that

[TABLE]

for all $\alpha,\beta\in\mathcal{I}$, see [20]. Here $\delta_{\alpha\beta}$ denotes the Kronecker delta. Thus, the sequence of the coefficients in (2.1), which is a sequence of real numbers, is obtained from $f_{\alpha}=\frac{1}{\alpha!}\,\mathbb{E}(F\,H_{\alpha})$, $\alpha\in\mathcal{I}$. Also, we have

[TABLE]

Property (2.2) is a consequence of the orthogonality of the Hermite polynomials

[TABLE]

for all $m,n\in\mathbb{N}$.

In [39] it was shown that the initial construction of the Wiener chaos which corresponds to the Gaussian measure and Hermite polynomials can be extended also to other types of measures, where instead of the Hermite polynomials other classes of orthogonal polynomials from the Askey scheme [36] are used. For example, the Gamma distribution corresponds to the Laguerre polynomials and thus to the Wiener–Laguerre chaos, while the Beta distribution is related to the Jacobi polynomials and thus to the Wiener–Jacobi chaos etc. Moreover, in [36] it was proven that the optimal exponential convergence rate for each Wiener-Askey chaos can be realized.

In this paper, we deal with stochastic evolution problems with non-Gaussian random inputs which are uniformly distributed. From the Askey scheme of orthogonal polynomials it follows that the uniform distribution, as a special case of the Beta distribution, corresponds to the special class of the Jacobi polynomials, the Legendre polynomials. Therefore, we are going to work with the Wiener–Legendre polynomial chaos.

2.2. Wiener–Legendre chaos representation

Denote by $\{p_{n}(x)\}_{n\in\mathbb{N}_{0}}$ the Legendre polynomials on $[-1,1]$. These polynomials are defined by the recursion

[TABLE]

They can be also obtained from Rodrigues’ formula [36]

[TABLE]

The Legendre polynomials satisfy the second order differential equation $(1-x^{2})p_{n}^{\prime\prime}(x)-2xp_{n}^{\prime}(x)+n(n+1)p_{n}(x)=0$, which appears in physics when solving the Laplace equation in spherical coordinates [36]. These polynomials are orthogonal and it holds

[TABLE]

The previous property (2.4) is equivalent to the orthogonality relation with respect to the uniform measure, i.e., the measure with the constant weighting function $w(x)=\frac{1}{2}$.

We consider square integrable random variables and stochastic processes on a probability space $(\Omega,\mathcal{F},\mathbb{P})$ with the measure $\mathbb{P}$ generated by the uniform distribution. Let $(L)^{2}=L^{2}(\Omega,\mathcal{F},\mathbb{P})$ be the Hilbert space of square integrable random variables with respect to the measure $\mathbb{P}$.

We define the $\alpha$th Fourier–Legendre polynomial as the product

[TABLE]

where $\{p_{n}\}_{n\in\mathbb{N}_{0}}$ are the Legendre polynomials and $\boldsymbol{\xi}=(\xi_{1},\xi_{2},...)$ is a sequence of independent uniformly distributed random variables with zero mean and unit variance. Note that the product in (2.5) is finite since each $\alpha\in\mathcal{I}$ has only finitely many nonzero components. Particularly,

[TABLE]

We also have

[TABLE]

since $\boldsymbol{\xi}(\omega)$ has zero mean. Moreover, from the orthogonality (2.4) of the Legendre polynomials we obtain that the family of the Fourier–Legendre polynomials $\{L_{\alpha}\}_{\alpha\in\mathcal{I}}$ is also orthogonal and

[TABLE]

for all $\alpha,\beta\in\mathcal{I}$.

Now we formulate the representation of a random variable in an analoguous way to Theorem 2.1.

Theorem 2.2** (Wiener–Legendre chaos expansion theorem).**

Each random variable $F\in(L)^{2}$ can be uniquely represented in the form

[TABLE]

where

[TABLE]

is the corresponding sequence of real coefficients. Moreover, it holds

[TABLE]

Remark 2.1*.*

We note here that the chaos representation (2.8) of a random variable with finite second moment with respect to the underlying probability measure $\mathbb{P}$ can be extended also to square integrable stochastic processes, where a family of real numbers $f_{\alpha}$ is replaced by an appropriate family of functions with values in a certain Banach space $X$. Particularly, an $X$-valued square integrable process $u=u(t,x,\omega)$ can be represented as

[TABLE]

In this context, the notation $u\in C([0,T],X)\otimes(L)^{2}$ means that the coefficients of the process $u$ given in the form (2.9) satisfy $u_{\alpha}\in C([0,T],X)$ for all $\alpha\in\mathcal{I}$. Additionally, the estimate

[TABLE]

holds, where the expectation $\mathbb{E}L_{\alpha}^{2}$ is given by (2.7). Similarly, a process $u\in C^{1}([0,T],X)\otimes(L)^{2}$ can be represented in the form (2.9), where its coefficients $u_{\alpha}\in C^{1}([0,T],X)$ for all $\alpha\in\mathcal{I}$. Moreover, it holds

[TABLE]

2.3. Karhunen–Loève expansion

The Karhunen–Loève expansion gives a way to represent a stochastic process as an infinite linear combination of orthogonal functions on a bounded interval. It is used to represent spatially varying random inputs in stochastic models. Various applications of the Karhunen–Loève expansion can be found in uncertainty propagation through dynamical systems with random parameter functions [7], [11], [24].

Theorem 2.3** (Karhunen–Loève expansion theorem, [11]).**

Let $v(x,\omega)$ be a spatially varying square integrable random field defined over the spatial domain ${\rm D}$ and a given probability space $(\Omega,\mathcal{F},\mathbb{P})$, with mean $\bar{v}(x)$ and continuous covariance function $C_{v}(x_{1},x_{2})$. Then, $v(x,\omega)$ can be represented in the form

[TABLE]

where $\lambda_{k}$ and $e_{k}$, $k\in\mathbb{N}$ are the eigenvalues and eigenfunctions of the covariance function, i.e., they solve the integral equation

[TABLE]

and $Z_{k}$ are uncorrelated zero mean random variables that have unit variance.

For some particular covariance functions $C_{v}$, the eigenpairs $(\lambda_{k},e_{k})_{k\in\mathbb{N}}$ are known a priory, and the eigenvalues $\lambda_{k}$ decay as $k$ increases. In general, the eigenvalues and eigenvectors of the covariance function have to be calculated numerically, i.e., by solving the discrete version of (2.11). This constitutes the bottleneck of the method as it requires a large number of calculations.

In practical applications, the series are truncated, i.e., the random field is approximated by

[TABLE]

which is the finite representation with the minimal mean square error over all such finite representations.

Remark 2.2*.*

Comparing the representation (2.10) with the form (2.8) we conclude that the random field $v$ is represented in terms of the Wiener–Askey polynomial chaos of orders zero and one, i.e., it is equivalent to the representation

[TABLE]

since $Z_{k}(\omega)=L_{\varepsilon^{(k)}}({\mathbf{Z}}(\omega))$, $k\in\mathbb{N}$ with ${\mathbf{Z}}(\omega)=(Z_{1}(\omega),Z_{2}(\omega),\dots)$ being a sequence of uncorrelated uniformly distributed zero mean random variables that have unit variance. The truncated version of the representation (2.13) is given by

[TABLE]

There, $n$ corresponds to the finite number of random variables of the sequence $\mathbf{Z}=(Z_{1},Z_{2},\dots Z_{n})$ that are applied in the approximation. This is used in Section 4.

More details on methods based on stochastic polynomial representations can be found, for example, in [1], [6], [11], [24], [39].

3. Splitting methods for SPDEs

In this section, we introduce a new numerical method which combines the Wiener–Askey polynomial chaos expansion [39] with deterministic splitting methods [10]. The method is then applied to problem (1.1) with non-Gaussian random inputs. First, we are going to state a theorem on the existence and uniqueness of the solution of (1.2). Then, we recall some convergence results of splitting methods in the deterministic setting. Finally, we provide a convergence analysis of our approach which is the main result of this section. Thorough this section we denote $\mathcal{L}=A+B$.

3.1. Existence and uniqueness of the solution

Recall that a solution of the considered stochastic evolution problem (1.2) belongs to the space of square integrable stochastic processes whose coefficients are continuously differentiable deterministic functions with values in $X$.

Definition 3.1**.**

A process $u$ is a (classical) solution of (1.2) if $u\in C([0,T],X)\otimes(L)^{2}\,\cap\,C^{1}((0,T],X)\otimes(L)^{2}$ and if $u$ satisfies (1.2) pointwise.

Let the following assumptions hold:

1. $(A1)$

Let $\mathcal{L}$ be a coordinatewise operator defined on some domain $\mathcal{D}(\mathcal{L})$ dense in $X$, i.e.,

[TABLE]

for $u$ of the form (2.9). Moreover, let $\mathcal{L}$ be the infinitesimal generator of a $C_{0}$ semigroup $(S_{t})_{t\geq 0}$ of type $(M,w)$, i.e.,

[TABLE]

for some $M>0$ and $w\in\mathbb{R}$. 2. (A2)

Let $u^{0}\in X\otimes(L)^{2}$ and $\mathcal{L}u^{0}\in X\otimes(L)^{2}$, i.e.,

[TABLE] 3. (A3)

The noise process is given in the form $G(t,x,\omega)=\sum\limits_{\alpha\in\mathcal{I}}g_{\alpha}(t,x)\,L_{\alpha}\in C^{1}([0,T],X)\otimes(L)^{2}$, i.e., it holds

[TABLE]

We note here that the derivative is a coordinatewise operator, i.e., for a process $u\in C^{1}([0,T],X)\otimes(L)^{2}$ it holds

[TABLE]

Theorem 3.1** (Existence and uniqueness of the solution).**

If the assumptions $(A1)$-$(A3)$ hold, then the stochastic Cauchy problem

[TABLE]

has a unique solution

[TABLE]

in $C^{1}([0,T],X)\otimes(L)^{2}$.

Proof.

We present the main steps of the proof. We are looking for a solution in chaos representation form

[TABLE]

Then, by applying the chaos expansion method, the stochastic equation (3.1) is transformed to the infinite system of deterministic problems

[TABLE]

for all $\alpha\in\mathcal{I}$ that can be solved in parallel. Since $g_{\alpha}\in C^{1}([0,T],X)$ the inhomogeneous initial value problem (3.3) has a solution $u_{\alpha}(t)\in C^{1}((0,T],X)$ for all $\alpha\in\mathcal{I}$. Moreover, the solution $u_{\alpha}$ is given by

[TABLE]

see [35]. Thus, for all fixed $\alpha\in\mathcal{I}$ the solution $u_{\alpha}(t)$ exists for all $t\in[0,T]$, and it is a unique classical solution on the whole interval $[0,T]$. Also,

[TABLE]

Moreover, the series $\sum_{\alpha\in\mathcal{I}}u_{\alpha}(t)\,L_{\alpha}$ converges in $C^{1}([0,T],X)\otimes(L)^{2}$. Namely, from the assumptions $(A1)$-$(A3)$ we obtain

[TABLE]

where $c=c(M,w,T)$ is a constant depending on $M,w$ and $T$. ∎

Remark 3.1*.*

If an operator $A$ is the infinitesimal generator of a $C_{0}$ semigroup and $B$ is a bounded operator then the operator $\mathcal{L}=A+B$ is also the infinitesimal generator of a $C_{0}$ semigroup and Theorem 3.1 holds. In particular, Theorem 3.1 also holds for analytic semigroups.

3.2. Splitting methods for deterministic problems

We briefly recall the convergence of two operator resolvent splitting methods: resolvent Lie splitting (a first-order method) and trapezoidal resolvent splitting (a second-order method). Resolvent splitting methods for the time integration of abstract evolution equations were studied in [17]. The convergence properties of splitting methods for inhomogeneous evolution equations were analyzed in [34]. Other splitting methods were also considered in the literature. For example, exponential splitting methods for homogeneous problems with unbounded operators were presented in [14], [16]. The inhomogeneous case was studied in [10]. Error bounds for exponential operator splittings were further discussed in [23].

3.2.1. Analytic setting

Let $X$ be an arbitrary Hilbert space with norm denoted by $\|\cdot\|$. Let $X^{*}$ be the dual space of $X$. For $t\in[0,T]$ we consider the inhomogeneous evolution equation

[TABLE]

where $(\mathcal{D}(\mathcal{L}),\mathcal{L})$, $(\mathcal{D}(A),A)$ and $(\mathcal{D}(B),B)$ are linear unbounded operators in $X$ such that $\mathcal{D}(\mathcal{L})\subseteq\mathcal{D}(A)\cap\mathcal{D}(B)$ and $g:[0,T]\to X$. We recall the main results from [17] and [34].

Let the following assumptions hold:

1. (a1)

The operators $(\mathcal{D}(\mathcal{L}),\mathcal{L})$, $(\mathcal{D}(A),A)$ and $(\mathcal{D}(B),B)$ are maximal dissipative and densely defined in $X$. 2. (a2)

$\mathcal{D}(\mathcal{L}^{2})\subseteq\mathcal{D}(AB)$ 3. (a3)

Let $0\in\rho(\mathcal{L})$, let $\mathcal{L}^{-1}g(t)\in\mathcal{D}(AB)$ for all $t\in[0,T]$ and

[TABLE]

with a moderate constant $c$.

Recall that an operator $(\mathcal{D}(G),G)$ is maximal dissipative in $X$ if the following conditions hold:

• (i)

for every $x\in\mathcal{D}(G)$ there exists an element $f\in F(x)=\{h\in X^{*}:h(x)=\|x\|^{2}=\|h\|^{2}\}\subseteq X^{*}$ such that $\text{Re\;}f(Gx)\leq 0$ and

• (ii)

$\text{range\;}(I-G)=X$.

Since we assumed that $X$ is a Hilbert space, every maximal dissipative operator in $X$ is densely defined. The assumption $(a1)$ is equivalent to claiming that the operators generate $C_{0}$ semigroups of contractions on $X$, see [35]. Additionally, from $(a1)$ the following estimates hold

[TABLE]

We recall briefly the results from regularity theory for analytic semigroups needed in the following sections.

Theorem 3.2**.**

([31])* Let $\mathcal{L}$ be the generator of an analytic semigroup and let the data of problem (3.4) satisfy*

[TABLE]

for some $\theta>0$. Then, the exact solution of problem (3.4) is given by the variation of constants formula

[TABLE]

It possesses the regularity

[TABLE]

The same regularity is obtained if $g$ is only continuous but has a slightly improved spatial regularity, see [31, Corollary 4.3.9].

Theorem 3.3**.**

([34])* Let $\mathcal{L}$ be the generator of an analytic semigroup. Under the further assumptions*

[TABLE]

for some $\theta>0$, the solution (3.5) of the evolution equation (3.4) possesses the improved regularity

[TABLE]

In the following we present two deterministic resolvent splitting methods [22], the resolvent Lie splitting and the resolvent trapezoidal splitting, that were both applied to inhomogeneous evolution equations (3.4) in [34].

3.2.2. Resolvent Lie splitting

The exact solution of the evolution equation (3.4) is given by the variation of constants formula (3.5). Then, at time $t_{n+1}=t_{n}+h$, with a positive step size $h$, the solution can be written as

[TABLE]

After expanding $g(t_{n}+s)$ in Taylor form we obtain

[TABLE]

see [34]. For resolvent Lie splitting, the numerical solution of (3.4) at time $t_{n+1}$ is denoted by $u^{n+1}$ and it is given by

[TABLE]

Theorem 3.4** (Resolvent Lie splitting, [34]).**

Let the assumptions $(a1)$, $(a2)$ and $(a3)$ be fulfilled and let the solution satisfy (3.7). Then the resolvent Lie splitting (3.8) is first-order convergent, i.e., the global error satisfies the bound

[TABLE]

with a constant $C$ that can be chosen uniformly on $[0,T]$ and, in particular, independently of $n$ and $h$.

Remark 3.2*.*

The constant $C$ in (3.9) depends on derivatives of the solution $u$ and on $AB\mathcal{L}^{-1}\,g(t)$, which are uniformly bounded on $[0,T]$ due to the asumptions of Theorem 3.4. A detailed proof is given in [34].

In particular, for a homogeneous evolution problem ($g=0$) the global error (3.9) can be estimated as

[TABLE]

where the positive constant $c$ is independent on $n$ and $h$, see [17].

We note that the full-order convergence of Lie resolvent splitting only requires additional smoothness in space of the inhomogeneity $g$.

3.2.3. The trapezoidal splitting

For a trapezoidal splitting method, the numerical solution of (3.4) at time $t_{n+1}=t_{n}+h$ with a positive time step size $h$ is given by

[TABLE]

with $u^{0}=u(0)$.

As we are considering a second-order method, we need more regularity of the solution. For analytic semigroups, this requirement can be expressed in terms of the data. The following modification of the assumption $(a3)$ is needed:

1. (a4)

Let $0\in\rho(\mathcal{L})$, let $\mathcal{L}^{-1}g^{\prime}(t)\in\mathcal{D}(AB)$ for all $t\in[0,T]$ and

[TABLE]

with a moderate constant $c$.

Since we assumed $X$ to be a Hilbert space, it follows from assumption (a1) that the estimates

[TABLE]

hold for all $h>0$.

Theorem 3.5**.**

([34])*

Let $\mathcal{L}$ be the generator of an analytic semigroup. If*

[TABLE]

for some $\theta>0$, then the exact solution (3.5) of the inhomogeneous evolution equation (3.4) satisfies

[TABLE]

Theorem 3.6**.**

(The trapezoidal splitting method, [34])*

Let the assumptions $(a1)$, $(a2)$ and $(a4)$ be fulfilled and let the solution satisfy (3.12). Then the trapezoidal splitting method (3.10) is second-order convergent, i.e., the global error satisfies the bound*

[TABLE]

with a constant $C$ that can be chosen uniformly on $[0,T]$ and, in particular, independently of $n$ and $h$.

Remark 3.3*.*

The constant $C$ in (3.13) depends on derivatives of the solution $u$ and on $AB\mathcal{L}^{-1}\,g^{\prime}(t)$, which are uniformly bounded on $[0,T]$ due to the asumptions of Theorem 3.6. More details are given in [34].

3.3. Convergence analysis

In order to solve problem (3.1) numerically, we approximate the solution $u$ by the truncated chaos representation form

[TABLE]

where $\mathcal{I}_{m,K}=\{\alpha\in\mathcal{I}:\,\alpha=(\alpha_{1},\dots,\alpha_{m},0,0,\dots),\,|\alpha|\leq K\}$. Here, $K\in\mathbb{N}$ is the highest degree of Legendre polynomials and $m\in\mathbb{N}$ is the number of random variables we want to use in the approximation (3.14). The $m$-dimensional random vector $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{m})$ has independent and identically distributed components $\xi_{i}\sim\mathcal{U}([-1,1])$ for $i=1,\ldots,m$. The choice of $m$ and $K$ influences the accuracy of the approximation. They can be chosen so that the norm of the approximation remainder $u-\tilde{u}$ is smaller than a given tolerance. The sum in (3.14) has

[TABLE]

terms, which means that $P$ coefficients of the solution will be computed. Thus, only the first $P$ equations of the system (3.3) are solved and in this way the approximation of the solution of the system is obtained. The global error of the proposed numerical scheme depends on the error generated by the truncation of the chaos expansion and the error of the discretisation method. Also, the statistics $\mathbb{E}\tilde{u}$ and $\operatorname{Var}\tilde{u}$ of the approximated solution can be calculated in terms of the obtained discretized coefficients. For more details on the truncation (3.14) see for instance [39]. In the following, we consider the two numerical resolvent splitting methods, Lie splitting and trapezoidal splitting, and provide error analysis for both of them.

Theorem 3.7** (Error generated by the truncation of the Wiener–Legendre chaos expansion).**

Let $\tilde{u}$ denote the truncated chaos representation of the solution $u$ of the stochastic evolution problem (3.1) given in the form (3.14). Let the assumptions $(A1)$-$(A3)$ hold. Then, $\tilde{u}$ approximates the solution $u$ and the approximation error satisfies the a priori bound

[TABLE]

Proof.

The approximation error due to the elimination of the higher order components of the Wiener–Legendre chaos expansion and the truncation of the noise term is obtained by

[TABLE]

which is finite by the assumptions (A1)-(A3). In the last estimate, we employed the bound derived in the proof of Theorem 3.1. ∎

Theorem 3.8** (Discretization error).**

Let $\tilde{u}$ denote the truncated chaos representation of the solution $u$ of the stochastic evolution problem (3.1) given in the form (3.14). Let a square integrable process $\tilde{u}^{n}_{dis}$ be given in the form

[TABLE]

where its coefficients ${u}^{n}_{\alpha,dis}$, $\alpha\in\mathcal{I}_{m,K}$ are numerical approximations of $u_{\alpha}$ for $\alpha\in\mathcal{I}_{m,K}$ at time $t_{n}=nh$ with a positive step size $h$. Assume that the coefficients $u_{\alpha}$ are sufficiently regular and the approximation

[TABLE]

holds for the particular numerical method applied. Then, the difference between $\tilde{u}$ evaluated at $t_{n}$ and $\tilde{u}^{n}_{dis}$ can be estimated by the a priori bound

[TABLE]

Proof.

From Parseval’s identity and the orthogonality of the polynomial basis $\{L_{\alpha}\}$, and using that the error (3.17) for a concrete numerical method, we obtain

[TABLE]

which completes the proof. ∎

In order to apply the splitting methods in the setting of [34], we are going to consider the analytic case and adapt Theorem 3.8. We replace the assumption $(A1)$ with the assumption:

1. (B1)

Let $(A,\mathcal{D}(A))$, $(B,\mathcal{D}(B))$ and $(\mathcal{L},\mathcal{D}(\mathcal{L}))$ be coordinatewise operators that generate analytic semigroups of contractions on $X$. Let $\mathcal{D}(\mathcal{L}^{2})\subseteq\mathcal{D}(AB)$.

Further, for the case of the resolvent Lie splitting we replace the assumptions $(A2)$ and $(A3)$ by:

1. (B2)

The noise process given by

[TABLE]

belongs to $C^{1+\theta}([0,T],X)\otimes(L)^{2}$ for some $\theta>0$, i.e.,

[TABLE]

holds. 2. (B3)

Let $u^{0}\in\mathcal{D}(\mathcal{L})\otimes(L)^{2}$ and $\mathcal{L}u^{0}+G(0)\in\mathcal{D}(\mathcal{L})\otimes(L)^{2}$, i.e.,

[TABLE] 3. (B4)

Let $0\in\rho(\mathcal{L})$, let $\mathcal{L}^{-1}G(t)\in\mathcal{D}(AB)\otimes(L)^{2}$ for all $t\in[0,T]$ and let the coefficients $g_{\alpha}$ of $G$ given by (3.18), satisfy the estimate

[TABLE]

with a moderate constant $c_{\alpha}$ for each $\alpha\in\mathcal{I}$.

Note that, under these assumptions, the existence theorem, Theorem 3.1, still holds. Particularly, for the resolvent Lie splitting it reads:

Theorem 3.9**.**

Let $\mathcal{L}$ be the generator of an analytic semigroup. Under the assumptions $(B2)$ and $(B3)$, the solution (3.2) of the stochastic evolution problem (1.2) posseses the improved regularity

[TABLE]

Proof.

By the method of chaos expansion, the stochastic evolution problem (1.2) transforms to the system of deterministic problems (3.3). From $(B2)$ and $(B3)$ it follows that $u_{\alpha}^{0}$ and $g_{\alpha}$ for each $\alpha\in\mathcal{I}$ satisfy the assumptions (3.6). After applying Theorem 3.3 we obtain the improved regularity $u_{\alpha}\in C^{2}([0,T],X)\cap C^{1}([0,T],D(\mathcal{L}))$, $\alpha\in\mathcal{I}$. ∎

Theorem 3.10** (Discretization error, the resolvent Lie splitting).**

Let the assumptions $(B1)$-$(B4)$ be fulfilled. Then, for the resolvent Lie splitting, Theorem 3.8 holds with

[TABLE]

The constants $c_{\alpha}$ can be chosen uniformly on $[0,T]$ and, in particular, independently of $n$ and $h$.

Proof.

The coefficients $u_{\alpha}$, for each $\alpha\in\mathcal{I}_{m,K}$ are the exact solutions of the deterministic initial value problems (3.3) and $u^{n}_{\alpha,dis}$ are their numerical approximations obtained by the resolvent Lie splitting (3.8). Moreover, $u_{\alpha}$ satisfy the assumptions (3.6) for all $\alpha\in\mathcal{I}$. Thus, we can apply Theorem 3.4 to each initial value problem (3.3) and obtain the global estimate (3.9) for each $\alpha\in\mathcal{I}_{m,K}$, i.e. $e_{\alpha}\leq c_{\alpha}h$, for $\alpha\in\mathcal{I}_{m,K}$. This leads to the desired result. ∎

In the case of the trapezoidal resolvent splitting, we need the following additional assumptions:

1. (B5)

The noise process $G$ given by (3.18) belongs to $C^{2+\theta}([0,T],X)\otimes(L)^{2}$ for some $\theta>0$. 2. (B6)

Let $\mathcal{L}^{2}u^{0}+\mathcal{L}G(0)+G^{\prime}(0)\in\mathcal{D}(\mathcal{L})\otimes(L)^{2}$, i.e.,

[TABLE] 3. (B7)

Let $0\in\rho(\mathcal{L})$, let $\mathcal{L}^{-1}G^{\prime}(t)\in\mathcal{D}(AB)\otimes(L)^{2}$ for all $t\in[0,T]$ and let the coefficients $g_{\alpha}$ of $G$ given by (3.18), satisfy the estimate

[TABLE]

with a moderate constant $c_{\alpha}$ for each $\alpha\in\mathcal{I}$.

Theorem 3.11**.**

Let $\mathcal{L}$ be the generator of an analytic semigroup. Under the assumptions $(B3)$, $(B5)$ and $(B6)$, the solution (3.2) of the stochastic evolution problem (1.2) posseses the improved regularity

[TABLE]

Proof.

The method of chaos expansion transforms the stochastic evolution problem (1.2) to the system of deterministic problems (3.3). From $(B3)$, $(B5)$ and $(B6)$ it follows that $u_{\alpha}^{0}$ and $g_{\alpha}$ for each $\alpha\in\mathcal{I}$ satisfy the assumptions (3.11). Then, the improved regularity $u_{\alpha}\in C^{3}([0,T],X)\cap C^{2}([0,T],D(\mathcal{L}))$, for $\alpha\in\mathcal{I}$ follows from Theorem 3.5. ∎

Theorem 3.12** (Discretization error, the trapezoidal resolvent splitting).**

Let the assumptions $(B1)$, $(B3)$ and $(B5)$-$(B7)$ be fulfilled. Then, for the trapezoidal resolvent splitting, Theorem 3.8 holds with

[TABLE]

The constants $c_{\alpha}$ can be chosen uniformly on $[0,T]$ and, in particular, independently of $n$ and $h$.

Proof.

From the assumptions it follows that the coefficients $u^{0}_{\alpha}$ and $g_{\alpha}$ satisfy (3.11) for each $\alpha\in\mathcal{I}_{m,K}$. We apply the trapezoidal resolvent splitting (3.10) in order to obtain the approximation $u^{n}_{\alpha,dis}$ of the exact solution $u_{\alpha}(t_{n})$ evaluated at $t_{n}$ of the initial value problem (3.3) for each $\alpha\in\mathcal{I}_{m,K}$. Thus, by Theorem 3.6 we obtain the global error estimate (3.13), i.e. $e_{\alpha}\leq c\,h^{2}$ for each $\alpha\in\mathcal{I}_{m,K}$. ∎

Denote by $\frac{1}{2}\Delta$ the constant on the right hand side of the estimate (3.16) obtained in Theorem 3.7. The full error estimates of the Wiener–Legendre chaos expansion combined with the two splitting methods are given in the following theorem.

Theorem 3.13** (Full error estimate).**

**

1. $(1)$

Let the assumptions of Theorem 3.10 hold. Then, the full error estimate of the Wiener–Legendre chaos expansion combined with the resolvent Lie splitting satisfies the following bound

[TABLE] 2. $(2)$

Let the assumptions of Theorem 3.12 hold. Then, the full error estimate of the Wiener–Legendre chaos expansion combined with the trapezoidal resolvent splitting satisfies the bound

[TABLE]

Proof.

The full error estimate reads

[TABLE]

by the triangle inequality and the orthogonality property (2.7). We apply Theorem 3.7 to the first term. In the case of the resolvent Lie splitting, the estimate (3.21) follows after applying Theorem 3.10, while in case of the trapezoidal resolvent splitting, Theorem 3.12 leads to the desired estimate (3.22). ∎

4. Numerical Results

In this section, we validate the proposed method and the convergence analysis presented in the previous section. For this purpose, we consider the two-dimensional problem

[TABLE]

where the operator $\mathcal{L}$ is defined by $\mathcal{L}u=(A+B)u=(au_{x})_{x}+(bu_{y})_{y}$ over the spatial domain ${\rm D}=[-1,1]^{2}$ with state variables $x$ and $y$, spatial non-Gaussian noise $v$ given in the form (2.10) and $t\in[0,T]$ for some $T>0$. This problem is an example of the problem class (1.2) with zero initial and boundary conditions. The solution $u$ of the considered problem (4.1) is given in its polynomial chaos representation (2.9) and approximated by a truncated expansion (3.14) in terms of Fourier–Legendre polynomials. The truncation procedure is explained in detail in Section 3.3.

Consider the set of multiindices $\mathcal{I}_{m,K}\subset\mathcal{I}$, i.e.,

[TABLE]

In this section, elements $\alpha\in\mathcal{I}_{m,K}$ will be denoted as $m$-tuples $\alpha=(\alpha_{1},\dots,\alpha_{m})$, omitting the components $\alpha_{j}=0$, $j\geq m+1$. Moreover, we set

[TABLE]

For fixed $m\in\mathbb{N}$ we consider an index function

[TABLE]

which enumerates multi-indices $\alpha=(\alpha_{1},\alpha_{2},\dots,\alpha_{m})\in\mathcal{I}_{m,K}$. The function $K_{m}$ is a bijection and each $\alpha\in\mathcal{I}_{m,K}$ corresponds to a unique $K_{m}(\alpha)=p\in\{0,1,\dots P-1\}$. For our purpose, we define the function $K_{m}$ by

[TABLE]

We use the index function $K_{m}$ to enumerate the Fourier–Legendre polynomials $L_{\alpha}$ for each $\alpha\in\mathcal{I}_{m,K}$. Thus, we denote by $(\Phi_{p})_{p\in\{0,1,\dots,P-1\}}$ the ordered Fourier–Legendre polynomials

[TABLE]

for $p=K_{m}(\alpha)$, $\alpha\in\mathcal{I}_{m,K}$, where we use the definition (2.5) of the Fourier–Legendre polynomials. For example, following the just introduced notation, we have $\Phi_{0}(\boldsymbol{\xi}(\omega))=L_{(0,0,\dots,0)}(\boldsymbol{\xi}(\omega))=1$ and

[TABLE]

Also, by applying the definition of the Legendre polynomials (2.3) we have

[TABLE]

as well as

[TABLE]

Moreover, it holds

[TABLE]

In the next step, we represent the solution $u$ of problem (4.1) by its truncated polynomial chaos expansion (3.14) and the noise term by its representation (2.14). Inserting the representations in (4.1) gives

[TABLE]

By performing a Galerkin projection we obtain

[TABLE]

for $\beta\in\mathcal{I}_{m,K}$.

Then, by applying the properties of the Fourier–Legendre polynomials (2.6) and (2.7), we obtain a system of deterministic equations (3.3). Particularly,

1. (i)

for $|\alpha|=0$:

[TABLE] 2. (ii)

for $|\alpha|=1$, i.e., $\alpha=\varepsilon^{(k)}$, $1\leq k\leq m$:

[TABLE] 3. (iii)

for $|\alpha|>1$:

[TABLE]

From (4.4) we clearly deduce that $u_{\alpha}\equiv 0$ for $|\alpha|>1$. In the calculations we also used $\mathbb{E}(L_{(0,0,\ldots,0)}Z_{j})=\mathbb{E}Z_{j}=0$ for $j\geq 1$ and

[TABLE]

This particularly implies

[TABLE]

which was used in equation (4.3).

The obtained system (4.2), (4.3) and (4.4) can be represented in terms of the index function $K_{m}$, i.e., in the form

[TABLE]

for $0\leq p\leq P-1$, where each $p$ corresponds to an $\alpha\in\mathcal{I}_{m,K}$ Each equation in (4.5) has the form of an inhomogeneous deterministic initial value problem, where the inhomogeneities $g_{p}$ are given by: $g_{0}=\bar{v}+1$ and $g_{p}=\sqrt{\lambda_{p}}\,e_{p}$ for $1\leq p\leq m$ and $g_{p}=0$ for $m<p\leq P-1$.

One way to approximate numerically a problem of the form

[TABLE]

with ${\rm D}=[-1,1]^{2}$ is to define a grid consisting of $N\times N$ equidistant computational points and define the discrete operators ${A}_{s}$ and ${B}_{s}$ by

[TABLE]

where

[TABLE]

for $i,j=1,\ldots,N$ and $s=2/(N+1)$. Due to the homogeneous Dirichlet boundary conditions we have:

[TABLE]

for all $i,j=0,\ldots,N+1$. By setting $\mathcal{L}_{s}=A_{s}+B_{s}$ we obtain the discretized problem

[TABLE]

where $g_{s}$ denotes the discretization of the inhomogeneity $g$.

Note that the number $P$ of partial differential equations one has to solve in (4.5) increases fast due to the factorials occurring in (3.15). Since $g_{p}=0$ for all $m<p\leq P-1$, $u_{p}=0$ is consequently the solution of the $p$th partial differential equation of (4.5). Therefore, we only have to solve the first $m+1$ partial differential equations instead of all $P$. Further, we see that the solution does not depend on the highest degree $K$ of the $m$-dimensional Legendre polynomials.

Let $u_{p}^{n}$ denote the numerical solution $u_{p}$ at time $t_{n}=hn$ and $g_{p}^{n}$ the function $g_{p}$ evaluated at time $t_{n}$. By setting

[TABLE]

the Lie resolvent splitting method is defined, see (3.8).

The trapezoidal splitting method is given by

[TABLE]

see (3.10).

In our numerical experiment, we consider (4.1) with constant coefficients $a(x,y)=b(x,y)=1$ for all $(x,y)\in{\rm D}=[-1,1]^{2}$ and set $T=1$. Note that for some $p\in\{0,\ldots,m\}$ the inhomogeneities $g_{p}$ might be incompatible with the boundary conditions at the corners of the spatial domain ${\rm D}$. Such an incompatibility results in order reduction, see [18]. This in particular leads to large errors near the corners of ${\rm D}$. To overcome this problem, we apply the modified Lie resolvent splitting [18] in this situation.

For $p\in\{0,\ldots,m\}$, let $u_{p}$ be the solution of the partial differential equation (4.5). Let $I=\{1,2,3,4\}$ be the set of indices of the corners of the spatial domain ${\rm D}$. They are enumerated from 1 to 4 counter-clockwise starting from the corner with coordinates $(-1,-1)$. Suppose that the inhomogeneity $g_{p}$ does not vanish at the corners $I_{p}\subset I$. Let $g_{p,i}(t)$ denote the value of the function $g_{p}$ at corner $i\in I_{p}$ and time $t\geq 0$. For $g_{p,i}(0)\neq 0$ we set

[TABLE]

where the polynomials $P_{i}$ are given by

[TABLE]

These four polynomials form a partition of unity.

Let $v_{i}$ be the solution of the stationary problem

[TABLE]

for $i\in I_{p}$. Note that $v_{i}$ can be computed once and for all. Then, let

[TABLE]

and apply the resolvent Lie splitting to the problem

[TABLE]

By setting

[TABLE]

we obtain the modified splitting scheme. Note that in our case $g_{p,i}^{\prime}(t)=0$ for all $i\in I_{p}$ and for all $p=0,\ldots,m$ since none of the inhomogeneities $g_{p}$ is time dependent.

In the implementation, the set $I_{p}$ for $p=0,\ldots,m$ is constructed by checking the values of the inhomogeneities $g_{p}$ at the corners, i.e.,

[TABLE]

for a user chosen tolerance TOL. If $I_{p}=\emptyset$, the standard Lie resolvent splitting given in (4.6) is applied.

In the following, we consider problem (4.1) with $v$ given by (2.10) with covariance function

[TABLE]

The reference solution $u^{\text{ref}}_{p}$ at time $t$ is calculated according to

[TABLE]

where $\varphi_{1}(z)=\frac{\mathrm{exp}(z)-1}{z}$ and $\mathrm{exp}(\cdot)$ denotes the matrix exponential. In all the examples shown we fix the highest degree of ordered Fourier–Legendre polynomials to $K=3$ and use a maximal number of $m=120$ uncorrelated zero-mean random variables $Z_{j}$ used in the truncated Karhunen–Loève expansion (2.12). If not stated explicitly, we fix the number of computational points to $N\times N=40\times 40$.

Figure 4.1 illustrates the impact of the modification of the Lie resolvent splitting method. The figure shows the pointwise error of the numerical solution at time $T=1$, i.e., $|u_{0}(T)-u^{\text{ref}}_{0}(T)|$ over the spatial domain ${\rm D}=[-1,1]^{2}$ when calculated with the Lie splitting and the modified Lie splitting given in (4.6) and (4.7), respectively. The pointwise error of the solution $u_{0}$ is not only reduced at all the four corners of the domain ${\rm D}$ but also approximately decreases by an order of magnitude.

Figure 4.2 shows the discrete $L^{2}$ error of $u_{p}$, $p=0,\ldots,7$ calculated with different time step sizes $h$. The time step sizes are set to $h_{q}=2^{q}$ for $q=-13,\ldots,-4$. The blue line denotes the error of the modified Lie splitting scheme of order 1. The red line and the green line illustrate the error of the Crank–Nicolson scheme and the trapezoidal splitting method, both of order two. The black dashed lines have slope 1 and 2, respectively. We see that for each $m$, the order plots confirm the respective orders of the methods which can be derived from theory.

The empirical variance $\mathrm{Var}(u)$ of $u$ is given by

[TABLE]

where we used the linearity of $\mathbb{E}$ and the orthogonality of the Fourier–Legendre polynomials. Furthermore, since $u_{\alpha}\equiv 0$ for $|\alpha|>1$, i.e., $u_{p}\equiv 0$ for $p>m$, the number of non-zero summands in the sum is $m$ and since $\mathbb{E}(\Phi_{p}^{2})=\frac{1}{3}$ for $1<p\leq m$, $\mathrm{Var}(u)$ reduces to

[TABLE]

Figure 4.3 shows the discrete $L^{2}$ error of the empirical variance of $u$ at time $T=1$ where the summation is truncated at different $n$. The time step $h$ used for the calculations is $h=2^{-10}$. Here, we clearly see the superiority of

the methods of order two compared to the modified Lie splitting for which the numerical approximation error prevails over the error induced by the truncation of the sum.

Finally, we report the computational work which is needed to solve the system of partial differential equations given in (4.2)-(4.4). Table 4.1 summarizes the computational time needed to obtain one solution of the system of partial differential equations as a function of the number of spatial grid points $N\times N=2^{k}\times 2^{k}$ for $k=2,3,\ldots 7$. The highest number of grid points we are able to use (16 384) is quite low due to the fact that the calculation of the eigenvalues and eigenfunctions of the integral equation given in (2.11) requires the storage of a dense matrix of the size $N^{2}\times N^{2}$. We clearly see that the Crank–Nicolson method is by far the slowest. Both splitting methods perform approximately the for smaller $N$, while for $N=2^{7}$, Lie splitting starts to clearly outperform trapezoidal splitting in terms of computational time.

5. Acknowledgements

This work was partially supported by a Research grant for Austrian graduates granted by the Office of the Vice Rector for Research of University of Innsbruck. The computational results presented have been partially achieved using the HPC infrastructure LEO of the University of Innsbruck. A. Kofler was supported by the program Nachwuchsförderung 2014 at University of Innsbruck. H. Mena was supported by the Austrian Science Fund – project id: P27926.