A fully adaptive multilevel stochastic collocation strategy for solving elliptic PDEs with random data

TL;DR

This paper introduces a fully adaptive multilevel stochastic collocation method for elliptic PDEs with random data, combining adaptive spatial mesh refinement with anisotropic sparse grid stochastic discretization to optimize computational efficiency.

Contribution

It presents a novel hierarchical adaptive strategy that optimizes both spatial and stochastic discretizations based on sample-dependent hierarchies, with rigorous analysis and optimal error control.

Findings

Significant reduction in computational complexity compared to single-level methods.

Demonstrated reliability of error control through numerical examples.

Effective combination of adaptive spatial and stochastic discretizations.

Abstract

We propose and analyse a fully adaptive strategy for solving elliptic PDEs with random data in this work. A hierarchical sequence of adaptive mesh refinements for the spatial approximation is combined with adaptive anisotropic sparse Smolyak grids in the stochastic space in such a way as to minimize the computational cost. The novel aspect of our strategy is that the hierarchy of spatial approximations is sample dependent so that the computational effort at each collocation point can be optimised individually. We outline a rigorous analysis for the convergence and computational complexity of the adaptive multilevel algorithm and we provide optimal choices for error tolerances at each level. Two numerical examples demonstrate the reliability of the error control and the significant decrease in the complexity that arises when compared to single level algorithms and multilevel algorithms that employ adaptivity solely in the spatial discretisation or in the collocation procedure.

Figures

Tables

Table 1. Table 1: Algorithm to approximate solution functionals ψ ​ ( u ) 𝜓 𝑢 \psi(u) by an adaptive multilevel stochastic collocation method.

Algorithm: Adaptive Multilevel Stochastic Collocation Method
1.	Given $ϵ$ and $q$,  estimate $C_X$, $C_Y$, $s^{\ast }$, ${\mu }^{\ast }$ and $K=K\left(ϵ\right)$.
2.	Set coarsest spatial tolerance:   ${\eta }_{X_0}:=ϵ/\left(2C_X{q}^K\right)$.
3.	Set other spatial tolerances:  ${\eta }_{X_k}:=q^k{\eta }_{X_0}$, $k=1,\dots ,K$.
4.	Compute ${\eta }_{X_{-1}}:=𝔼\left[{ℐ}_{M_0}\left[\psi \left(u_0\right)\right]\right]$ with tolerances ${\eta }_{X_0}$ and ${\eta }_{Y_0}:={\eta }_{X_0}$.
5.	Set stochastic tolerances:
	${\eta }_{Y_{K-k}}:={\left(2C_Y{G}_K\left({\mu }^{\ast }\right)\right)}^{-1}{\left(F_k\left(s^{\ast }\right)\right)}^{{\mu }^{\ast }/\left({\mu }^{\ast }+1\right)}{\eta }_{X_{k-1}}ϵ$, $k=0,\dots ,K$,
	with $F_k$ and $G_K$ defined in (20) and (21), respectively.
6.	Coarsest level $k=0$ (reusing samples from Step 4):
	Compute $E_0:=𝔼\left[{ℐ}_{M_K}\left[\psi \left(u_0\right)\right]\right]$ with tolerances ${\eta }_{X_0}$ and ${\eta }_{Y_K}$.
7.	Multilevel differences $k=1,\dots ,K$ (reusing samples from level $k-1$):
	Compute $E_k:=𝔼\left[{ℐ}_{M_{K-k}}\left[\psi \left(u_k\right)-\psi \left(u_{k-1}\right)\right]\right]$ with tolerances ${\eta }_{X_k}$, ${\eta }_{X_{k-1}}$ and ${\eta }_{Y_{K-k}}$.
8.	Compute $𝔼\left[{\psi }_K^{\left(\mathrm{ML}\right)}\right]:={\sum }_{k=0,\dots ,K}{E}_k$.

Table 2. Table 2: Adaptive multilevel collocation method: Sequences of spatial and stochastic tolerances, η X k subscript 𝜂 subscript 𝑋 𝑘 {\eta_{X_{k}}} and η Y k subscript 𝜂 subscript 𝑌 𝑘 {\eta_{Y_{k}}} , for overall tolerances ϵ = 10 − 2 , … , 2.5 × 10 − 4 italic-ϵ superscript 10 2 … 2.5 superscript 10 4 \epsilon=10^{-2},\ldots,2.5\times 10^{-4} and certain approximations of η X − 1 subscript 𝜂 subscript 𝑋 1 {\eta_{X_{-1}}} (above). Number of collocation points taken by the anisotropic Smolyak algorithm for K = 2 𝐾 2 K\!=\!2 (below).

$ϵ$	${10}^{-2}$	$5\times {10}^{-3}$	$2.5\times {10}^{-3}$	${10}^{-3}$	$5\times {10}^{-4}$	$2.5\times {10}^{-4}$
${\eta }_{X_{-1}}$	$4.7125$	$4.7346$	$4.7465$	$4.7525$	$4.7554$	$4.7566$
${\eta }_{X_0}$	$1.3\times {10}^{-1}$	$6.3\times {10}^{-2}$	$3.1\times {10}^{-2}$	$1.3\times {10}^{-2}$	$6.3\times {10}^{-3}$	$3.1\times {10}^{-3}$
${\eta }_{X_1}$	$2.5\times {10}^{-2}$	$1.3\times {10}^{-2}$	$6.3\times {10}^{-3}$	$2.5\times {10}^{-3}$	$1.3\times {10}^{-3}$	$6.3\times {10}^{-4}$
${\eta }_{X_2}$	$5.0\times {10}^{-3}$	$2.5\times {10}^{-3}$	$1.3\times {10}^{-3}$	$5.0\times {10}^{-4}$	$2.5\times {10}^{-4}$	$1.3\times {10}^{-4}$
${\eta }_{Y_2}$	$7.2\times {10}^{-3}$	$4.1\times {10}^{-3}$	$2.3\times {10}^{-3}$	$1.1\times {10}^{-3}$	$6.2\times {10}^{-4}$	$3.5\times {10}^{-4}$
${\eta }_{Y_1}$	$1.3\times {10}^{-2}$	$6.1\times {10}^{-3}$	$3.0\times {10}^{-3}$	$1.1\times {10}^{-3}$	$5.5\times {10}^{-4}$	$2.6\times {10}^{-4}$
${\eta }_{Y_0}$	$3.0\times {10}^{-2}$	$1.5\times {10}^{-2}$	$7.2\times {10}^{-3}$	$2.8\times {10}^{-3}$	$1.3\times {10}^{-3}$	$6.4\times {10}^{-4}$
$M_0$	$35$	$51$	$65$	$131$	$161$	$321$
$M_1$	$33$	$33$	$33$	$33$	$33$	$33$
$M_2$	$33$	$33$	$33$	$33$	$33$	$33$

Equations

$$a_y(u(\cdot ,y),v)$$

$$u_{M,h}^{(\mathrm{SL})}(x,y)=\mathcal{I}[u_h](x,y)=\sum _{m=1}^M{u}_m(x){\varphi }_m(y),$$

$$L_{\rho }^2(\Gamma ;H_0^1(D))=\{v:\Gamma \to H_0^1(D)\text{mbox}measurable:{\int }_{\Gamma }\Vert v(y,\cdot ){\Vert }_{H_0^1(D)}^2\rho (y)dy<\infty \}$$

$$\Vert v{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}^2={\int }_{\Gamma }\Vert v(y,\cdot ){\Vert }_{H_0^1(D)}^2\rho (y)dy=:\mathbb{E}\left[\Vert v(y,\cdot ){\Vert }_{H_0^1(D)}^2\right].$$

$$\mathsf{SOLVE}⟹\mathsf{ESTIMATE}⟹\mathsf{MARK}⟹\mathsf{REFINE}$$

$$1\ge {\eta }_{X_0}>{\eta }_{X_1}>\cdots >{\eta }_{X_k}>\cdots >{\eta }_{X_K}>0.$$

$$V_0(y)\subset V_1(y)\subset \cdots \subset V_K(y)\subset H_0^1(D).$$

$$\Vert u(\cdot ,y)-u_k(y){\Vert }_{H_0^1(D)}\le C_k(y)\cdot {\eta }_{X_k},k=0,\dots ,K,$$

$$\Vert u-u_k{\Vert }_{L_{\rho }^2(\Gamma ,H_0^1(D))}\le C_X\cdot {\eta }_{X_k},k=0,\dots ,K,$$

$$C_X:=\underset{k=0,\dots ,K}{\max }{\left({\int }_{\Gamma }{C}_k^2(y)\rho (y)dy\right)}^{1/2}$$

$${\mathcal{I}}_{M_k}:C^0(\Gamma )\to L_{\rho }^2(\Gamma )$$

$$\Vert (u_k-u_{k-1})-{\mathcal{I}}_{M_{K-k}}[u_k-u_{k-1}]{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}\le C_Y\cdot {\eta }_{Y_{K-k}},k=0,\dots ,K,$$

$$u_K^{(\mathrm{ML})}=\sum _{k=0}^K{\mathcal{I}}_{M_{K-k}}[u_k-u_{k-1}]=\sum _{k=0}^K\left(u_{M_{K-k},k}^{(\mathrm{SL})}-u_{M_{K-k},k-1}^{(\mathrm{SL})}\right).$$

$$\Vert u-u_K^{(\mathrm{ML})}{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}\le \Vert u-u_K{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}+\Vert u_K-u_K^{(\mathrm{ML})}{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}.$$

$$\begin{array}[]{rll}\|u_{K}-u_{K}^{(\mathrm{ML})}\|_{L^{2}_{\rho}(\Gamma;H^{1}_{0}(D))}&\!=\!&\displaystyle\left\|\sum_{k=0}^{K}(u_{k}-u_{k-1})-{\cal I}_{M_{K-k}}[u_{k}-u_{k-1}]\right\|_{L^{2}_{\rho}(\Gamma;H^{1}_{0}(D))}\\[14.22636pt] &&\leq\displaystyle\sum_{k=0}^{K}\left\|(u_{k}-u_{k-1})-{\cal I}_{M_{K-k}}[u_{k}-u_{k-1}]\right\|_{L^{2}_{\rho}(\Gamma;H^{1}_{0}(D))}\leq\sum_{k=0}^{K}C_{Y}\cdot{\eta_{Y_{K-k}}}.\end{array}$$

$$\Vert u-u_K^{(\mathrm{ML})}{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}\le 2C_X\cdot {\eta }_{X_K},$$

$$\begin{array}[]{lrll}\text{(A1)}&\quad A_{k}&\leq&C_{A}\cdot{\eta_{X_{k}}}\!\!{}^{-s},\\[5.69054pt] \text{(A2)}&\quad C_{Y}\cdot{\eta_{Y_{K-k}}}&=&C_{I}(N)\,M_{K-k}^{-\mu}\,{\eta_{X_{k-1}}}\end{array}$$

$$\Vert u_0{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}\le {\eta }_{X_{-1}}:=const.$$

$$C^{(\mathrm{ML})}\le \sum _{k=0}^K{M}_{K-k}(A_k+A_{k-1}),$$

$$\Vert u-u_K^{(\mathrm{ML})}{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}\le ϵ$$

$$\displaystyle C_{\epsilon}^{(\mathrm{ML})}\leq C\cdot\big{(}G_{K}(\mu)\big{)}^{\frac{\mu+1}{\mu}}\,\epsilon^{-\frac{1}{\mu}}+C_{A}\,\sum_{k=0}^{K}F_{k}(s)\,{\eta_{X_{k-1}}}$$

$$F_k(s)$$

$$G_K(\mu )$$

$${\eta_{Y_{K-k}}}=\big{(}2\,C_{Y}\,G_{K}(\mu)\big{)}^{-1}\big{(}F_{k}(s)\big{)}^{\frac{\mu}{\mu+1}}\,{\eta_{X_{k-1}}}\,\epsilon.$$

$$\underset{M_0,\dots ,M_K}{\min }\sum _{k=0}^K{C}_A\cdot M_{K-k}\left({\eta }_{X_k}{}^{-s}+{\eta }_{X_{k-1}}{}^{-s}\right)s.t.\sum _{k=0}^K{C}_I\cdot M_{K-k}^{-\mu }{\eta }_{X_{k-1}}=\frac{ϵ}{2}.$$

$$M_{K-k}=\big{(}2\,C_{I}\,G_{K}(\mu)\big{)}^{\frac{1}{\mu}}\big{(}F_{k}(s)\big{)}^{-\frac{1}{\mu+1}}\,\epsilon^{-\frac{1}{\mu}}$$

$$C_{ϵ}^{(\mathrm{ML})}$$

$$\displaystyle=\,\sum_{k=0}^{K}C_{A}\cdot\left(\big{(}2\,C_{I}\,G_{K}(\mu)\big{)}^{\frac{1}{\mu}}\,\big{(}F_{k}(s)\big{)}^{-\frac{1}{\mu+1}}\,\epsilon^{-\frac{1}{\mu}}+1\right)F_{k}(s)\,{\eta_{X_{k-1}}}$$

$$\displaystyle\leq C\cdot\big{(}G_{K}(\mu)\big{)}^{\frac{1}{\mu}}\epsilon^{-\frac{1}{\mu}}\sum_{k=0}^{K}\big{(}F_{k}(s)\big{)}^{\frac{\mu}{\mu+1}}\,{\eta_{X_{k-1}}}+C_{A}\,\sum_{k=0}^{K}F_{k}(s)\,{\eta_{X_{k-1}}}$$

$$\displaystyle=C\cdot\big{(}G_{K}(\mu)\big{)}^{\frac{\mu+1}{\mu}}\,\epsilon^{-\frac{1}{\mu}}+C_{A}\,\sum_{k=0}^{K}F_{k}(s)\,{\eta_{X_{k-1}}}$$

$${\eta_{Y_{K-k}}}=\big{(}2\,C_{Y}\,G_{K}(\mu)\big{)}^{-1}\,\big{(}F_{k}(s)\big{)}^{\frac{\mu}{\mu+1}}\,{\eta_{X_{k-1}}}\,\epsilon.$$

$$\Vert u-u_K^{(\mathrm{ML})}{\Vert }_{L_{\rho }^2(\Gamma ;H_0^1(D))}\le ϵ$$

Full text

Keywords. multilevel methods, hierarchical methods, adaptivity, stochastic collocation, PDEs with random data, sparse grids, uncertainty quantification, high-dimensional approximation.

AMS subject classifications. 65C20, 65C30, 65N35, 65M75

1 Introduction

A central task in computational science and engineering is the efficient numerical treatment of partial differential equations (PDEs) with uncertain input data (coefficients, source term, geometry, etc). This has been a very active area of research in recent years with a host of competing ideas. Our focus here is on elliptic PDEs with random input data in a standard (pointwise) Hilbert space setting: Find $u(\cdot,y)\in V$ such that

[TABLE]

where $a_{y}(\cdot,\cdot):V\times V\to\mathbb{R}$ is a parameter-dependent inner product in a Hilbert space $V$ and $\ell_{y}(v):V\to\mathbb{R}$ is a parameter-dependent bounded linear functional. In general, the solution will be represented as $u(x,y):D\times\Gamma\rightarrow\mathbb{R}$ where $D\subset\mathbb{R}^{d}$, $d=1,2,3$, is the deterministic, bounded (physical) domain and $\Gamma=\Gamma_{1}\times\Gamma_{1}\times\cdots\Gamma_{N}$ is a (stochastic) parameter space of finite dimension $N$ (finite noise assumption). The component parameters $y_{1},\ldots,y_{n}$ will be associated with independent random variables that have a joint probability density function $\rho(y)=\Pi_{n=1}^{N}\hat{\rho}(y_{n})\in L^{\infty}(\Gamma)$ such that $\hat{\rho}_{n}:[-1,1]\rightarrow\mathbb{R}$.

The most commonly used approach to solve (1) is the Monte Carlo (MC) method based on random samples $y^{(m)}\in\Gamma$. It has a dimension-independent convergence rate and leads to a natural decoupling of the stochastic and spatial dependency. While MC methods are well suited to adaptive spatial refinement, they are not able to exploit any smoothness or special structure in the parameter dependence. The standard, single-level stochastic collocation method (see, for example, [2, 31, 32, 37]) for (1) is similar to the MC method in that it involves independent, finite-dimensional spatial approximations $u_{h}(y^{(m)})\approx u(y^{(m)})$ at a set $\{y^{(m)}\}_{m=1,\ldots,M}$ of deterministic sampling points in $\Gamma$, so as to construct an interpolant

[TABLE]

in the polynomial space ${\cal P}_{M}=\mbox{span}\{\phi_{m}\}_{m=1,\ldots,M}\subset L^{2}_{\rho}(\Gamma)$ with basis functions $\{\phi_{m}\}_{m=1,\ldots,M}$. The coefficients $u_{m}(x)$ are determined by the interpolating condition ${\cal I}[u_{h}](x,y^{(m)})=u_{h}(x,y^{(m)})$, for $m=1,\ldots,M$. The quality of the interpolation process depends on the accuracy of the spatial approximations $u_{h}(y^{(m)})$ and on the number of collocation points $M$, which typically, grows rapidly with increasing stochastic dimension $N$.

Two concepts in the efficient treatment of (1) that have shown a lot of promise and are (to some extent) complementary are (i) adaptive sparse-grid interpolation methods that are able to exploit smoothness with respect to the uncertain parameters and (ii) multilevel approaches that aim at reducing the computational cost through a hierarchy of spatial approximations. Adaptive sparse-grid methods, where the set of sample points is adaptively generated, can be traced back to Gerstner & Griebel [23] and have been extensively tested in collocation form; see [1, 11, 25, 30, 33]. We note that parametric adaptivity has also been explored in a Galerkin framework; see [7, 9, 13, 16, 17], and that there are a number of recent papers aimed at proving dimension-independent convergence; see [3, 8, 29, 40, 41].

The multilevel Monte Carlo (MLMC) method was originally proposed as an abstract variance reduction technique by Heinrich [26], and independently for stochastic differential equations in mathematical finance by Giles [24]. In the context of uncertainty quantification it was developed by Barth et al. [4] and Cliffe et al. [12] and extended to stochastic collocation sampling methods soon thereafter by Teckentrup et al. [36]. The methodology in this paper can be viewed as an extension of this body of work.

Problems where the PDE solution develops singularities (or at least steep gradients) in random locations in the spatial domain are a challenging feature in many real-world applications. Spatial adaptivity driven by a posteriori error estimates has been investigated within single-level stochastic collocation in [33] and has been considered in a MLMC framework in [18]. To date, however, MC-based methods [14, 19, 27, 39] are the only methods where the spatial approximation is locally refined at each parametric sample point. In what follows, we will combine a pointwise adaptive mesh refinement for the approximation of the spatial approximations $u_{h}(y^{(m)})$ with an adaptive (anisotropic) sparse Smolyak grid, in order to improve the efficiency of the multilevel method to reach a user-prescribed tolerance for the accuracy of the multilevel interpolant, as well as for the computation of associated quantities of interest. By rigorous control of both error components, the fully adaptive algorithm proposed herein is able to exploit smoothness and structure in the parameter dependence while also exploiting fully the computational gains due to the spatially adapted hierarchy of spatial approximations at each collocation point. Under certain assumptions on the efficiency of the two-component adaptive schemes, the complexity of the new multilevel algorithm can be estimated rigorously. Two specifically chosen numerical examples will be looked at in the final section to demonstrate its efficiency.

2 Methodology

Following the paper of Teckentrup et al. [36], we consider steady-state diffusion problems with uncertain parameters. Thus, we will assume that $V:=H^{1}_{0}(D)$ and that the variational formulation (1) admits a unique solution in the weighted Bochner space

[TABLE]

with corresponding norm

[TABLE]

The structure of our adaptive algorithm is identical to that introduced in [36]. The distinctive aspects are discussed in this section.

2.1 Adaptive spatial approximation

We will focus on finite element approximation in the spatial domain $D$ using classical adaptive methodology. That is, we execute the loop

[TABLE]

until the estimate of the error in the second step is less than a prescribed tolerance. Next, we let $\{{\eta_{X_{k}}}\}_{k=0,...,K}$ be a decreasing sequence of tolerances with

[TABLE]

Then, for each fixed parameter $y\in\Gamma$, we run our adaptive strategy to compute approximate spatial solutions $u_{k}(y)\in V_{k}(y)$ on a sequence of nested subspaces

[TABLE]

This is the point of departure from the previous work [36], wherein the sequence of approximation spaces $V_{0}\subset V_{1}\subset\cdots\subset V_{K}$ could be chosen in a general way, but in a manner that was fixed for all collocation points $y\in\Gamma$.

Starting from the pointwise error estimate

[TABLE]

and supposing measurability of the discrete spaces $V_{k}(y)$ and bounded second moments of $C_{k}(y)$, we directly get the following error bound

[TABLE]

with a constant

[TABLE]

that does not depend on $y$ and $k$. Adaptive algorithms proposed by Dörfler [15] and Kreuzer[28] converge for fixed $y\in\Gamma$ and ${\eta_{X_{k}}}\rightarrow 0$. The constant $C_{X}$ in (10) is related to the effectivity of the a posteriori error estimation strategy. Values close to one can be obtained using hierarchical error estimators (see, for example, Bespalov et al. [10]) or using gradient recovery techniques in the asymptotic regime.

2.2 Adaptive stochastic interpolation

Let us assume $u\in C^{0}(\Gamma;H^{1}_{0}(D))$ and denote by $\{{\cal I}_{M_{k}}\}_{k=0,1,\ldots}$ a sequence of interpolation operators

[TABLE]

with $M_{k}$ points from the $N$-dimensional space $\Gamma$. We construct each of these operators by a hierarchical sequence of one-dimensional Lagrange interpolation operators with the anisotropic Smolyak algorithm, which was introduced by Gerstner & Griebel [23]. The method is dimension adaptive, using the individual surplus spaces in the multi-dimensional hierarchy as natural error indicators.

Let $\{{\eta_{Y_{k}}}\}_{k=0,...,K}$ be a second sequence of tolerances, a priori not necessarily decreasing. Under suitable regularity assumptions for the uncertain data (see e.g. Babuška et al. [2, Lem. 3.1, 3.2]), we can assume that there exist numbers $M_{k}$, $k=0,1,\ldots,K$, and a constant $C_{Y}>0$ not depending on $k$ such that

[TABLE]

where, for simplicity, we set $u_{-1}=0$.

Since (8) implies $\|u_{k}-u_{k-1}\|_{L^{2}_{\rho}(\Gamma;H^{1}_{0}(D))}\leq C\cdot{\eta_{X_{k-1}}}\,$, which is decreasing as $k\rightarrow\infty$, we can expect that for higher $k$ it suffices to use less accurate interpolation operators, i.e., smaller numbers $M_{K-k}$, to achieve a required accuracy. Indeed this is the main motivation to set up a multilevel interpolation approximation. Note that in this way, the tolerances ${\eta_{Y_{K-k}}}$ are strongly linked to the spatial tolerances ${\eta_{X_{k}}}$. We will make this connection more precise and give suitable values for the sequence of tolerances $\{{\eta_{Y_{k}}}\}_{k=0,\ldots,K}$ in the next section.

2.3 An adaptive multilevel strategy

Given the sequences $\{u_{k}\}$ and $\{{\cal I}_{M_{k}}\}$, we define the multilevel interpolation approximation in the usual way by

[TABLE]

Observe that the most accurate interpolation operator ${\cal I}_{M_{K}}$ is used on the coarsest spatial approximation $u_{0}$ whereas the least accurate interpolation operator ${\cal I}_{M_{0}}$ is applied to the difference of the finest spatial approximations $u_{K}-u_{K-1}$. The close relationship between the spatial and stochastic approximations at index $k$ is clearly visible in (13). We note that the use of adaptive interpolation operators is also mentioned in [36]. The key difference is that the hierarchy in [36] is based on the number of collocation points $M_{k}$ rather than the tolerances $\{{\eta_{Y_{k}}}\}$ employed here.

To show the convergence of the multilevel approximation $u_{K}^{(\mathrm{ML})}$ to the true solution $u$, we split the error into the sum of a spatial discretization error and a stochastic interpolation error using the triangle inequality

[TABLE]

Due to (9) the first term on the right hand side of (14) is bounded by $C_{X}\cdot{\eta_{X_{K}}}$. The aim is now to choose the tolerances $\{{\eta_{Y_{k}}}\}$ in an appropriate way to reach the same accuracy for the second term. From (12), we estimate the stochastic interpolation error as follows:

[TABLE]

To obtain an accuracy of the same size as the spatial discretization error, a first choice for the tolerances would be to simply demand ${\eta_{Y_{k}}}\leq C_{X}\cdot{\eta_{X_{K}}}/((K+1)C_{Y})$, for all $k=0,\ldots,K$, and conclude that

[TABLE]

so that the adaptive multilevel method converges for ${\eta_{X_{K}}}\rightarrow 0$. However, the values for ${\eta_{Y_{k}}}$ can be optimized by minimizing the computational costs while keeping the desired accuracy.

At this point it is appropriate to consider the computational cost, $C_{\epsilon}^{(\mathrm{ML})}$, of the multilevel stochastic collocation estimator $u_{K}^{(\mathrm{ML})}$ required to achieve an accuracy $\epsilon$. Thus, in order to quantify the contributions from the spatial discretization and the stochastic collocation, we will need to make two assumptions to link the cost with the error bounds in (9) and (12). Let $A_{k}$ denote an upper bound for the cost to solve the deterministic PDE at sample point $y\in\Gamma$ with accuracy ${\eta_{X_{k}}}$. Then, we assume that for all $k=0,\ldots,K$,

[TABLE]

as well as the special case

[TABLE]

Here, the constants $C_{A}>0$, $C_{I}(N)>0$ are independent of $y$, $k$, and the rates $s,\mu>0$. Note that we could also consider the exact cost per sample and introduce a sample dependent constant $C_{A}(y)$ in (A1), but that would complicate the subsequent analysis.

Assumption (A1) usually holds for first-order adaptive spatial discretization methods with $s=d$, when coupled with optimal linear solvers such as multigrid. The factors on the right-hand side in (A2) reflect best the convergence of the sparse grid approximations in (12) with respect to the total number $M_{K-k}$ of collocation points, see Nobile et al. [31, 32] or [36, Theorem 5.5]. To estimate the difference $u_{k}-u_{k-1}$, we use the fact that $\|u_{k}-u_{k-1}\|_{L^{2}_{\rho}(\Gamma;H^{1}_{0}(D))}\leq C\cdot{\eta_{X_{k-1}}}$ with a constant $C>0$ close to $C_{X}$. It follows from $\|u-u_{k-1}\|_{H^{1}_{0}(D)}\approx\|u_{k}-u_{k-1}\|_{H^{1}_{0}(D)}$, which is the basis for the very good performance of hierarchical error estimators. We will absorb $C$ into $C_{I}$ in the sequel.

The rate $\mu$ depends in general on the dimension $N$. Theoretical results for the anisotropic classical Smolyak algorithm are given in [31, Thm. 3.8]. However, it has recently been shown in [40, 41] that under certain smoothness assumptions $\mu$ can be independent of $N$.

An upper bound for the total computational cost of the approximation $u_{K}^{(\mathrm{ML})}$ can then be defined as

[TABLE]

with $A_{-1}:=0$. In a first step, we will consider a general sequence $\{{\eta_{X_{k}}}\}_{k=0,\ldots,K}$ without defining a decay rate a priori. Following the argument in [36], with a priori estimates for the spatial error replaced by tolerances, then leads to an estimate of $\epsilon$-cost $C_{\epsilon}^{(\mathrm{ML})}$ and identifies a set of optimal tolerances ${\eta_{Y_{k}}}$ in (12).

Theorem 2.1**.**

Let ${\eta_{X_{0}}},{\eta_{X_{1}}},\ldots,$ be a decreasing sequence of spatial tolerances satisfying (6). Suppose assumptions (A1) and (A2) hold. Then, for any $\epsilon$, there exist an integer $K=K(\epsilon)$ and a sequence of tolerances $\{{\eta_{Y_{k}}}\}_{k=0,\ldots,K}$ in (12) such that

[TABLE]

and

[TABLE]

with $C=C_{A}\,(2\,C_{I})^{\frac{1}{\mu}}$ and

[TABLE]

where for ease of notation we set ${\eta_{X_{-1}}}=\infty$. The optimal choice for the tolerances ${\eta_{Y_{k}}}$ is then given by

[TABLE]

The utility of (22) is that near-optimal tolerances can be readily computed if estimates of the constants and the rates in (A1) and (A2) are available.

Proof: As in the convergence analysis above, we split the error and make sure that both the spatial discretization error and the stochastic interpolation error are bounded by $\epsilon/2$. First, we choose an appropriate $K\geq 0$ and ${\eta_{X_{k}}}$ such that $C_{X}\cdot{\eta_{X_{k}}}<\epsilon/2$. This is, of course, always possible and fixes the number $K=K(\epsilon)$ as a function of $\epsilon$. Next we determine the set $\{M_{k}\}_{k=0,\ldots,K}$ so that the computational cost in (17) is minimized subject to the requirement that the stochastic interpolation error is bounded by $\epsilon/2$. Using assumptions (A1) and (A2), this reads

[TABLE]

Application of the Lagrange multiplier method with all $M_{k}$ treated as continuous variables as in Giles [24] gives the optimal choice for the number of samples

[TABLE]

with $F_{k}$ and $G_{K}$ defined in (20) and (21). To ensure that $M_{K-k}$ is an integer, we round up to the next integer. The $\epsilon$–cost of the multilevel approximation can then be estimated as follows

[TABLE]

with C=C_{A}\,\big{(}2\,C_{I}\big{)}^{\frac{1}{\mu}}.

The optimal tolerances ${\eta_{Y_{k}}}$ can be directly determined from assumption (A2). Thus with $M_{k}$ defined above we get

[TABLE]

Note that with these values $\sum_{k=0}^{K}C_{Y}\cdot{\eta_{Y_{k}}}=\epsilon/2$, which gives the desired accuracy in (15). $\Box$

Observe that the function $G_{K}(\mu)$ as well as the second term in (19) still depend on $\epsilon$, because $K$ is a function of $\epsilon$. In this way, the choice of the tolerances ${\eta_{X_{k}}}$ has an influence on the rate $-1/\mu$, which could be further optimized. However, in the following we will restrict our attention to a typical geometric design with ${\eta_{X_{k}}}=q^{k}\,{\eta_{X_{0}}},\;k=1,2,\ldots,$ with a positive reduction factor $q<1$. The overall cost can then be estimated using a standard construction, leading to the following result.

Theorem 2.2**.**

Let the sequence of spatial tolerances $\{{\eta_{X_{k}}}\}_{k=0,1,\ldots,K}$ in (6) be defined by ${\eta_{X_{k}}}=q^{k}\,{\eta_{X_{0}}}$ with a reduction factor $q<1$. Suppose assumptions (A1) and (A2) hold. Then, for any $\epsilon<1$, there exists an integer $K=K(\epsilon)$ such that

[TABLE]

and

[TABLE]

Proof: See [35, 36]. $\Box$

In typical applications, it is usually more natural to consider a functional $\psi$ of the solution $u$ instead of the solution itself. Thus, suppose a (possibly nonlinear) functional $\psi:H^{1}_{0}(D)\rightarrow\mathbb{R}$ with $\psi(0)=0$ is given. In this case, we define the following single-level and multi-level stochastic collocation approximations:

[TABLE]

with $u_{-1}:=0$. As in (9) and (12), we can ensure that for the adaptive error control of the expected values, for all $k=0,\ldots,K,$ we have

[TABLE]

In practice, the following analogue of Theorem 2.2 holds for the expected value of the error of the multilevel approximation of functionals.

Proposition 2.1**.**

Let the sequence of spatial tolerances $\{{\eta_{X_{k}}}\}_{k=0,1,\ldots,K}$ in (6) be defined by ${\eta_{X_{k}}}=q^{k}\,{\eta_{X_{0}}}$ with a reduction factor $q<1$. Suppose assumptions (A1) and (A2) hold with convergence rates $\mu^{*}$ and $s^{*}$. Then, for any $\epsilon<1$, there exists an integer $K(\epsilon)$ such that

[TABLE]

and

[TABLE]

Proof: see the discussion of Proposition 4.6 in [36].

Note that in analogy to strong and weak convergence of solutions of stochastic differential equations, one would expect the error in mean associated with approximating a functional $\psi(u)$ to decrease at a faster rate than the error in norm associated with approximating the entire solution $u$. Thus, we anticipate that $\mu^{*}>\mu$. Moreover, the convergence rate of the error in the functional $\psi(u_{k})$ is in general larger than the convergence rate of the error in the $H_{1}$-norm, which leads to a smaller value $s^{*}<s$ in assumption (A1). Specifically, in the case of $H^{2}$-regularity in space and using an optimal linear solver, such as multigrid, we anticipate that ${s^{*}}=d/2$.

We conclude this section with a complete algorithmic description of our adaptive multilevel stochastic collocation method. Table 1 illustrates the main steps. The strategy is self-adaptive in nature. Thus, once the tolerances $\{{\eta_{X_{k}}}\}_{k=0,\ldots,K}$ and $\{{\eta_{Y_{k}}}\}_{k=0,\ldots,K}$ are set, the algorithm delivers an approximate functional $\psi_{K}^{(\mathrm{ML})}$ with accuracy close to the user-prescribed tolerance $\epsilon$. An analogous algorithm can be defined to deliver the numerical solution $u_{K}^{(\mathrm{ML})}$ close to a user-prescribed tolerance $\epsilon$. A priori, to obtain optimal results, the reliability of the estimation for the adaptive spatial discretization and the adaptive Smolyak algorithm needs to be studied in order to provide values for the constants $C_{X}$ and $C_{Y}$. More specifically, while the spatial tolerances ${\eta_{X_{k}}}$ can be freely chosen (by fixing the number of levels $K$ and the reduction factor $q$), the optimal choice of the tolerances ${\eta_{Y_{k}}}$ in (22) requires estimates of the parameters $s^{*}$ and $\mu^{*}$. A discussion of how to effect this in a preprocessing step is given in the next section. A crucial point to note is that, even without this information, the adaptive anisotropic Smolyak algorithm will automatically detect the importance of various directions in the parameter space $\Gamma\subset\mathbb{R}^{N}$.

An important point already mentioned in [36] is that the optimal rounded values for the number of samples, $M_{k}$, will not be used by the algorithm, because they do not necessarily correspond to an adaptive sparse grid level. However, for each level $k$, the tolerance ${\eta_{Y_{k}}}$ can be ensured by choosing $\widetilde{M}_{k}\geq M_{k}$ slightly larger, resulting in a slight inefficiency of the sparse grid approximation. Note that, in practice, the same behaviour is observed for adaptive spatial discretizations. In any case, no restart is needed; cf. [36, Section 6.3].

3 Numerical examples

First, we provide general information on the adaptive components used. Then, numerical results are presented for two isotropic diffusion problems with uncertain source term and uncertain geometry, respectively. All calculations have been done with Matlab version R2017a on a Latitude $7280$ with an i5-7300U Intel processor running at 2.7 GHz.

For the spatial approximation considered in the two examples, we use the adaptive piecewise linear finite element method implemented by Funken, Praetorius and Wissgott in the Matlab  package p1afem.222The p1afem software package can be downloaded from the author’s webpage [20]. Using Matlab built-in functions and vectorization for an efficient realization, the code performs with almost linear complexity in terms of degrees of freedom with respect to the runtime. A general description of the underlying ideas can be found in [22] and the technical report [21] provides detailed documentation. The code is easy to modify. In order to control the accuracy of solution functionals, the dual weighted residual method (DWRM) introduced by Becker & Rannacher [6] is adopted. In what follows, we will give a short summary of the underlying principles that are relevant for our implementation.

Let $y\in\Gamma$ be fixed. For ease of presentation we write $u=u(\cdot,y)$ and consider the variational problem (1) with $a(u,v)=\int_{D}\nabla u\cdot\nabla v\,dx$. We also suppose the solution functional $\psi(u)$ takes the specific form

[TABLE]

as discussed in [36, Example 5.13]. Let $u_{k}\in V_{k}$ be the finite element solution computed on the (adaptive) mesh $\mathcal{T}_{k}$. Then the DWRM provides a representation of the error in the solution functional in the form

[TABLE]

by simply neglecting the higher-order term. Letting $w$ be the exact solution of the linearised dual problem

[TABLE]

and letting $w_{k}\in V_{k}$ be the finite element approximation of the dual solution on the same mesh, we have, using Galerkin orthogonality,

[TABLE]

In practice, when solving the dual problem we compute an approximation $\phi_{k}\approx w-w_{k}$ of its error in the hierarchical surplus space of quadratic finite elements. Hierarchical error estimators using this approach are implemented in p1afem for the primal solution $u_{k}$.

Putting this all together, we obtain

[TABLE]

We use $|\eta_{T}|$ as refinement indicators and mark elements $T\in\mathcal{T}_{k}$ for refinement using the standard Dörfler criterion from [15], which determines the minimal set $\mathcal{M}\subset\mathcal{T}_{k}$ such that $\theta\sum_{T\in\mathcal{T}_{k}}|\eta_{T}|\leq\sum_{T\in\mathcal{M}}|\eta_{T}|$. We typically set a value of $\theta=0.6$ in our calculations. Refinement by newest vertex bisection is applied to guarantee nested finite element spaces and the optimal convergence of the adaptive finite element method, see [28]. The adaptive process is terminated when the absolute value of $\sum_{T\in\mathcal{T}_{k}}\eta_{T}$ is less than a prescribed tolerance.

Turning to the adaptive anisotropic Smolyak algorithm, the main idea for the construction of the sparse grid interpolation operators in (11) is to use the hierarchical decomposition

[TABLE]

with multi-indices ${\bf i}=(i_{1},\ldots,i_{N})\in I\subset\mathbb{N}^{N}_{+}$, $m({\bf i})=(m(i_{1}),\ldots,m(i_{N}))$, and univariate polynomial interpolation operators ${\cal I}_{n}^{m(i_{n})}:C^{0}(\Gamma_{n})\rightarrow\mathbb{P}_{m(i_{n})-1}$, which use $m(i_{n})$ collocation points to construct a polynomial interpolant in $y_{n}\in\Gamma_{n}$ of degree at most $m(i_{n})-1$. The operators $\triangle^{m({\bf i})}$ are often referred to as hierarchical surplus operators. The function $m$ has to satisfy $m(0)=0$, $m(1)=1$, and $m(i)<m(i+1)$. We set ${\cal I}^{0}_{n}=0$ for all $n=1,\ldots,N$ and use the nested sequence of univariate Clenshaw–Curtis nodes with $m(i)=2^{i-1}+1$ if $i>1$. In (40), $M_{k}$ is then the number of all explored quadrature points in $\Gamma$ determined by $m({\bf i})$. To get good approximation properties, the index set $I$ should satisfy the downward closed set property, i.e.,

[TABLE]

As usual, we ensure that ${\bf 1}\in I$ to also recover constant functions.

The hierarchical structure in (40) allows to interpret updates that are derived by adding further differences $\triangle^{m({\bf i}_{a})}$, i.e., enhancing the index set $I$ by an admissible index ${\bf i}_{a}$ that satisfies (41), as error indicators for already computed approximations. There are several adaptive strategies available. One could explore the whole margin of $I$ defined by

[TABLE]

Generally, this approach is computationally challenging and yields a fast increase of quadrature points. Instead, as suggested by Gerstner & Griebel [23], the margin is reduced to the set

[TABLE]

In each step, the adaptive Smolyak algorithm computes the profits $\triangle^{m({\bf i}_{a})}$ for all ${\bf i}_{a}\in R_{I}$ – reusing already computed profits – and replaces the index in $R_{I}$ with the highest profit, say ${\bf i}_{max}$, by its admissible neighbours taken from the set $\{{\bf i}_{max}+{\bf e}_{j},j=1,\ldots,N\}$. These neighbours are then explored next. The algorithm stops if the absolute value of the highest profit is less than a prescribed tolerance. This adaptive strategy is implemented in Matlab in the Sparse Grid Kit333This package can be downloaded from the CSQI website [34]. and its numerical performance is discussed in the review paper [5].

We adopt this refinement strategy with one minor deviation in our numerical experiments: namely, at the final iteration step, instead of adding the new profits to the interpolant, the previous approximate value is returned. In that way, the last highest profit can be considered as a more realistic error indicator. This adaptation allows a better understanding of the convergence behaviour of the multilevel approach. In practical calculations one would almost certainly use the final value.

3.1 Uncertain source and boundary conditions: $d=2$, $N=2$

We will refer to this example as the one peak test problem. It was introduced by Kornhuber & Youett [27, Sect. 5.1] in order to assess the efficiency of adaptive multilevel Monte Carlo methods. The challenge is to solve the Poisson equation $-\nabla^{2}u=f$ in a unit square domain $D=(-1,1)\times(-1,1)$ with Dirichlet boundary data $u=g$ on $\partial D$ (points in $D$ are represented by ${{x}}=(x_{1},x_{2})$). The source term $f$ and boundary data are uncertain and are parameterised by $y=(y_{1},y_{2})$, representing the image of a pair of independent random variables with $y_{j}\sim{\cal U}[-1/4,1/4]$. In the isotropic case studied in [27], the source term $f$ and boundary data $g$ are chosen so that the uncertain PDE problem has a specific pathwise solution given by

[TABLE]

where the scaling factor $\beta>0$ is chosen to be large—so as to generate a highly localised Gaussian profile centered at the uncertain spatial location $(y_{1},y_{2})$. The value of $\beta$ that we will take in this study is $\beta=50$. (The other values discussed in  [27] are $\beta=10$ and $\beta=150$.)

In this work, the test problem is made anisotropic by scaling the solution in the first coordinate direction by a linear function $\alpha(y_{1})=18y_{1}+11/2$ so that the $\alpha$ takes values in the interval $[-1,10]$. The corresponding pathwise solution is then given by

[TABLE]

The goal is to approximate the following quantity of interest (QoI)

[TABLE]

In our experiments we solve the simplified problem with $u=0$ on $\partial D$. (This has no impact on accuracy for the error tolerances that we will considered; see the discussion in the appendix). The initial mesh with $81$ vertices is generated by taking three uniform refinement steps, starting from two triangles with common edge from $(-1,-1)$ to $(1,1)$. For the specific parameter choice $y=(-0.22,-0.22)^{T}$ and spatial tolerance ${\eta_{X}}=2\times 10^{-3}$, the adaptive algorithm terminated after 5 steps giving the numerical solution and corresponding mesh shown in Fig. 1.

Running the adaptive algorithm with a tighter tolerance of ${\eta_{X}}=5\times 10^{-6}$ generated a concentrated mesh with $158\,734$ points after 13 refinement steps. The left part of Fig. 2 shows the efficiency of the locally adaptive procedure in comparison with a uniform refinement strategy. For smaller tolerances, the gain of efficiency in terms of overall cpu time is a factor $10$. Comparing with a reference value $\psi(u)=0.025315675\ldots$, generated by running the adaptive algorithm with a very small tolerance, we see that the effectivity index, i.e., the ratio between estimator and true error, tends asymptotically to $1$ for both uniform and adaptive refinement.

The very good quality of the spatial error estimation process is still maintained if we sample over the whole parameter space using an isotropic sparse grid of $145$ collocation points – generated using Sparse Grid Kit by doubling of points in each stochastic dimension. In this case we have a reference value of $\mathbb{E}[\psi(u)]=0.015095545\ldots$, generated analytically using an asymptotic argument—details are given in the appendix. The error estimates for adaptive meshes associated with different spatial error tolerances are plotted in the right part of Fig. 2. Observe that the estimators deliver upper bounds for the numerical errors and the tolerances are always satisfied. As expected from the theory, the convergence rates for $\psi(u)$ and $\mathbb{E}[\psi(u)]$ in terms of computing time are both close to $-1$. So we have $C_{X}\!=\!1$ in (30) and $s^{*}\!=\!1$ in assumption (A1) in Section 2.

Next we consider the convergence behaviour and the quality of the error estimates for the adaptive anisotropic Smolyak algorithm. A reference solution generated by solving the problem with a small spatial error tolerance on a $33\times 33$ tensor-product grid of Clenshaw–Curtis points is shown in the left part of Fig. 3. Looking at the solution structure, it is evident that the QoI has inherent structure that can be exploited by an adaptive strategy: more specifically, very few sample points are needed to accurately compute it ! The right part of Fig. 3 shows the results for a decreasing sequence of stochastic tolerances ${\eta_{Y}}=10^{-3},\ldots,10^{-6}$, where we have fixed the spatial tolerance ${\eta_{X}}=10^{-7}$. The corresponding numbers of collocation points are $5,7,11$, and $11$. Looking at these more closely we find that all the quadrature rules generated by the adaptive procedure are one-dimensional rules in the $y_{1}$ direction matching the variation in the reference. Other features of note are that the prescribed tolerances are always satisfied and that the errors (with respect to the reference solution) are nicely estimated by the error indicators. A least-squares fit gives an averaged value of $\mu^{*}=9.75$ for the convergence order so that assumption (A2) in Section 2 is satisfied. An approximate value $\mu^{*}\approx 10.74$ can be computed with a much lower spatial tolerance $\eta_{X}=1.25\times 10^{-4}$ and $19$ collocation points in a few seconds.

Having estimated the parameters in the second step of the algorithm in Table 1, we now run the adaptive multilevel approach with overall accuracy requirements of $\epsilon=10^{-5},5\times 10^{-6},2.5\times 10^{-6}$ and $10^{-6}$. To illustrate the performance gains, we simply set $K\!=\!2$ (so that three levels are used) and assign a spatial error reduction factor of $q=0.2$. The spatial tolerances are then given by ${\eta_{X_{k}}}=\epsilon q^{k-2}/2$ with $k=0,1,2$. To calculate, in a first step, a sufficiently accurate approximation for the tolerance ${\eta_{X_{-1}}}=\mathbb{E}[\psi(u_{0})]$ at reasonable cost, we apply the anisotropic Smolyak algorithm with ${\eta_{X}}={\eta_{Y}}={\eta_{X_{0}}}$. Note that these samples can be reused later in the first level of the multilevel scheme. Eventually, the stochastic tolerances are derived from (22) with $C_{Y}=0.1$, $\mu^{*}=9.75$ and $s^{*}=1$.

Results for the three-level and single-level approach are summarised in Fig. 4. For adaptive spatial meshes, we observe that the errors of the expected values $\mathbb{E}[\psi_{2}^{(\mathrm{ML})}]$ are very close to the prescribed tolerances. This is not always the case for the values $\mathbb{E}[\psi_{2}^{(\mathrm{SL})}]$ computed by the single-level approach. The three-level approach performs reliably and outperforms the single-level version clearly. The orders of convergence, $p_{\mathrm{ML}}\!=\!-1/s^{*}\!=\!-1$ and $p_{\mathrm{SL}}\!=\!-1/(s^{*}+1/{\mu^{*}})\!\approx\!-0.91$, predicted by Theorem 2.1 for the accuracy in terms of computational complexity are also visible. Already for coarse tolerances, the three-level approach with uniform spatial meshes shown plotted with inverted blue triangles) is not competitive. For higher tolerances, the calculation exceeded memory requirements (more than $1.6\times 10^{7}$ FE degrees of freedom).

These CPU timings are impressive. We have also solved the one peak test problem using an efficient adaptive stochastic Galerkin (SG) approximation strategy. While the linear algebra associated with the Galerkin formulation is decoupled in this case, the computational overhead of evaluating the right-hand-side vector is a big limiting factor in terms of the relative efficiency. Using SG the quantity

[TABLE]

where

[TABLE]

and

[TABLE]

must be computed in every element $\triangle$ in the current subdivision and for every parametric function $\psi_{\nu_{k}}(y)$ in the active index set. This is an extremely demanding quadrature problem!

We have also tested an adaptive multilevel Monte Carlo method (see, [12, 27]) on this problem. For the lowest tolerance, $\epsilon=10^{-5}$, and an average over $5$ independent realizations, the three-level algorithm achieves an accuracy of $2.77\times 10^{-6}$ in $7.76\times 10^{4}$ sec. The numbers of averaged samples for each level are $M_{0}=519634$, $M_{1}=6153$, and $M_{2}=243$. Obviously, the slow Monte Carlo convergence rate of $\mu=0.5$ is prohibitive for higher tolerances here.

3.2 Uncertain geometry: $d=2$, $N=16$

In our second example, we again consider a two-dimensional Poisson problem, but now with geometry of the computational domain being uncertain. We will refer to this example as the two hole test problem. The uncertain domain is defined by

[TABLE]

where the holes $D_{1}(y)=\overline{P_{1}P_{2}P_{3}P_{4}}$ and $D_{2}(y)=\overline{P_{5}P_{6}P_{7}P_{8}}$ are taken as quadrilaterals with the uncertain vertices

[TABLE]

Here, the random vector $y=(y_{1},\ldots,y_{16})^{T}$ consists of sixteen uniformly distributed random variables $y_{i}\sim{\cal U}[-1,1]$, $i=1,\ldots,16$. We set

[TABLE]

to represent different strength of uncertainty and hence anisotropy in the stochastic space. We impose a constant volume force $f\!\equiv\!1$ and fix the component by homogeneous Dirichlet boundary conditions on the whole boundary $\partial D(y)$, including $\partial D_{1}(y)\cup\partial D_{2}(y)$ with stochastically varying positions in space. Our goal is then to study the effect of this uncertainty on the expectation of the overall displacement calculated by

[TABLE]

This allows an assessment of the desired averaged load capacity of the component, which takes into account uncertainties in the manufacturing process.

Applying a parameter-dependent map, each domain $D(y)$ can be mapped to the fixed nominal domain $D_{0}=D(0)$ with $0\in\mathbb{R}^{16}$. Such a domain mapping approach was introduced by Xiu & Tartakovsky [38] and allows us to reformulate the problem in the form (1) with parameter-dependent coefficients on $D_{0}$. The well established theory for elliptic partial differential equations with random input data can then be applied without modifications to show the well-posedness of the setting with random domains.

All calculations start with an initial criss-cross structured mesh consisting of $1920$ triangles, which are adjusted to the random holes. In Fig. 5, the numerical solution and its corresponding adaptive mesh for a spatial tolerance ${\eta_{X}}=10^{-1}$ and the random vector $y\!=\!(0.5,0.5,-0.5,-0.5,-0.5,-0.5,1,-1,-1,-1,1,1,-1,1,-1,-1)^{T}$ are shown. The adaptive algorithm refines the mesh at the eight reentrant corners due to the fact that the exact solution contains a loss of regularity there. Exemplarily, we study the convergence rates of adaptive and uniform refinements for the nominal domain $D_{0}$, where $u(0)$ is contained in $H^{5/3}(D_{0})$. While the correctly adapted grids still recover the optimal order $-1$ for the approximation of $\psi(u(0))$ in terms of CPU time, the use of a sequence of uniform meshes sees the order drop down to the theoretical value $-0.66$, as can be seen in the left part of Fig. 6. The DWR-estimators perform quite well and deliver accurate upper bounds of the error. So we again set $C_{X}=1$ in (30) and use $s^{*}\!=\!1$ for adaptive meshes and $s^{*}=1.6$ for uniform meshes in assumption (A1). Note that the latter choice takes into account that due to the larger interior angle at the random holes for some parameter values $y\not=0$, the regularity of the solution, and thus also the spatial convergence rate, is slightly lower.

In order to estimate the parameter $\mu^{*}$, we apply the anisotropic Smolyak algorithm and calculate samples at reasonable costs with tolerances ${\eta_{X}}={\eta_{Y}}=10^{-2}/2^{i},i=0,\ldots,3$. The value of $\mathbb{E}[\psi(u)]$ for $i=3$ is taken as the reference value, leading to an approximate convergence order of $\mu^{*}=3.4$. From Fig. 6, we detect that the errors are significantly smaller than the prescribed tolerances, as it was also the case in the first test example. Therefore, we again choose $C_{Y}=0.1$ and set $\mu^{*}=3.4$ in assumption (A2). For later comparison, we compute a reference value $\mathbb{E}[\psi(u)]=4.758057\ldots$ with ${\eta_{X}}=10^{-4}$ and ${\eta_{Y}}=5\times 10^{-4}$. The required number of collocation points is $233$ and the vector of maximum polynomial degree is $(3\,2\,3\,3\,3\,3\,3\,3\,3\,3\,2\,2\,3\,2\,2\,2)^{T}$.

For the adaptive multilevel approach, we use again three levels and set the reduction factor to $q=0.2$. The spatial tolerances are then computed from ${\eta_{X_{k}}}=\epsilon q^{k-2}/2$, $k=0,1,2$, with the six overall accuracy requirements $\epsilon=10^{-2},5\times 10^{-3},2.5\times 10^{-3},10^{-3},5\times 10^{-4},2.5\times 10^{-4}$. Following the algorithm given in Table 1, we first compute ${\eta_{X_{-1}}}=\mathbb{E}[\psi(u_{0})]$ with ${\eta_{X}}={\eta_{Y}}={\eta_{X_{0}}}$ and then derive the stochastic tolerances ${\eta_{Y_{k}}}$ from (22) with $C_{Y}=0.1$, $\mu^{*}=3.4$, and $s^{*}=1$. The corresponding values are given in Table 2. For the multilevel approach with uniform spatial meshes, we use $s^{*}=1.6$ to determine the stochastic tolerances.

In Fig. 7, errors for the expected value $\mathbb{E}[\psi_{2}^{(\mathrm{ML})}]$ versus computing time are plotted for the three-level approach with adaptive and uniform spatial meshes. For comparison, we also show results for the single-level approach with adaptive spatial meshes, computed with ${\eta_{X}}=\epsilon/2$ and ${\eta_{Y}}=\epsilon/(2C_{Y})$. In all cases, except the first two for the single-level method, the overall tolerances are satisfied. The schemes work remarkably reliably. However, the achieved accuracy is often significantly better than the prescribed one – never more than a factor $10$ though. One possible reason could be that cancellation effects are typically overlooked if the overall error is split into two parts, which are then individually controlled.

The multilevel method with adaptive spatial meshes outperforms the single-level approach. The number of collocation points taken by the anisotropic Smolyak algorithm are listed in Table 2. We observe that the number of samples for the differences keeps constant, showing that the increasing samples in the zeroth level always catch enough information to eventually reach the tolerance. The corresponding numbers of collocation points for the single-level method are $(1,1,33,41,51,105)$. Note that the computing time also includes the effort for the estimation process, which is not visible in these numbers; see also the discussion for our realization of the adaptive anisotropic Smolyak algorithm above. It is also clear that the averaged constant $C_{Y}=0.1$ is too optimistic for the first runs with $\epsilon=10^{-2},5\times 10^{-3}$, which leads to the algorithm taking only one collocation point. The approximate orders of convergence for both methods, $p_{\mathrm{ML}}=-1/s^{*}=-1$ and $p_{\mathrm{SL}}=-1/(s^{*}+1/\mu^{*})\approx-0.77$, predict the observed asymptotic rates for the computing times quite well. However, the actual estimates can only serve as rough indicators for the achieved accuracy.

The multilevel method with uniform spatial meshes performs better than the single-level one for the first two tolerances, but becomes quickly inefficient for higher tolerances. Due to the larger value $s^{*}=1.6$, it needs significantly more samples for coarser meshes: $M_{0}=(65,131,567)$ and $M_{1}=(33,33,105)$, compared to the multilevel approach with adaptive spatial meshes, see Table 2. $M_{2}$ remains the same. The observed convergence order $-0.54$ is close to the predicted value $-1/s^{*}=-0.625$. Also for this example, it becomes obvious that uniform meshes cannot compete with adaptive meshes for higher tolerances.

This conclusion is also valid for tests with adaptive multilevel Monte Carlo. Results for tolerances $\epsilon=10^{-2},5\times 10^{-3}$ are shown in Fig. 7 for three levels. The numbers of optimized samples are $M_{0}=(4140,17013)$, $M_{1}=(48,108)$, $M_{2}=(5,6)$, respectively. All values are calculated from an average over 5 independent realizations. Although the variance reduction is quite high, leading to surprisingly small numbers $M_{2}$, the fast increasing numbers for $M_{0}$ are still too challenging, especially for higher tolerances.

4 Concluding remarks

Adaptive methods have the potential to drastically reduce the number of degrees of freedom and to realise optimal complexity in terms of computing time, even in cases where the exact solutions have spatial singularities. As well as providing a posteriori error estimates, adaptive methods usually outperform methods based on uniform mesh refinement. When coupled with multilevel stochastic algorithms, they can increase the computational efficiency in the sampling process when the stochastic dimension increases and hence provide a general tool to further delay the curse of dimensionality.

State-of-the-art adaptive methods implemented in the open-source Matlab packages p1afem and Sparse Grid Kit have been used in off-the-shelf fashion in this study. The numerical results for the adaptive multilevel collocation method demonstrate the reliability of the error control and a significant decrease in complexity compared to uniform spatial refinement methods and single-level stochastic sampling methods. While these advantages are already known for the methods involved, this study shows that the adaptive components can be combined to give an efficient algorithm for solving PDEs with random data to arbitrary levels of accuracy.

Appendix.

The pathwise solution of the one peak problem is given by

[TABLE]

with $\alpha(y_{1})=18y_{1}+11/2$. The quantity of interest is given by

[TABLE]

The specific choice $\beta=50$ taken in the numerical experiments leads to two simplifications.

• (S1)

The Dirichlet boundary condition ($u$ satisfying (A.1) on $\partial D$) may be replaced without significant loss of accuracy by the numerical approximation

[TABLE]

Justification: By direct computation, the maximum value of $u$ on the boundary occurs for $\alpha=1$, when the standard deviation $\sigma$ of the Gaussian peak is minimised (thus $\sigma_{*}=1/10$ and $y=(-1/4,y_{2})$ for all values $-1/4\leq y_{2}\leq 1/4$). For all such values of $y_{2}$ the nearest point on the boundary is ${{x}}=(-1,y_{2})$, and the shortest distance to the boundary is $d=3/4$ so that the maximum value of $u$ on the boundary is given by u_{\max}=\exp\big{(}-\frac{9\beta}{16}\big{)}\approx 6\cdot 10^{-13}.

• (S2)

Thus, we may readily compute a reference value (accurate to more than 10 digits):

[TABLE]

Justification: the integral in (A.2) is separable, thus

[TABLE]

The key point is that the inner integrals in the above expression are Gaussians with small variance. In both cases the integrand is close to zero at the end points $-1$ and $1$. Thus, in both cases the integral can be extended by zero to the range $(-\infty,\infty)$ (this could be made rigorous using probabilistic arguments) to give the following reference value approximation

[TABLE]

Acknowledgements. The first author is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the collaborative research center TRR154 “Mathematical modeling, simulation and optimisation using the example of gas networks” (Project-ID 239904186, TRR154/2-2018, TP B01), the Graduate School of Excellence Computational Engineering (DFG GSC233), and the Graduate School of Excellence Energy Science and Engineering (DFG GSC1070). The second author was supported by the UK Engineering and Physical Sciences Research Council (EPSRC) under grant number EP/K031368/1. The third author was funded by a Romberg visiting scholarship at the University of Heidelberg in 2019. The authors would also like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme “Uncertainty quantification for complex systems: theory and applications” (Jan–Jul 2018), where the work on this paper was initiated.