Multilevel Methods for Uncertainty Quantification of Elliptic PDEs with Random Anisotropic Diffusion

TL;DR

This paper develops multilevel methods for efficiently quantifying uncertainty in elliptic PDEs with random anisotropic diffusion, showing how the decay of the Karhunen-Loève expansion influences solution regularity and computational complexity.

Contribution

It establishes a link between the decay of the Karhunen-Loève expansion and solution regularity, enabling multilevel methods to efficiently approximate statistical quantities.

Findings

Decay of the Karhunen-Loève expansion determines solution regularity.

Multilevel collocation and quadrature methods achieve expected convergence rates.

Numerical examples validate the theoretical results in three dimensions.

Abstract

We consider elliptic diffusion problems with a random anisotropic diffusion coefficient, where, in a notable direction given by a random vector field, the diffusion strength differs from the diffusion strength perpendicular to this notable direction. The Karhunen-Lo\`eve expansion then yields a parametrisation of the random vector field and, therefore, also of the solution of the elliptic diffusion problem. We show that, given regularity of the elliptic diffusion problem, the decay of the Karhunen-Lo\`eve expansion entirely determines the regularity of the solution's dependence on the random parameter, also when considering this higher spatial regularity. This result then implies that multilevel collocation and multilevel quadrature methods may be used to lessen the computation complexity when approximating quantities of interest, like the solution's mean or its second moment, while still yielding the expected rates of convergence. Numerical examples in three spatial dimensions are provided to validate the presented theory.

Tables

Table 1. Table 1. The number of samples for the first six levels and the respective parameter dimensions.

$l$	$0$	$1$	$2$	$3$	$4$	$5$
$N_l$	$10$	$24$	$57$	$135$	$320$	$762$
$M_1$	$17$	$26$	$30$	$36$	$44$	$52$
$M_2$	$14$	$26$	$30$	$36$	$43$	$52$

Equations

$$\mathbf{A}[\omega]\mathrel{\mathrel{\mathop{:}}=}a\mathbf{I}+\Big{(}\big{\lVert}\mathbf{V}[\omega]\big{\rVert}_{2}-a\Big{)}\frac{\mathbf{V}[\omega]\mathbf{V}^{\mathsf{T}}[\omega]}{\mathbf{V}^{\mathsf{T}}[\omega]\mathbf{V}[\omega]},$$

$$\mathbf{A}[\omega ]:=$$

$$\left[\begin{array}{ccc}\cos {\alpha }_1[\omega ]& 0& \sin {\alpha }_1[\omega ]\\ 0& 1& 0\\ -\sin {\alpha }_1[\omega ]& 0& \cos {\alpha }_1[\omega ]\end{array}\right]\left[\begin{array}{ccc}{a}_{\mathbf{f}}[\omega ]& 0& 0\\ 0& a_{\mathbf{t}}[\omega ]& 0\\ 0& 0& a_{\mathbf{s}}[\omega ]\end{array}\right]{\left[\begin{array}{ccc}\cos {\alpha }_1[\omega ]& 0& \sin {\alpha }_1[\omega ]\\ 0& 1& 0\\ -\sin {\alpha }_1[\omega ]& 0& \cos {\alpha }_1[\omega ]\end{array}\right]}^{\mathsf{T}}$$

$${\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\alpha }_2[\omega ]& \sin {\alpha }_2[\omega ]\\ 0& -\sin {\alpha }_2[\omega ]& \cos {\alpha }_2[\omega ]\end{array}\right]}^{\mathsf{T}}{\left[\begin{array}{ccc}\cos {\alpha }_3[\omega ]& -\sin {\alpha }_3[\omega ]& 0\\ \sin {\alpha }_3[\omega ]& \cos {\alpha }_3[\omega ]& 0\\ 0& 0& 1\end{array}\right]}^{\mathsf{T}}{\left[\begin{array}{ccc}|& |& |\\ \mathbf{f}& \mathbf{t}& \mathbf{s}\\ |& |& |\end{array}\right]}^{\mathsf{T}}$$

$$\text{for almost every $\omega\in\Omega$: }\left\{\begin{aligned} -\operatorname{div}_{\mathbf{x}}\big{(}\mathbf{A}[\omega]\operatorname{\nabla}_{\mathbf{x}}u[\omega]\big{)}&=f&\quad&\text{in }D,\\ u[\omega]&=0&\quad&\text{on }\partial D,\end{aligned}\right.$$

$$\lVert v\rVert_{L_{\mu}^{p}(\mathcal{M};\mathcal{X})}\mathrel{\mathrel{\mathop{:}}=}\begin{cases}\bigg{[}\displaystyle\int_{\mathcal{M}}\big{\lVert}v(x)\big{\rVert}_{\mathcal{X}}^{p}\operatorname{d}\,\mu(x)\bigg{]}^{1/p},&p<\infty,\\ \displaystyle\operatorname*{ess\,sup}_{x\in\mathcal{M}}\big{\lVert}v(x)\big{\rVert}_{\mathcal{X}},&p=\infty.\end{cases}$$

$$\mathcal{B}(\mathcal{X}_{1},\ldots,\mathcal{X}_{r};\mathcal{Y})\mathrel{\mathrel{\mathop{:}}=}\mathcal{B}\big{(}\mathcal{X}_{1};\mathcal{B}(\mathcal{X}_{2},\ldots,\mathcal{X}_{r};\mathcal{Y})\big{)}$$

$$\mathcal{B}^{0}(\mathcal{X};\mathcal{Y})\mathrel{\mathrel{\mathop{:}}=}\mathcal{Y}\quad\text{and}\quad\mathcal{B}^{r+1}(\mathcal{X};\mathcal{Y})\mathrel{\mathrel{\mathop{:}}=}\mathcal{B}\big{(}\mathcal{X};\mathcal{B}^{r}(\mathcal{X};\mathcal{Y})\big{)}.$$

$$\lvert\mathbf{v}^{\mathsf{T}}\mathbf{w}\rvert=\big{\lvert}\langle\mathbf{v},\mathbf{w}\rangle\big{\rvert}\leq\lVert\mathbf{v}\rVert_{2}\lVert\mathbf{w}\rVert_{2},$$

$$\Vert \mathbf{v}{\mathbf{w}}^{\mathsf{T}}{\Vert }_2=\Vert \mathbf{v}{\Vert }_2\Vert \mathbf{w}{\Vert }_2.$$

$$\text{for almost every $\omega\in\Omega$: }\left\{\begin{aligned} -\operatorname{div}_{\mathbf{x}}\big{(}\mathbf{A}[\omega]\operatorname{\nabla}_{\mathbf{x}}u[\omega]\big{)}&=f&\quad&\text{in }D,\\ u[\omega]&=0&\quad&\text{on }\partial D,\end{aligned}\right.$$

$$\underline{a}\leq\operatorname*{ess\,inf}_{\mathbf{x}\in D}\lambda_{\min}\big{(}\mathbf{A}[\omega](\mathbf{x})\big{)}\leq\operatorname*{ess\,sup}_{\mathbf{x}\in D}\lambda_{\max}\big{(}\mathbf{A}[\omega](\mathbf{x})\big{)}\leq\overline{a}\quad\text{$\mathbb{P}$-almost surely}$$

$$\mathbf{A}[\omega](\mathbf{x})\mathrel{\mathrel{\mathop{:}}=}a\mathbf{I}+\Big{(}\big{\lVert}\mathbf{V}[\omega](\mathbf{x})\big{\rVert}_{2}-a\Big{)}\frac{\mathbf{V}[\omega](\mathbf{x})\mathbf{V}^{\mathsf{T}}[\omega](\mathbf{x})}{\mathbf{V}^{\mathsf{T}}[\omega](\mathbf{x})\mathbf{V}[\omega](\mathbf{x})},$$

$$\underline{a}\leq\operatorname*{ess\,inf}_{\mathbf{x}\in D}\big{\lVert}\mathbf{V}[\omega](\mathbf{x})\big{\rVert}\leq\operatorname*{ess\,sup}_{\mathbf{x}\in D}\big{\lVert}\mathbf{V}[\omega](\mathbf{x})\big{\rVert}\leq\overline{a}\quad\text{$\mathbb{P}$-almost surely}$$

$$\mathbf{V}[\mathbf{y}](\mathbf{x})=E[\mathbf{V}](\mathbf{x})+\sum _{k=1}^{\infty }{\sigma }_k{\mathbf{\psi }}_k(\mathbf{x})y_k,$$

$$\text{for almost every $\mathbf{y}\in\square$: }\left\{\begin{aligned} -\operatorname{div}_{\mathbf{x}}\big{(}\mathbf{A}[\mathbf{y}]\operatorname{\nabla}_{\mathbf{x}}u[\mathbf{y}]\big{)}&=f&\quad&\text{in }D,\\ u[\mathbf{y}]&=0&\quad&\text{on }\partial D.\end{aligned}\right.$$

$$\left\{\begin{aligned} &\text{Find $u\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;H_{0}^{1}(D)\big{)}$ such that}\\ &\quad\big{(}\mathbf{A}[\mathbf{y}]\operatorname{\nabla}_{\mathbf{x}}u[\mathbf{y}],\operatorname{\nabla}_{\mathbf{x}}v\big{)}_{L^{2}(D;\mathbb{R}^{d})}=\big{(}f,v\big{)}_{L^{2}(D;\mathbb{R}^{d})}\\ &\quad\text{for almost every $\mathbf{y}\in\square$ and all $v\in H_{0}^{1}(D)$.}\end{aligned}\right.$$

$$\big{\lVert}u[\mathbf{y}]\big{\rVert}_{L_{\mathbb{P}_{\mathbf{y}}}^{\infty}(\square;H^{1}(D))}\leq\frac{1}{\underline{a}c_{V}^{2}}\Big{(}\lVert f\rVert_{H^{-1}(D)}\Big{)},$$

$$\gamma_{\kappa,k}\mathrel{\mathrel{\mathop{:}}=}\big{\lVert}\sigma_{k}{\bm{\psi}}_{k}\big{\rVert}_{W^{\kappa,\infty}(D;\mathbb{R}^{d})},$$

$$c_{{\bm{\gamma}}_{\kappa}}=\max\big{\{}\lVert{\bm{\gamma}}_{\kappa}\rVert_{\ell^{1}(\mathbb{N})},2\big{\}}.$$

$$\mathbf{V}[\mathbf{y}](\mathbf{x})=E[\mathbf{V}](\mathbf{x})+\sum _{k=1}^M{\sigma }_k{\mathbf{\psi }}_k(\mathbf{x})y_k,$$

$$\mathbf{A}\in W^{\kappa ,\infty }(D;{\mathbb{R}}^{d\times d})\cap {\mathcal{R}}_{\kappa }$$

$$\Big{(}\mathbf{A}\operatorname{\nabla}_{\mathbf{x}}u,\operatorname{\nabla}_{x}v\Big{)}_{L^{2}(D;\mathbb{R}^{d})}=(h,v)_{L^{2}(D)}$$

$$\lVert u\rVert_{\kappa+1,2,D}\leq C_{\kappa,er}\big{(}D,\lVert\mathbf{A}\rVert_{\mathcal{R}_{\kappa}}\big{)}\lVert h\rVert_{\kappa-1,2,D},$$

$$\lVert v\rVert_{\eta,p,D;\mathcal{X}}\mathrel{\mathrel{\mathop{:}}=}\lVert v\rVert_{W^{\eta,p}(D;\mathcal{X})}\mathrel{\mathrel{\mathop{:}}=}\sum_{\lvert{\bm{\alpha}}\rvert\leq\eta}\frac{1}{{\bm{\alpha}}!}\big{\lVert}\operatorname{\partial}_{\mathbf{x}}^{\bm{\alpha}}v\big{\rVert}_{p,D;\mathcal{X}},$$

$$\Vert \cdot {\Vert }_{p,D;\mathcal{X}}:=\Vert \cdot {\Vert }_{L^p(D;\mathcal{X})}.$$

$$\Vert {\mathrm{div}}_{\mathbf{x}}\mathbf{v}{\Vert }_{\eta -1,p,D}\le \eta d\Vert \mathbf{v}{\Vert }_{\eta ,p,D}$$

$$\Vert D_{\mathbf{x}}\mathbf{v}{\Vert }_{\eta -1,p,D}\le \eta d\Vert \mathbf{v}{\Vert }_{\eta ,p,D}.$$

$$\Vert {\nabla }_{\mathbf{x}}v{\Vert }_{\eta -1,p,D}\le \eta d\Vert v{\Vert }_{\eta ,p,D}.$$

$$\Vert {\mathrm{div}}_{\mathbf{x}}\mathbf{v}{\Vert }_{\eta -1,p,D}$$

$$\le \sum _{|\mathbf{\alpha }|\le \eta -1}\frac{1}{\mathbf{\alpha }!}\sum _{i=1}^d\Vert {\partial }_{\mathbf{x}}^{\mathbf{\alpha }}{\partial }_{x_i}{v}_i{\Vert }_{p,D}\le \eta \sum _{|\mathbf{\alpha }|\le \eta }\frac{1}{\mathbf{\alpha }!}\sum _{i=1}^d\Vert {\partial }_{\mathbf{x}}^{\mathbf{\alpha }}{v}_i{\Vert }_{p,D}$$

$$\le \eta d\sum _{|\mathbf{\alpha }|\le \eta }\frac{1}{\mathbf{\alpha }!}\Vert {\partial }_{\mathbf{x}}^{\mathbf{\alpha }}\mathbf{v}{\Vert }_{p,D}\le \eta d\Vert \mathbf{v}{\Vert }_{\eta ,p,D}.$$

Full text

Funding: The work of the authors was supported by the Swiss National Science Foundation (SNSF) through the project “Multilevel Methods and Uncertainty Quantification in Cardiac Electrophysiology” (grant 205321_169599).

1. Introduction

The numerical approximation of quantities of interest, such as expectation, variance, or more general output functionals, of the solution of a diffusion problem with a scalar random diffusion coefficient with multilevel collocation or multilevel quadrature methods has been considered previously, see e.g. [2, 8, 13, 14, 20, 23, 27, 34] and the references therein; in this isotropic case, the mixed smoothness required for the use of such multilevel methods has been provided in [9] for uniformly elliptic diffusion coefficients and in [26] for log-normally distributed diffusion coefficients. However, in simulations of certain diffusion phenomena in science and engineering, the diffusion that needs to be modeled may not necessarily be isotropic. One specific application we have in mind here stems from cardiac electrophysiology, where the electrical activation of the human heart is considered. It is known that the fibrous structure of the heart plays a major role when considering the electrical and mechanical properties of the heart. And while the fibres have a complex and generally well-organised structure, see e.g. [11, 30, 31, 32], the exact fibre orientation may vary between individuals and also over time in an individual, for example due to the presence of scaring of the heart.

More generally, we wish to be able to model diffusion in a fibrous media, where fibre direction and diffusion strength in fibre direction are subject to uncertainty. For this setting, the following random anisotropic diffusion coefficient was defined in [24]:

[TABLE]

where $a$ is a given value and $\mathbf{V}$ is a random vector-valued field, over a given spatial domain $D$ and a given probability space $(\Omega,\mathcal{F},\mathbb{P})$. The fibre direction is hence given by $\mathbf{V}/\lVert\mathbf{V}\rVert_{2}$ with the diffusion strength in the fibre direction being $\lVert\mathbf{V}\rVert_{2}$ and the diffusion strength perpendicular to the fibre direction is defined by $a$. While we only consider this model hereafter, the techniques we use may also be applied straightforwardly to other models of random anisotropic diffusion coefficients, such as, for example, the following model in three spatial dimensions:

[TABLE]

Here, $\mathbf{f}$, $\mathbf{t}$ and $\mathbf{s}$ are vector fields describing the fibre, the transverse sheet and the sheet normal directions in the heart, which at each point in $D$ yields an orthonormal basis. These vector fields could, for example, be derived from measurements or be generated by an algorithm such as the Laplace-Dirichlet Rule-Based algorithm described in [3]. Note that, in this model, the diffusion strengths are random fields $a_{\mathbf{f}}$, $a_{\mathbf{t}}$ and $a_{\mathbf{s}}$, and the fibre and sheet directions are locally angularily perturbed by random fields $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$.

We shall consider the second order diffusion problem with this uncertain diffusion coefficient $\mathbf{A}$ given by

[TABLE]

with the known function $f$ as a source. The result of this article is then as follows. Having spatial $H^{s}$-regularity of the underlying diffusion problem, given by sufficient smoothness of the right hand $f$ side and the domain $D$, then the random solution $u$ admits analytic regularity with respect to the stochastic parameter also in the $H^{s}(D)$-norm provided that the random vector-valued field offers enough spatial regularity. This mixed regularity is the essential ingredient in order to apply multilevel collocation or multilevel quadrature methods without deteriorating the rate of convergence, see [20] for instance.

The rest of the article is organised as follows: In Section 2, we provide basic definitions and notation for the functional analytic framework to be able to state and then also reformulate the model problem, by using the Karhunen-Loève expansion of the diffusion describing random vector-valued field $\mathbf{V}$, into its stochastically parametric and spatially weak formulation. Section 3 then deals with the regularity of the solution of the stochastically parametric and spatially weak formulation of the model problem with respect to the stochastic parameter and some given higher spatial regularity in the model problem. We then use the fact that the higher spatial regularity can be kept, when considering the regularity of the solution with respect to the stochastic parameter, to arrive at convergence rates when considering multilevel quadrature, such as multilevel quasi-Monte Carlo quadrature, to approximate the solution’s mean and second moment. Numerical examples are provided in Section 4 as validation; specifically we use multilevel quasi-Monte Carlo quadrature to approximate the solution’s mean and second moment in a setting with three spatial dimensions. Lastly, we give our conclusions in Section 5.

2. Problem formulation

2.1. Notation and precursory remarks

For a given Banach space $\mathcal{X}$ and a complete measure space $\mathcal{M}$ with measure $\mu$ the space $L_{\mu}^{p}(\mathcal{M};\mathcal{X})$ for $1\leq p\leq\infty$ denotes the Bochner space, see [25], which contains all equivalence classes of strongly measurable functions $v\colon\mathcal{M}\to\mathcal{X}$ with finite norm

[TABLE]

A function $v\colon\mathcal{M}\to\mathcal{X}$ is strongly measurable if there exists a sequence of countably-valued measurable functions $v_{n}\colon\mathcal{M}\to\mathcal{X}$, such that for almost every $m\in\mathcal{M}$ we have $\lim_{n\to\infty}v_{n}(m)=v(m)$. Note that, for finite measures $\mu$, we also have the usual inclusion $L_{\mu}^{p}(\mathcal{M};\mathcal{X})\supset L_{\mu}^{q}(\mathcal{M};\mathcal{X})$ for $1\leq p<q\leq\infty$.

Let $\mathcal{X}$, $\mathcal{X}_{1},\ldots,\mathcal{X}_{r}$ and $\mathcal{Y}$ be Banach spaces, then we denote the Banach space of bounded, linear maps from $\mathcal{X}$ to $\mathcal{Y}$ as $\mathcal{B}(\mathcal{X};\mathcal{Y})$; furthermore, we recursively define

[TABLE]

and the special case

[TABLE]

For $\mathbf{T}\in\mathcal{B}(\mathcal{X}_{1},\ldots,\mathcal{X}_{r};\mathcal{Y})$ and $\mathbf{v}_{j}\in\mathcal{X}_{j}$ we use the notation $\mathbf{T}\mathbf{v}_{1}\cdots\mathbf{v}_{r}\mathrel{\mathrel{\mathop{:}}=}\mathbf{T}(\mathbf{v}_{1},\ldots,\mathbf{v}_{r})\in\mathcal{Y}$.

Subsequently, we will always equip $\mathbb{R}^{d}$ with the norm $\lVert\cdot\rVert_{2}$ induced by the canonical inner product $\langle\cdot,\cdot\rangle$ and $\mathbb{R}^{d\times d}$ with the induced norm $\lVert\cdot\rVert_{2}$. Then, for $\mathbf{v},\mathbf{w}\in\mathbb{R}^{d}$, the Cauchy-Schwartz inequality gives us

[TABLE]

where equality holds only if $\mathbf{v}=\mathbf{w}$, and we also have, by straightforward computation, that

[TABLE]

We also note that to avoid the use of generic but unspecified constants in certain formulas we use $C\lesssim D$ to mean that $C$ can be bounded by a multiple of $D$, independently of parameters which $C$ and $D$ may depend on. Obviously, $C\gtrsim D$ is defined as $D\lesssim C$ and we write $C\eqsim D$ if $C\lesssim D$ and $C\gtrsim D$. Lastly, note that for the natural numbers $\mathbb{N}$ denotes them including [math] and $\mathbb{N}^{*}$ excluding [math].

2.2. Model problem

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a separable, complete probability space. Then, we consider the following second order diffusion problem with a random anisotropic diffusion coefficient

[TABLE]

where $D\subset\mathbb{R}^{d}$ is a Lipschitz domain with $d\geq 1$ and the function $f\in H^{-1}(D)$ describes the known source. The random anisotropic diffusion coefficient is given as the random matrix field \mathbf{A}\in L_{\mathbb{P}}^{\infty}\big{(}\Omega;L^{\infty}(D;\mathbb{R}^{d\times d})\big{)}, which satisfies the uniform ellipticity condition

[TABLE]

for some constants $0<\underline{a}\leq\overline{a}<\infty$ and is almost surely symmetric almost everywhere. Without loss of generality, we assume $\underline{a}\leq 1\leq\overline{a}$.

We specifically consider diffusion coefficients that are of form

[TABLE]

where $a\in\mathbb{R}$ is a given positive number and \mathbf{V}\in L_{\mathbb{P}}^{\infty}\big{(}\Omega;L^{\infty}(D;\mathbb{R}^{d})\big{)} is a random vector-valued field. We note that such a field $\mathbf{A}$ accounts for a medium that has homogeneous diffusion strength $a$ perpendicular to $\mathbf{V}$ and has diffusion strength \big{\lVert}\mathbf{V}[\omega](\mathbf{x})\big{\rVert}_{2} in the direction of $\mathbf{V}$. The randomness of the specific direction and length of $\mathbf{V}$ therefore quantifies uncertainty of this notable direction and its diffusion strength. To guarantee the uniform ellipticity condition (2), we require that

[TABLE]

as well as $\underline{a}\leq a\leq\overline{a}$.

It is assumed that the spatial variable $\mathbf{x}$ and the stochastic parameter $\omega$ of the random field have been separated by the Karhunen-Loève expansion of $\mathbf{V}$, yielding a parametrised expansion

[TABLE]

where $\mathbf{y}=(y_{k})_{k\in\mathbb{N}^{*}}\in\square\mathrel{\mathrel{\mathop{:}}=}[{-1},1]^{\mathbb{N}^{*}}$ is a sequence of uncorrelated random variables, see e.g. [24]. In the following, we will denote the pushforward of the measure $\mathbb{P}$ onto $\square$ as $\mathbb{P}_{\mathbf{y}}$. Then, we also view $\mathbf{A}[\mathbf{y}](\mathbf{x})$ and $u[\mathbf{y}](\mathbf{x})$ as being parametrised by $\mathbf{y}$ and restate (1) as

[TABLE]

We now impose some common assumptions, which make the Karhunen-Loève expansion computationally feasible.

Assumption 2.1**.**

The random variables $(y_{k})_{k\in\mathbb{N}^{*}}$ are independent and identically distributed. Moreover, they are uniformly distributed on \big{[}{-1},1\big{]}.

Lastly, we note that the spatially weak form of (6) is given by

[TABLE]

This also entails the well known stability estimate.

Lemma 2.2**.**

There is a unique solution u\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;H_{0}^{1}(D)\big{)} of (7), which fulfils

[TABLE]

where $c_{V}$ is the Poincaré-Friedrichs constant of $H_{0}^{1}(D)$.

3. Parametric regularity and multilevel quadrature

We now derive regularity estimates for the solution $u$ of (7) and apply multilevel quadrature for approximating the mean of $u$. The regularity estimates are based on the following assumption on the decay of the expansion of $\mathbf{V}$.

Assumption 3.1**.**

We assume that the ${\bm{\psi}}_{k}$ are elements of $W^{\kappa,\infty}(D;\mathbb{R}^{d})$ for a $\kappa\in\mathbb{N}$ and that the sequence ${\bm{\gamma}}_{\kappa}=(\gamma_{\kappa,k})_{k\in\mathbb{N}}$, given by

[TABLE]

is at least in $\ell^{1}(\mathbb{N})$, where we have defined ${\bm{\psi}}_{0}\mathrel{\mathrel{\mathop{:}}=}\operatorname{\mathbb{E}}[\mathbf{V}]$ and $\sigma_{0}\mathrel{\mathrel{\mathop{:}}=}1$. Furthermore, we define

[TABLE]

We furthermore assume that the vector field $\mathbf{V}$ is given by a finite rank Karhunen-Loève expansion, i.e.

[TABLE]

where $\square\mathrel{\mathrel{\mathop{:}}=}[{-1},1]^{M}$. We note that the regularity estimates however will not depend on the rank $M$ and therefore, if necessary, a finite rank can be attained by appropriate truncation.

Furthermore, for the regularity estimates we also require an elliptic regularity result.

Assumption 3.2**.**

Let $\mathcal{R}_{\kappa}$ be a Banach space with norm $\lVert\cdot\rVert_{\mathcal{R}_{\kappa}}$ such that, for all

[TABLE]

that fulfil (2), we have that the problem of solving

[TABLE]

for any $h\in H^{\kappa-1}(D)$ has a unique solution $u\in H_{0}^{1}(D)$, which also lies in $H^{\kappa+1}(D)$, with

[TABLE]

where $C_{\kappa,er}$ only depends on $D$ and continuously on $\lVert\mathbf{A}\rVert_{\mathcal{R}_{\kappa}}$.

We assume from here on that $\mathbf{A}$ also lies in $L_{\mathbb{P}_{\mathbf{y}}}^{\infty}(\square;\mathcal{R}_{\kappa})$.

Note, that for $\kappa=0$, this reduces to the stability estimate, for which the parametric regularity may be found in [24]. Therefore, we will hereafter only consider the case where $\kappa\geq 1$. Such an elliptic regularity estimate for example is known for $\kappa=1$, when the domain is convex and bounded and $\mathcal{R}_{\kappa}=C^{0,1}(\overline{D};\mathbb{R}^{d\times d})$, see [16, Theorems 3.2.1.2 and 3.1.3.1]. The elliptic regularity estimate is also known to hold for $\kappa\geq 1$ and $d=2$, when the domain’s boundary is smooth and $\mathcal{R}_{\kappa}=W^{\kappa,\infty}(D;\mathbb{R}^{d\times d})$, see [7].

The rest of this section is now split into three subsections. The first subsection is dedicated to introducing some useful norms as well as deriving Corollary 3.11 and Theorem 3.12. These two results are then used in the second subsection to step-by-step provide regularity estimates for the different terms that make up the diffusion coefficient, based on Assumption 3.1 on the decay of the expansion of $\mathbf{V}$, which then yields regularity estimates for the diffusion coefficient. By using this decay as well as Assumption 3.2, we derive the regularity estimates of the solution $u$ in Theorem 3.18. In the third subsection, we then briefly discuss what kind of convergence rates and computational complexity this regularity of $u$ yields, for the example of approximating $\operatorname{\mathbb{E}}[u]$ using multilevel quadrature methods.

3.1. Precursory remarks

For the Sobolev–Bochner spaces $W^{\eta,p}(D;\mathcal{X})$ with $\eta\in\mathbb{N}$ and $1\leq p\leq\infty$, we introduce the norms given by

[TABLE]

for $v\in W^{\eta,p}(D;\mathcal{X})$ where $\mathcal{X}$ is a Banach space with norm $\lVert\cdot\rVert_{\mathcal{X}}$ and where we make use of the shorthand

[TABLE]

We may omit the specification of the Banach space $\mathcal{X}$, for example when $\mathcal{X}$ is the space $\mathbb{R}$, $\mathbb{R}^{d}$ or $\mathbb{R}^{d\times d}$.

For these norms, we have the following lemma giving a bound after applying a $\operatorname{div}_{\mathbf{x}}$, $\operatorname{D}_{\mathbf{x}}$ or $\operatorname{\nabla}_{\mathbf{x}}$:

Lemma 3.3**.**

Let $\eta\in\mathbb{N}^{*}$, $1\leq p\leq\infty$. For $\mathbf{v}\in W^{\eta,p}(D;\mathbb{R}^{d})$ we have that $\operatorname{div}_{\mathbf{x}}\mathbf{v}\in W^{\eta-1,p}(D;\mathbb{R})$ with

[TABLE]

and $\operatorname{D}_{\mathbf{x}}\mathbf{v}\in W^{\eta-1,p}(D;\mathbb{R}^{d\times d})$ with

[TABLE]

For $v\in W^{\eta,p}(D)$ we have that $\operatorname{\nabla}_{\mathbf{x}}v\in W^{\eta-1,p}(D;\mathbb{R}^{d})$ with

[TABLE]

Proof.

We calculate

[TABLE]

We also may compute

[TABLE]

as well as

[TABLE]

The Leibniz rule also yields the following kind of submultiplicativity for these norms:

Lemma 3.4**.**

Let $\eta\in\mathbb{N}$, $1\leq p,p^{\prime}\leq\infty$, $\mathcal{X}$ and $\mathcal{Y}$ be Banach spaces and

[TABLE]

with $q=(p^{-1}+p^{\prime-1})^{-1}\geq 1$. Then, we have

[TABLE]

Proof.

Let ${\bm{\alpha}},{\bm{\beta}}\in\mathbb{N}^{d}$ be two multi-indices, then we have

[TABLE]

We now can calculate

[TABLE]

By a change of variables, i.e. replacing ${\bm{\alpha}}$ with ${\bm{\alpha}}+{\bm{\beta}}$, and remarking that

[TABLE]

we find the identity

[TABLE]

Consequently, we arrive at the desired estimate:

[TABLE]

By induction, we thus arrive at the following corollary:

Corollary 3.5**.**

Let $\eta\in\mathbb{N}$, $1\leq p,p_{1},\ldots,p_{r}\leq\infty$, $\mathcal{X}_{1},\ldots,\mathcal{X}_{r}$ and $\mathcal{Y}$ be Banach spaces and

[TABLE]

with $q=(p^{-1}+p_{1}^{-1}+\cdots+p_{r}^{-1})^{-1}\geq 1$. Then, we have

[TABLE]

As we will need the Faà di Bruno formula, see [10], we just restate it here — slightly adapted to our notation and usage — for reference:

Remark 3.6**.**

Given $\mathbf{W}\colon X\to\mathcal{Y}$ and $\mathbf{v}\colon D\to\mathcal{X}$, where $\mathcal{X}$ and $\mathcal{Y}$ are Banach spaces, $X\subset\mathcal{X}$ is open with $\operatorname{img}_{D}\mathbf{v}\subset X$ and $\mathbf{W},\mathbf{v}$ are both sufficiently differentiable for the formula to make sense, then

[TABLE]

for ${\bm{\alpha}}\in\mathbb{N}^{M}$ with ${\bm{\alpha}}\neq{\bm{0}}$, where $P({\bm{\alpha}},r)$ is the set of integer partitions of a multiindex ${\bm{\alpha}}$ into $r$ non-vanishing multiindices, given by

[TABLE]

It also follows from [10] that, for ${\bm{\alpha}}\in\mathbb{N}^{d}$ and $r\in\mathbb{N}$, we have

[TABLE]

where $S_{n,r}$ denotes the Stirling numbers of the second kind, see [1].

Lastly, as we know that $\sum_{r=1}^{\lvert{\bm{\alpha}}\rvert}r!S_{\lvert{\bm{\alpha}}\rvert,r}$ equals the $\lvert{\bm{\alpha}}\rvert$-th ordered Bell number, we can bound it, see [4], giving

[TABLE]

Lastly, the Faà di Bruno formula now yields the following lemma:

Lemma 3.7**.**

Let $\eta\in\mathbb{N}$, $1\leq p\leq\infty$, $\mathcal{X}$ and $\mathcal{Y}$ be Banach spaces,

[TABLE]

$X\subset\mathcal{X}$* be open with $\operatorname{img}_{D}\mathbf{v}\subset X$, $\mathbf{W}\colon X\to\mathcal{Y}$ $\eta$-times Fréchet differentiable and, for $0\leq r\leq\eta$,*

[TABLE]

Then, we have

[TABLE]

with

[TABLE]

Proof.

The Faà di Bruno formula leads to

[TABLE]

Now, with Corollary 3.5,

[TABLE]

leads to

[TABLE]

Furthermore, we also introduce the shorthand notations

[TABLE]

for v\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;W^{\eta,p}(D;\mathcal{X})\big{)} where $\mathcal{X}$ is a Banach space with norm $\lVert\cdot\rVert_{\mathcal{X}}$. We may also omit the specification of the Banach space $\mathcal{X}$ in this shorthand, for example when $\mathcal{X}$ is the space $\mathbb{R}$, $\mathbb{R}^{d}$ or $\mathbb{R}^{d\times d}$.

With this notation we still carry over the previous results, yielding the following lemmata and corollaries. Lemma 3.3 directly transforms to:

Corollary 3.8**.**

Let $\eta\in\mathbb{N}^{*}$, $1\leq p\leq\infty$. For \mathbf{v}\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;W^{\eta,p}(D;\mathbb{R}^{d})\big{)} we have that \operatorname{div}_{\mathbf{x}}\mathbf{v}\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;W^{\eta-1,p}(D;\mathbb{R})\big{)} with

[TABLE]

and \operatorname{D}_{\mathbf{x}}\mathbf{v}\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;W^{\eta-1,p}(D;\mathbb{R}^{d\times d})\big{)} with

[TABLE]

For v\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;W^{\eta,p}(D)\big{)} we have that \operatorname{\nabla}_{\mathbf{x}}v\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;W^{\eta-1,p}(D;\mathbb{R}^{d})\big{)} with

[TABLE]

We use Lemma 3.4 to arrive at the following result:

Lemma 3.9**.**

Let $\eta,\in\mathbb{N}$, $1\leq p,p^{\prime}\leq\infty$, ${\bm{\nu}}\in\mathbb{N}^{M}$, $\mathcal{X}$ and $\mathcal{Y}$ be Banach spaces and

[TABLE]

for all ${\bm{\alpha}}\leq{\bm{\nu}}$ with $q=(p^{-1}+p^{\prime-1})^{-1}\geq 1$. Then, we have

[TABLE]

with

[TABLE]

for all ${\bm{\alpha}}\leq{\bm{\nu}}$.

Proof.

The Leibniz rule and the application of Lemma 3.4 yield:

[TABLE]

Again, similar as before by induction, we arrive at:

Corollary 3.10**.**

Let $\eta,\in\mathbb{N}$, $1\leq p,p_{1},\ldots,p_{r}\leq\infty$, ${\bm{\nu}}\in\mathbb{N}^{M}$, $\mathcal{X}_{1},\ldots,\mathcal{X}_{r}$ and $\mathcal{Y}$ be Banach spaces and

[TABLE]

for all ${\bm{\alpha}}\leq{\bm{\nu}}$ with $q=(p^{-1}+p_{1}^{-1}+\cdots+p_{r}^{-1})^{-1}\geq 1$. Then, we have

[TABLE]

with

[TABLE]

for all ${\bm{\alpha}}\leq{\bm{\nu}}$.

Moreover, we also will require the following corollary:

Corollary 3.11**.**

Let $\eta,\in\mathbb{N}$, $1\leq p_{1},p_{2}\leq\infty$, ${\bm{\nu}}\in\mathbb{N}^{M}$, $\mathcal{X}_{1},\mathcal{X}_{2}$ and $\mathcal{Y}$ be Banach spaces and

[TABLE]

for all ${\bm{\alpha}}\leq{\bm{\nu}}$ with $q=(p_{1}^{-1}+p_{2}^{-1})^{-1}\geq 1$. Then, we have

[TABLE]

with

[TABLE]

for all ${\bm{\alpha}}\leq{\bm{\nu}}$.

Lastly, from Theorem 3.12 we can the derive the following:

Theorem 3.12**.**

Let $\eta\in\mathbb{N}$, $1\leq p\leq\infty$, ${\bm{\nu}}\in\mathbb{N}^{M}$, $\mathcal{X}$ and $\mathcal{Y}$ be Banach spaces,

[TABLE]

$X\subset\mathcal{X}$* be open with $\operatorname{img}_{\square}\operatorname{img}_{D}\mathbf{v}\subset X$, $\mathbf{W}\colon X\to\mathcal{Y}$ be $\eta+\lvert{\bm{\nu}}\rvert$-times Fréchet differentiable and, for ${\bm{\alpha}}\leq{\bm{\nu}}$ and $0\leq t\leq\eta+\lvert{\bm{\nu}}\rvert$,*

[TABLE]

Then, we have

[TABLE]

with

[TABLE]

and, for ${\bm{\alpha}}\neq{\bm{0}}$,

[TABLE]

Proof.

The application of Lemma 3.7 leads to

[TABLE]

from which, with $s=0$, the assertion follows directly for ${\bm{\alpha}}={\bm{0}}$.

For ${\bm{\alpha}}\neq{\bm{0}}$, we remark that the Faà di Bruno formula and Corollary 3.10 yield

[TABLE]

3.2. Parametric regularity

We first consider the two terms

[TABLE]

which appear in the diffusion coefficient. For these we have:

Lemma 3.13**.**

We have for all ${\bm{\alpha}}\in\mathbb{N}^{M}$ that

[TABLE]

with $k_{\mathbf{B}}\mathrel{\mathrel{\mathop{:}}=}k_{C}\mathrel{\mathrel{\mathop{:}}=}c_{{\bm{\gamma}}_{\kappa}}^{2}$. Furthermore, we know that $\operatorname{img}_{\square}\operatorname{img}_{D}C\subset{[\underline{a}^{2},\overline{a}^{2}]}$.

Proof.

More verbosely, $\mathbf{B}$ is given by

[TABLE]

from which we can derive the first order derivatives, yielding

[TABLE]

and from those also the second order derivatives. They are given by

[TABLE]

Since the second order derivatives with respect to $\mathbf{y}$ are constant, all higher order derivatives with respect to $\mathbf{y}$ vanish.

We obviously have

[TABLE]

From (8) we can now derive the bound

[TABLE]

and (9) leads us to

[TABLE]

Therefore, we have

[TABLE]

and are finished since $c_{{\bm{\gamma}}_{\kappa}}\geq 2$.

Starting from the equation

[TABLE]

the bound for $C$ clearly follows analogously as for $\mathbf{B}$. Lastly, we note that (4) directly implies that $\operatorname{img}_{\square}\operatorname{img}_{D}C\subset{[\underline{a}^{2},\overline{a}^{2}]}$. ∎

Next, we consider the terms

[TABLE]

for which we have:

Lemma 3.14**.**

We know for all ${\bm{\alpha}}\in\mathbb{N}^{M}$ that

[TABLE]

with

[TABLE]

Proof.

Because $D=v\circ C$ with $v(x)=x^{-1}$ and $E=w\circ C$ with $w(x)=\sqrt{x}$ are composite functions, we can employ Theorem 3.12 to bound their derivatives. For this, we remark that the $t$-th derivative of $v$ and $w$ are given by

[TABLE]

and

[TABLE]

where c_{t}\mathrel{\mathrel{\mathop{:}}=}\prod_{i=0}^{t-1}\big{(}\frac{1}{2}-i\big{)}. For $x\in{[\underline{a}^{2},\overline{a}^{2}]}$ we therefore have

[TABLE]

which implies

[TABLE]

using Lemma 3.13 — with which we then furthermore arrive at the estimates

[TABLE]

and, analogously,

[TABLE]

as well as, for ${\bm{\alpha}}\neq{\bm{0}}$,

[TABLE]

and, again analogously,

[TABLE]

Combining these estimates finally yields the assertions. ∎

We can now consider the two terms

[TABLE]

For these we have:

Lemma 3.15**.**

We have for all ${\bm{\alpha}}\in\mathbb{N}^{M}$ that

[TABLE]

with

[TABLE]

Proof.

We use Corollary 3.11 with the bounds from Lemma 3.13 and Lemma 3.14 to arrive at

[TABLE]

and

[TABLE]

Lastly, the combinatorial identity

[TABLE]

yields the bounds

[TABLE]

and

[TABLE]

with which the assertion follows. ∎

We now can consider the complete diffusion coefficient, yielding:

Theorem 3.16**.**

The derivatives of the diffusion matrix $\mathbf{A}$, defined in (3), satisfy

[TABLE]

for all ${\bm{\alpha}}\in\mathbb{N}^{M}$ with

[TABLE]

Proof.

We can state $\mathbf{A}$ as

[TABLE]

which, by taking the norm of the derivative and using Lemma 3.15, yields

[TABLE]

We now define the modified sequence ${\bm{\mu}}_{\kappa}=(\mu_{\kappa,k})_{k\in\mathbb{N}}$ as

[TABLE]

and also

[TABLE]

thus, we have

[TABLE]

Now, Assumption 3.2 directly implies the following result:

Lemma 3.17**.**

The unique solution u\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;H_{0}^{1}(D)\big{)} of (7) fulfils u\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;H^{\kappa+1}(D)\big{)}, with

[TABLE]

where

[TABLE]

However, by also leveraging the higher spatial regularity in the Karhunen-Loève expansion of the random vector-valued field, we can show that the solution $u$ admits analytic regularity with respect to the stochastic parameter $\mathbf{y}$ also in the $H^{\kappa+1}(D)$-norm. This mixed regularity is then the essential ingredient when applying multilevel methods.

Theorem 3.18**.**

The derivatives in $\mathbf{y}$ of the solution $u$ of (7) satisfy

[TABLE]

Proof.

By differentiation of the variational formulation (7) with respect to $\mathbf{y}$ we arrive, for arbitrary $v\in H_{0}^{1}(D)$, at

[TABLE]

Applying the Leibniz rule on the left-hand side yields

[TABLE]

Then, by rearranging and using the linearity of the gradient, we find

[TABLE]

Using Green’s identity, we can then write

[TABLE]

Thus, we arrive at

[TABLE]

from which we derive

[TABLE]

where

[TABLE]

We note that, by definition of $c$, we have $c\geq 2$ and furthermore, because of Lemma 3.17, we also have that \big{\ltvert}u\big{\rtvert}_{H^{1}(D)}\leq c, which means that the assertion is true for $\lvert{\bm{\alpha}}\rvert=0$. Thus, we can use an induction over $\lvert{\bm{\alpha}}\rvert$ to prove the hypothesis

[TABLE]

for $\lvert{\bm{\alpha}}\rvert>0$.

Let the assertions hold for all ${\bm{\alpha}}$, which satisfy $\lvert{\bm{\alpha}}\rvert\leq n-1$ for some $n\geq 1$. Then, we know for all ${\bm{\alpha}}$ with $\lvert{\bm{\alpha}}\rvert=n$ that

[TABLE]

Making use of the combinatorial identity (10) yields

[TABLE]

Now, since $c\geq 2$, we have $c\leq 2(c-1)$ and hence also

[TABLE]

This completes the proof. ∎

3.3. Numerical quadrature in the parameter

Coming from the solution of (7), that is u\in L_{\mathbb{P}_{\mathbf{y}}}^{\infty}\big{(}\square;H^{\kappa+1}(D)\big{)}, we now wish to know the moments of $u$. In this section, we will therefore consider the approximation of the mean of $u$.

The mean of $u$ is given by the Bochner integral

[TABLE]

Therefore, we may proceed to approximate it by considering a generic quadrature method $Q_{N}$; that is

[TABLE]

where

[TABLE]

are the weight and evaluation point pairs. We assume that the quadrature chosen fulfils

[TABLE]

for some constants $C>0$ and $r>0$.

We will employ the quasi-Monte Carlo quadrature based on the Halton sequence, i.e. $\omega_{i}^{(N)}=1/N$ and ${\bm{\xi}}_{i}^{(N)}=2{\bm{\chi}}_{i}-{\bm{1}}$, where ${\bm{\chi}}_{i}$ denotes the $i$-th $M$-dimensional Halton point, cf. [18]. Then, we know that, given that there exists an $\varepsilon>0$ such that $\gamma_{\kappa,k}\leq ck^{-3-\varepsilon}$ holds for some $c>0$, for every $\delta>0$ there exists a constant $C=C(\delta)>0$ such that (12) holds for $r=1-\delta$, see e.g. [22] which is a consequence of [35]. Clearly, other, possibly more sophisticated, quadrature methods may also be considered, for example, other quasi-Monte Carlo quadratures, such as those based on the Sobol sequence or other low-discrepancy sequences as well as their higher-order adaptations, and anisotropic sparse grid quadratures, see e.g. [12, 17, 29, 33].

To approximate the mean of $u$ as in (11), we require the values $u[\mathbf{y}]$ for $\mathbf{y}={\bm{\xi}}_{i}$. These values can be approximated by $u_{l}[\mathbf{y}]$, where $u_{l}$ is the Galerkin approximation of the spatially weak formulation on a finite dimensional subspace $V_{l}$ of $H_{0}^{1}(D)$; that is, $u_{l}$ is the solution of

[TABLE]

We assume that a sequence of $V_{l}$ can be chosen for $l\in\mathbb{N}$ such that there is a constant $K$ with

[TABLE]

For example, we can consider $V_{l}$ to be the spaces of continuous finite elements of order $\kappa$ coming from a sequence of quasi-uniform meshes $\mathcal{T}_{l}$ using isoparametric elements, where the mesh size behaves like $2^{-l}$. Then, it is known from finite element theory that we have (13) with $K\sim\lVert u\rVert_{L_{\mathbb{P}_{\mathbf{y}}}^{\infty}(\square;H^{\kappa+1}(D))}$, see e.g. [5, 6].

The combination of the error estimates (12) and (13) then leads to

[TABLE]

Thus, choosing N_{l}\mathrel{\mathrel{\mathop{:}}=}\big{\lceil}2^{\kappa l/r}\big{\rceil} finally yields

[TABLE]

In contrast, the mixed regularity, shown before in Theorem 3.18, allows us to consider a multilevel adaptation, which may be given as

[TABLE]

where

[TABLE]

Indeed, this is the sparse grid combination technique as introduced in [15], see also [14, 20]. It thus follows that

[TABLE]

For complexity considerations, we shall consider a quadrature that is nested, i.e. we may set ${\bm{\xi}}_{i}={\bm{\xi}}_{i}^{(N)}$ as it does not depend on $N$. Then, we note that $Q_{l}^{\mathrm{ML}}[u_{0},\ldots,u_{l}]$ may explicitly be stated as

[TABLE]

Computing $Q_{N_{l}}[u_{l}]$ requires thus the values $u_{l,i}(\mathbf{x})\mathrel{\mathrel{\mathop{:}}=}u(\mathbf{x},{\bm{\xi}}_{i})$, which can be derived by solving

[TABLE]

Generally, when considering a sequence of finite element spaces $V_{l}$ as described above, the number of degrees of freedom behaves like \mathcal{O}\big{(}2^{ld}\big{)} and computing one $u_{l,i}$ using state of the art methods will have a complexity that is \mathcal{O}\big{(}2^{ld}\big{)}. As this has to be done $N_{l}$ times for the calculation of $Q_{N_{l}}[u_{l}]$, a complexity scaling is obtained that is \mathcal{O}\big{(}2^{l(\kappa/r+d)}\big{)}. Therefore, for the computation of the multilevel quadrature $Q_{l}^{\mathrm{ML}}[u_{0},\ldots,u_{l}]$, we arrive at an over-all complexity of

[TABLE]

We mention that also non-nested quadrature formulae can be used but lead to a somewhat larger constant in the complexity estimate, see [14] for the details.

Remark 3.19**.**

If we redefine the $N_{l}$ as N_{l}\mathrel{\mathrel{\mathop{:}}=}\big{\lceil}l^{(1+\varepsilon)/r}2^{\kappa l/r}\big{\rceil} for an $\varepsilon>0$, then we have

[TABLE]

and, as proposed in [2], we arrive at

[TABLE]

So, the logarithmic factor, which shows up in the convergence rate, can be removed by increasing the quadrature accuracy slightly faster. Note that this modification increases the hidden constant with a dependance on $\varepsilon$.

Remark 3.20**.**

In the particular situation of a standard quasi-Monte Carlo method, we can consider $\delta^{\prime}$ such that $\delta>\delta^{\prime}>0$. Then, the quadrature error satisfies the estimate

[TABLE]

With a similar argument as in [2], it follows that

[TABLE]

That is, the logarithmic factor, which shows up in the convergence rate, is removed at the cost of replacing the constant $C_{\delta}$ with $C_{\delta^{\prime}}$ and adds a constant with a dependance on $\delta^{\prime}$ yielding an increased hidden constant.

While we have exclusively considered the case of the mean of the solution $u$ here, we do note that analogous statements may also be shown for example for the higher-order moments, see [20] for instance.

4. Numerical Results

We will now consider two examples of the model problem (1) with a diffusion coefficient of form (3) using the unit cube $D\mathrel{\mathrel{\mathop{:}}=}(0,1)^{3}$ as the domain of computations. Therefore, in view of $H^{2}$-regularity of the spatial problem under consideration, we are only considering the situation with $\kappa=1$. In both examples, we set the global strength $a$ to $a\mathrel{\mathrel{\mathop{:}}=}0.12$ and choose the right hand side $f\equiv 1$. For convenience, we define

[TABLE]

Example 1**.**

In this first example, we choose the description of $\mathbf{V}$ to be defined by

[TABLE]

and

[TABLE]

We note that for $j\in\{2,3\}$ the covariance in the normal direction on the parts of the boundary with $x_{j}\in\{0,1\}$ is suppressed.

Example 2**.**

For this second example we choose the description of $\mathbf{V}$ to be defined by

[TABLE]

and

[TABLE]

Here, the covariance in the normal direction on all of the boundary is suppressed.

The numerical implementation is performed with aid of the problem-solving environment DOLFIN [28], which is a part of the FEniCS Project [28]. The Karhunen-Loève expansion of the vector field $\mathbf{V}$ is computed by the pivoted Cholesky decomposition, see [19, 21] for the details. For the finite element discretisation, we employ the sequence of nested triangulations $\mathcal{T}_{l}$, yielded by successive uniform refinement, i.e. cutting each tetrahedron into $8$ tetrahedra. The base triangulation $\mathcal{T}_{0}$ consists of $6\cdot 2^{3}=48$ tetrahedra. Then, we use interpolation with continuous element-wise linear functions and the truncated pivoted Cholesky decomposition for the Karhunen-Loève expansion approximation and continuous element-wise linear functions in space. The truncation criterion for the pivoted Cholesky decomposition is that the relative trace error is smaller than $10^{-4}\cdot 4^{-l}$.

Since the exact solutions of the examples are unknown, the errors will have to be estimated. Therefore, in this section, we will estimate the errors for the levels [math] to $5$ by substituting the exact solution with the approximate solution computed on the level $6$ triangulation $\mathcal{T}_{6}$ using the quasi-Monte Carlo quadrature based on Halton points with $10^{4}$ samples.

For every level, we also define the number of samples used by the quasi-Monte Carlo method based on Halton points (QMC); we choose

[TABLE]

with $\delta\mathrel{\mathrel{\mathop{:}}=}0.2$; see Table 1 for the resulting values of $N_{l}$. This then also implies the amount of samples used on the different levels when using the multilevel quasi-Monte Carlo method based on Halton points (MLQMC). Based on these choices, we expect to see an asymptotic rate of convergence of $2^{-l}$ in the $H^{1}$-norm for the mean and in the $W^{1,1}$-norm for the variance.

Figures 1 and 2 show the estimated errors of the solution’s first moment on the left hand side and of the solution’s second moment on the right hand side, each versus the discretisation level for the QMC and MLQMC quadrature for the two different examples. As expected, the QMC quadrature methods achieves the predicted rate of convergence in both examples, and this rate of convergence also carries over to its multilevel adaptation (MLQMC).

5. Conclusion

In this article, we have considered the second order diffusion problem

[TABLE]

with the uncertain diffusion coefficent given by

[TABLE]

This models anisotropic diffusion, where the diffusion strength in the direction given by $\mathbf{V}/\lVert\mathbf{V}\rVert_{2}$ is $\lVert\mathbf{V}\rVert_{2}$ and perpendicular to it is $a$, which can be used to model both diffusion in media that consist of thin fibres or thin sheets.

After having restated the problem in a parametric form by considering the Karhunen-Loève expansion of the random vector field $\mathbf{V}$, we have shown that, given regularity of the elliptic diffusion problem, the decay of the Karhunen-Loève expansion of $\mathbf{V}$ entirely determines the regularity of the solution’s dependence on the random parameter, also when considering this higher regularity in the spatial domain.

We then leverage this result to reduce the complexity of the approximation of the solution’s mean, by using the multilevel quasi-Monte Carlo method instead of the quasi-Monte Carlo method, while still retaining the same error rate. Indeed, while the QMC method yields a scheme, where the uncertainty added increases the complexity, this is not the case, when considering two or more spatial dimensions and the MLQMC method. That is, given elliptic regularity and up to a constant in the complexity, adding uncertainty comes for free. The numerical experiments corroborate these theoretical findings.

While we considered the use of QMC and its multilevel adaptation, one can clearly also consider other quadrature methods, such as the anisotropic sparse grid quadrature, and then reduce the complexity by passing to their multilevel adaptations. Likewise, multilevel collocation is also applicable.