A domain mapping approach for elliptic equations posed on random bulk and surface domains

TL;DR

This paper develops a domain mapping method to efficiently approximate statistical moments of solutions to elliptic PDEs on random geometries, including surfaces and bulk-surface systems, with proven error estimates and numerical validation.

Contribution

It introduces a geometric analysis framework for reformulating elliptic equations on random surfaces onto a fixed domain, enabling finite element and Monte Carlo methods with optimal error bounds.

Findings

Theoretical convergence rates are established.

Numerical experiments confirm the accuracy of the method.

Applicable to both surface and bulk-surface elliptic problems.

Abstract

In this article, we analyse the domain mapping method approach to approximate statistical moments of solutions to linear elliptic partial differential equations posed over random geometries including smooth surfaces and bulk-surface systems. In particular, we present the necessary geometric analysis required by the domain mapping method to reformulate elliptic equations on random surfaces onto a fix deterministic surface using a prescribed stochastic parametrisation of the random domain. An abstract analysis of a finite element discretisation coupled with a Monte-Carlo sampling is presented for the resulting elliptic equations with random coefficients posed over the fixed curved reference domain and optimal error estimates are derived. The results from the abstract framework are applied to a model elliptic problem on a random surface and a coupled elliptic bulk-surface system and the theoretical convergence rates are confirmed by numerical experiments.

Figures

Tables

Table 1. Table 1 . Error in L 2 ​ ( Ω M ; L 2 ​ ( Γ 0 ) ) superscript 𝐿 2 superscript Ω 𝑀 superscript 𝐿 2 subscript Γ 0 L^{2}(\Omega^{M};L^{2}(\Gamma_{0})) .

$h$	$M$	$E_{L^2\left({\Gamma }_0\right)}$	$eoc\left(h\right)$	$eoc\left(M\right)$
0.171499	1	0.776832	-	-
0.0877058	16	0.387486	1.03722	-0.250864
0.0441081	256	0.106022	1.88556	-0.467444
0.0220863	4096	0.0267303	1.99202	-0.496955

Table 2. Table 2 . Error in L 2 ​ ( Ω M ; H 1 ​ ( Γ 0 ) ) superscript 𝐿 2 superscript Ω 𝑀 superscript 𝐻 1 subscript Γ 0 L^{2}(\Omega^{M};H^{1}(\Gamma_{0})) .

$h$	$M$	$E_{H^1\left({\Gamma }_0\right)}$	$eoc\left(h\right)$	$eoc\left(M\right)$
0.171499	64	4.89172	-	-
0.0877058	256	3.68809	0.421176	-0.203734
0.0441081	1024	1.90402	0.961875	-0.476911
0.0220863	4096	0.961782	0.987348	-0.492633

Table 3. Table 3 . Errors in L 2 ​ ( Ω M ; L 2 ​ ( D 0 ) ) superscript 𝐿 2 superscript Ω 𝑀 superscript 𝐿 2 subscript 𝐷 0 L^{2}(\Omega^{M};L^{2}(D_{0})) and L 2 ​ ( Ω M ; L 2 ​ ( Γ 0 ) ) . superscript 𝐿 2 superscript Ω 𝑀 superscript 𝐿 2 subscript Γ 0 L^{2}(\Omega^{M};L^{2}(\Gamma_{0})).

$h$	$M$	Bulk error $E_{L^2\left(D_0\right)}$	$eoc\left(h\right)$	$eoc\left(M\right)$	Surface error $E_{L^2\left({\Gamma }_0\right)}$	$eoc\left(h\right)$	$eoc\left(M\right)$
0.27735	1	0.619144	-	-	5.0787	-	-
0.156174	16	0.198298	1.98249	-0.410651	1.06707	2.71654	-0.562702
0.0830455	256	0.0540441	2.05828	-0.468866	0.28356	2.0983	-0.477981
0.0428353	4096	0.0152612	1.91003	-0.456067	0.0723061	2.06414	-0.492866

Table 4. Table 4 . Errors in L 2 ​ ( Ω M ; H 1 ​ ( D 0 ) ) superscript 𝐿 2 superscript Ω 𝑀 superscript 𝐻 1 subscript 𝐷 0 L^{2}(\Omega^{M};H^{1}(D_{0})) and L 2 ​ ( Ω M ; H 1 ​ ( Γ 0 ) ) superscript 𝐿 2 superscript Ω 𝑀 superscript 𝐻 1 subscript Γ 0 L^{2}(\Omega^{M};H^{1}(\Gamma_{0})) .

$h$	$M$	Bulk error $E_{L^2\left(D_0\right)}$	$eoc\left(h\right)$	$eoc\left(M\right)$	Surface error $E_{L^2\left({\Gamma }_0\right)}$	$eoc\left(h\right)$	$eoc\left(M\right)$
0.27735	64	3.41133	-	-	15.5792	-	-
0.156174	256	2.17523	0.783494	-0.324584	7.85391	1.1926	-0.494068
0.0830455	1024	1.08874	1.09584	-0.499252	4.20041	0.990894	-0.451441
0.0428353	4096	0.55599	1.01511	-0.484767	2.12783	1.02727	-0.490574

Equations

$$\text{BODY}$$

$$\Vert f{\Vert }_{L^p(\Omega ;V)}=\{\begin{array}{ll}{\left({\int }_{\Omega }\Vert f(\omega ){\Vert }_V^{p}d\mathbb{P}(\omega )\right)}^{\frac{1}{p}}& p\in [1,\infty )\\ \underset{\omega }{\mathrm{esssup}}\Vert f(\omega ){\Vert }_V& p=\infty ,\end{array}$$

$$L^2(\Omega ;V)\cong L^2(\Omega )\otimes V.$$

$$-\Delta u(\omega )$$

$$u(\omega )$$

$$-\frac{1}{\sqrt{g(\omega )}}\nabla \cdot \left(\sqrt{g(\omega )}{G}^{-1}(\omega )\nabla (u\circ \varphi )(\omega )\right)$$

$$u(\omega )$$

$$G(\omega )=\nabla {\varphi }^{\top }(\omega )\nabla \varphi (\omega )g(\omega )=\text{det }G(\omega ).$$

$$\varphi (\omega ,x)=\mathbb{E}[\varphi ](x)+\sum _{k=1}^N{Y}_k(\omega ){\varphi }_k(x)$$

$$D(\omega )=\{\varphi (\omega ,x)|x\in D_0\}\subset {\mathbb{R}}^n,$$

$$\mathbb{E}[D]=\{\mathbb{E}[\varphi ](x)|x\in D_0\}.$$

$${\mathbb{E}}_V[D]=\{x|p(x)\ge \lambda \},$$

$$\mathbb{E}[\varphi ]:D_0\to \mathbb{E}[D],$$

$${\varphi }_e=\varphi \circ \mathbb{E}[\varphi {]}^{-1}.$$

$$\text{QoI}[u]=\mathbb{E}[u\circ {\varphi }_e]\text{on }\mathbb{E}[D].$$

$$U\cap \Gamma =\{x\in U|\phi (x)=0\}.$$

$${\nu }^{\Gamma }=\pm \frac{\nabla \phi }{|\nabla \phi |}$$

$${\nabla }_{\Gamma }f=\nabla \bar{f}-(\nabla \bar{f}\cdot {\nu }^{\Gamma }){\nu }^{\Gamma }={\mathcal{P}}_{\Gamma }\nabla \bar{f}$$

$${\nabla }_{\Gamma }f=({\underline{D}}_1^{\Gamma }f,...,{\underline{D}}_{n+1}^{\Gamma }f{)}^{\top }.$$

$${\Delta }_{\Gamma }f={\nabla }_{\Gamma }\cdot {\nabla }_{\Gamma }f=\sum _{i=1}^{n+1}{\underline{D}}_i^{\Gamma }{\underline{D}}_i^{\Gamma }f.$$

$$x=a^{\Gamma }(x)+d^{\Gamma }(x){\nu }^{\Gamma }(a^{\Gamma }(x)).$$

$$\nabla d^{\Gamma }(x)$$

$$\nabla a^{\Gamma }(x)$$

$$\Gamma =\{\varphi (x)|x\in {\Gamma }_0\},$$

$$\Gamma =\{\varphi (x)=x+h(x){\nu }^{{\Gamma }_0}(x)|x\in {\Gamma }_0\}$$

$$D=\{\varphi (x)|x\in D_0\}\Gamma =\{\varphi (x)|x\in {\Gamma }_0\}$$

$$\hat{f}=f\circ \varphi .$$

$$X:U\to W\cap \Gamma$$

$${\nabla }_{\Gamma }f\circ X$$

$${\Delta }_{\Gamma }f\circ X$$

$$G_{{\Gamma }_0}={\nabla }_{{\Gamma }_0}{\varphi }^{\top }{\nabla }_{{\Gamma }_0}\varphi +{\nu }^{{\Gamma }_0}\otimes {\nu }^{{\Gamma }_0},$$

$$X=\varphi \circ \sigma :U\to W\cap \Gamma$$

Full text

TODO

1. Introduction

In the mathematical characterization of numerous scientific and engineering systems, the topology of the domain may not be precisely described. The main sources of uncertainty are usually insufficient data, measurement errors or manufacturing variability. This uncertainty in the geometry often naturally appears in many applications including surface imaging, manufacturing of nano-devices, material science and biological systems. As a result, the analysis of uncertainty in the computational domain has become an interesting and rich mathematical field.

A comprehensive summary concerning the first directions in the treatment of elliptic partial differential equations (PDEs) in random domains can be found in [32, 18, 26, 28, 4] and recently [10] for a parabolic equation on a randomly evolving domain. Aside from the fictitious domain method [4, 27, 28], the main approaches utilize a probabilistic framework by describing the random boundary of the domain with a random field. This probabiltistic approach is usually proceeded with one of two main techniques: the perturbation approach and the domain mapping method. The perturbation approach [17, 20] exploits a shape Taylor expansion with respect to the boundary random field to represent the solution, however as a result is limited to consideration of only small random deformations. The domain mapping approach [32, 19, 5] on the other hand allows does not suffer the same limitations. The key idea behind this method is to define an extension of the random boundary process into the interior domain to form a complete random mapping for the whole domain and then to use this domain mapping to transform the original partial differential equation on the random domain onto the fixed deterministic reference domain resulting in an partial differential equations with random coefficients. For the latter formulation, there is a wealth of literature available on numerical techniques to compute any quantities of interest, see for example [24, 23, 16]. The aim of this paper, is to incorporate the domain mapping method with the well-developed field of surface PDEs [14, 13, 8] which has so far only considered uncertainty in the coefficients of the considered PDEs, see [11]. This will lead to more realistic geometric description of many of the situations previously dicussed. Note that while the domain mapping method will be applicable to domains with random rough surfaces, we will only choose to focus on sufficiently smooth random surfaces and leave the rough case for future considerations.

The layout of the article is as follows. In Section 2, we provide an overview of the domain mapping method for partial differential equations in flat random domains and furthermore discuss suitable notions for the expectation of a family of random domains. In Section 3, we introduce the necessary geometric analysis and computations required to apply the domain mapping method to elliptic partial differential equations posed on random surfaces. In Section 4, we present a model elliptic problem on a random surface and a coupled elliptic system on a random bulk-surface, and analyse weak formulations in both the stochastic and spatial variables for the reformulated equations with stochastic coefficients on the fixed reference surface and bulk-surface respectively. Section 5 provides an abstract analysis of a finite element discretisation incorporating a pertubation to the variational set-up due to a first order approximation of the curved reference domain, and couples with a Monte-Carlo sampling to approximate the first moment of the solution. An optimal error estimate is derived and subsequently applied in Section 6, to two discretisations of the proposed reformulated problems. We conclude in Section 7, by presenting numerical results confirming the theoretical rate of convergence.

2. The domain mapping method

We begin with a brief introduction on spaces of random fields. For further details on these spaces, we refer the reader to [24]. Note throughout this paper, we will let $(\Omega,\mathcal{F},\mathbb{P})$ denote a complete, separable probability space consisting of a sample space $\Omega$, a $\sigma-$algebra of events $\mathcal{F}$ and a probability measure $\mathbb{P}$.

2.1. Random field notation

For a given Banach space $V$ and $p\in[1,\infty]$, the Lebesgue-Bochner space $L^{p}(\Omega;V)$ consists of all strongly $\mathcal{F}-$measurable functions $f:\Omega\rightarrow V$ for which the norm

[TABLE]

is finite. For convenience, we will express the parameters of a given random field $(f(\omega))(x)$ by $f(\omega,x)$. In the case that $V$ is a separable Hilbert space, it follows that $L^{2}(\Omega;V)$ is also a separable Hilbert space and furthermore is isomorphic to the tensor product

[TABLE]

For details, see [29].

2.2. The domain mapping method

To illustrate the key concepts of the domain mapping method, consider the following boundary value problem

[TABLE]

posed on an open, connected, bounded domain $D(\omega)\subset\mathbb{R}^{2}$ with a random boundary $\Gamma(\omega)=\partial D(\omega)$. Here the prescribed random field $f(\omega):D(\omega)\rightarrow\mathbb{R}$ and additionally the boundary, will be assumed to be sufficiently regular to ensure well-posedness for $a.e.$ $\omega$. The first essential feature of the domain mapping method is the representation of the stochastic boundary via a random field. More precisely, in the above context we will assume that there exists a random field $\phi\in L^{\infty}(\Omega;C^{0}(\Gamma_{0};\mathbb{R}^{2})),$ that maps a fixed closed curve $\Gamma_{0}\subset\mathbb{R}^{2}$ onto realisations of the random boundary $\phi(\omega,\cdot):\Gamma_{0}\rightarrow\Gamma(\omega),$ see figure 1. The next step in the method is to define an extension of the boundary process into the interior to form a stochastic mapping $\phi(\omega,\cdot):\overline{D_{0}}\rightarrow\overline{D(\omega)}$ for the whole domain. For instance, [32] proposed an extension based on the solution of the Laplace equation over the unit square with boundary conditions prescribed by segments of the random boundary. However, alternative approaches may wished to be considered depending on the application in question and the geometry of the computational reference domain.

With a complete domain mapping at hand, the random domain problem (2.2) can now be reformulated as a partial differential equation with random coefficients over the fixed deterministic domain $D_{0}$,

[TABLE]

where the specific random coefficients for this particular problem are given by

[TABLE]

We now have access to a wide breadth of numerical techniques, including Monte-Carlo [6, 21] and the stochastic Galerkin method [2, 25], to compute any statistical quantities of interest.

Remark 2.1**.**

Note that the choice of the reference domain $D_{0}$, for the stochastic domain mapping $\phi$ describing the complete random geometry in question, is arbitrary and should be chosen in such a way that simplifies the computation at hand. Furthermore in practice, only statistical properties such as the expectation and two-point covariance function of the stochastic mapping $\phi$ will be known. As a result, an approximation of the true process will instead be used, commonly taking the form of a truncated series

[TABLE]

with centered, uncorrelated random coefficients $Y_{k}$ with unit variance, such as a truncated Karhunen-Loève expansion. Considerations of the induced error is beyond the scope of this paper and we instead refer the reader to [18].

2.3. Expected domain and quantity of interest

In order to give a precise definition of our quantity of interest, which for our purpose shall be some notion of a mean solution, we will first need to fix a suitable domain of definition. A natural choice would be the parametrisation based expected domain, introduced in [7] for random star-shaped domains, which we shall generalise as follows.

Definition 2.1** (Parametrisation based expected domain).**

Given a family of random Lipschitz domains

[TABLE]

parametrised over a fixed Lipchitz domain $D_{0}\subset\mathbb{R}^{n}$ under the Lipschitz continuous mapping $\phi(\omega,\cdot):D_{0}\rightarrow\mathbb{R}^{n}$. Assuming $\phi(\cdot,x)$ is integrable for all $x\in D_{0}$, the parametrisation based expected domain $\mathbb{E}[D]$ of the random domain $D(\omega)$ is given by

[TABLE]

Remark 2.2**.**

Note that there are other alternative methods in which to define the expected value of a family of random sets. For example, we could characterise the random set $D(\omega)$ as an indicator function $1_{D(\omega)}$ and then use its average, the so-called coverage function $p(x)=\mathbb{P}(x\in D(\omega))$ to define the expected value to be set

[TABLE]

where the parameter $\lambda>0$ is selected in a such a way that the volume of $\mathbb{E}_{V}[D]$ is close as possible to the expected volume of the random sets $D(\omega)$. This is known as the Vorob’ev expectation and was shown in [7] not to coincide with the parameterisation based expectation. Although there is no canonical definition of the expected value of a random domain, the parametrisation based expected domain fits naturally in the setting of the domain mapping method and thus will be adopted.

Assumption 2.1**.**

We will assume that the expected value of the stochastic mapping

[TABLE]

is bi-Lipschitz continuous. This will ensure that the parametrisation based expected domain $\mathbb{E}[D]$ is also Lipschitz continuous and furthermore of the same dimension as $D_{0}$ and $D(\omega)$.

We will denote the induced zero-mean stochastic mapping between the parametrisation based expected domain $\mathbb{E}[D]$ and realisations of the random domain $D(\omega)$ by

[TABLE]

See figure 2 for an illustration of the different mappings and domains. Our quantity of interest can now be defined on the expected domain as follows.

Definition 2.2** (QoI).**

Given a random field $u(\omega,\cdot):D(\omega)\rightarrow\mathbb{R}$ defined over the family of random Lipschitz domains given in (2.3), the expected value of the random field is given by

[TABLE]

As previously discussed our aim is to apply the domain mapping method for random domains which involve random surfaces. We will therefore now proceed with some preliminary computations of geometric quantities as well as tangential derivatives of functions given over parametrised hypersurfaces in terms of quantites of the reference surface and derivatives of the domain mapping and corresponding pull-back function. This will provide a basis for the domain mapping method to be employed to several model PDEs over random surfaces.

3. Computations for the pull-back of tangential differential operators and geometric quanitities of parametrised hypersurfaces

Let us first introduce some notation for hypersurfaces that will be adopted throughout this paper. For a more detailed introduction, see [14].

3.1. Hypersurface notation

A set $\Gamma\subset\mathbb{R}^{n+1}$ is said to be a $C^{k}$-hypersurface for $k\in\mathbb{N}\cup\{\infty\}$, provided that for every $x\in\Gamma$ there exists an open set $U\subset\mathbb{R}^{n+1}$ containing $x$ and a smooth function $\varphi\in C^{k}(U)$ such that $\nabla\varphi(x)\neq 0$ on $U\cap\Gamma$ and

[TABLE]

The unit normal vector field $\nu^{\Gamma}$ to the hypersurface $\Gamma$ can be computed via

[TABLE]

with a choice of orientation. For a differentiable function $f:\Gamma\rightarrow\mathbb{R}$, we define the tangential gradient by

[TABLE]

where $\mathcal{P}_{\Gamma}=I-\nu^{\Gamma}\otimes\nu^{\Gamma}$ is the projection operator mapping onto the tangent space $T\Gamma$ to the hypersurface $\Gamma$ and $\bar{f}$ is a smooth extension of $f$ to an open neighbourhood in $\mathbb{R}^{n+1}$. It can be shown that the tangential gradient is independent of the extension chosen [14, Lemma 2.4] and we shall denote its components by

[TABLE]

We shall further denote the tangential derivative of the unit normal by $\mathcal{H}^{\Gamma}=\nabla_{\Gamma}\nu^{\Gamma}$ and will refer to this matrix as the extended Weingarten map. It can be shown that $\mathcal{H}^{\Gamma}$ is symmetric with a zero eigenvalue corresponding to the unit normal vector $\nu^{\Gamma}$ and furthermore agrees with the Weingarten map when restricted to the tangent space $T\Gamma$, see [8] for details. The Laplace-Beltrami operator is then defined for twice differentiable functions as follows

[TABLE]

We next introduce the Fermi coordinates with the following well-known lemma [14, Lemma 2.8]. These are a global coordinate system defined in an open neighbourhood around $\Gamma$ in which every point can be uniquely expressed in terms of its signed distance $d^{\Gamma}(x)$ and its closest point $a^{\Gamma}(x)$ on the surface $\Gamma$ .

Lemma 3.1**.**

Let $d^{\Gamma}$ denote the signed distance function to $\Gamma$ oriented in the chosen direction of the unit normal vector field $\nu^{\Gamma}$. Then there exists $\delta>0$ such that for every $x\in U_{\delta}:=\{y\in\mathbb{R}^{n+1}\,|\,|d^{\Gamma}(y)|<\delta\}$ there exists a unique point $a^{\Gamma}(x)\in\Gamma$ that satisfies

[TABLE]

Furthermore, assuming $\Gamma\in C^{2}$ it follows that $d^{\Gamma}\in C^{2}(U_{\delta})$ and $a^{\Gamma}\in C^{1}(U_{\delta})$ with

[TABLE]

3.2. Geometric settings

As a point of reference, we will now describe the deterministic geometric settings that will be considered for the parametrised surfaces in the subsequent calculations. In each case, the reference surface $\Gamma_{0}\subset\mathbb{R}^{n+1}$ will be assumed to be of class at least $C^{2}$ and oriented by the unit normal vector field $\nu^{\Gamma_{0}}$. The general geometric setting for the parametrised surface will be as follows.

Geometric setting 1** (Parametrised surface).**

The hypersurface $\Gamma\subset\mathbb{R}^{n+1}$ will be given by

[TABLE]

for a given mapping $\phi:\Gamma_{0}\rightarrow\mathbb{R}^{n+1}$.

We will further consider the special case, where the parametrised surface has the following graph-like representation over the reference surface.

Geometric setting 2** (Graph-like surface).**

The surface $\Gamma\subset\mathbb{R}^{n+1}$ will be prescribed by

[TABLE]

for a given height function $h:\Gamma_{0}\rightarrow\mathbb{R}$ defined over the reference surface.

Additionally, we will consider the case where the surface is compact (and thus without a boundary) and encloses an open bulk domain.

Geometric setting 3** (Parametrised bulk-surface).**

The open bulk domain $D\subset\mathbb{R}^{n+1}$ and its boundary $\Gamma=\partial D$ which is a surface, will be given by

[TABLE]

for a given parametrisation $\phi:\overline{D_{0}}\rightarrow\mathbb{R}^{n+1}$ defined over an open bulk domain $D_{0}\subset\mathbb{R}^{n+1}$ with boundary $\Gamma_{0}=\partial D_{0}$.

Note that in each case, the given parametrisation $\phi$ will be assumed to be a sufficiently smooth diffeomorphism for the calculation in question. Furthermore, we shall denote the associated pull-back of a given function $f$ defined over the parametrised domain onto the reference domain by

[TABLE]

3.3. The tangential gradient and Laplace-Beltrami operator

Considering a general parametrised hypersurface $\Gamma$ as described in (3.8), we will now compute expressions for the pull-back of the tangential gradient $\nabla_{\Gamma}$ and Laplace-Beltrami operator $\Delta_{\Gamma}$ onto the reference surface $\Gamma_{0}$ under the domain mapping $\phi$. As a motivation for these calculations, let us first recall that for a given local parametrisation

[TABLE]

of the hypersurface $\Gamma$, where $U\subset\mathbb{R}^{n}$ and $W\subset\mathbb{R}^{n+1}$ denote open sets, we can express the tangential gradient and Laplace-Beltrami operator in local coordinates as follows

[TABLE]

where $F=f\circ X$ and the first fundamental form $G:U\rightarrow\mathbb{R}^{n\times n}$ is defined as $G=\nabla X^{\top}\nabla X$ with $g=\text{det}G$. In deriving expressions for the pull-back onto the reference surface $\Gamma_{0}$ instead of the local coordinates, we will see similar expresssions to (3.12) and (3.13) but with the first fundamental form replaced by the following tensor $G_{\Gamma_{0}}:\Gamma_{0}\rightarrow\mathbb{R}^{(n+1)\times(n+1)}$ defined by

[TABLE]

where we will similarly denote its determinant by $g_{\Gamma_{0}}=\text{det}G_{\Gamma_{0}}$. This tensor can be seen to arise by considering a local parametrisation $\sigma:U\rightarrow V\cap\Gamma_{0}$ of the reference surface $\Gamma_{0}$, with $V\subset\mathbb{R}^{n+1}$ denoting an open set, and the induced local parametrisation

[TABLE]

of the hypersurface $\Gamma$. By computing the first fundamental form with the chain rule, we observe that

[TABLE]

Since $\nabla_{\Gamma_{0}}\phi^{\top}\nabla_{\Gamma_{0}}\phi\nu^{\Gamma_{0}}=0$ and its restriction to the tangent space maps $\nabla_{\Gamma_{0}}\phi^{\top}\nabla_{\Gamma_{0}}\phi:T\Gamma_{0}\rightarrow T\Gamma_{0}$, we are able to extend in the normal direction as in (3.14) to form an invertible matrix. Furthermore as $\nabla\sigma\in T\Gamma_{0}$, it follows that we have

[TABLE]

Note that the given extension (3.14) in the normal direction is a natural choice since the surface measures $dA_{\Gamma}$ and $dA_{\Gamma_{0}}$ of the respective surfaces can be shown to satisfy the relation $dA_{\Gamma}=\sqrt{g_{\Gamma_{0}}}dA_{\Gamma_{0}}$ under the domain transformation mapping $\phi$. We now continue by proving that a similar expression to (3.12) holds for the pull-back of the tangential gradient.

Lemma 3.2** (Tangential gradient).**

Given any differentiable function $f:\Gamma\rightarrow\mathbb{R}$, the pull-back of the tangential gradient onto the reference surface $\Gamma_{0}$ is given by

[TABLE]

Proof.

Differentiating the associated pull-back function $\hat{f}=f\circ\phi$ and applying the chain rule for tangential derivatives gives

[TABLE]

Since the tangential gradient of the surface parametrisation bijectively maps $\nabla_{\Gamma_{0}}\phi:T_{(\cdot)}\Gamma_{0}\rightarrow T_{\phi(\cdot)}\Gamma$ and additionally has kernel equal to span$\{\nu^{\Gamma_{0}}\}$, we see that in order to invert the matrix $\nabla_{\Gamma_{0}}\phi$, we must first modify the corresponding linear map to bijectively map the space span$\{\nu^{\Gamma_{0}}\}$ into span$\{\nu^{\Gamma}\circ\phi\}$. One possible solution is to add the linear map $L:\mathbb{R}^{n+1}\rightarrow\mathbb{R}^{n+1}$ characterised by

[TABLE]

which translates to adding the following tensor product

[TABLE]

and thus leads to (3.15). For the second equality, we again use the property that the restriction $\nabla_{\Gamma_{0}}\phi:T_{(\cdot)}\Gamma_{0}\rightarrow T_{\phi(\cdot)}\Gamma$ is a bijective mapping to express $(\nabla_{\Gamma}f)\circ\phi=\nabla_{\Gamma_{0}}\phi\,\alpha$ for some $\alpha\in T\Gamma_{0}$. Substituting into (3.16) then gives

[TABLE]

Hence we deduce $\alpha=G_{\Gamma_{0}}^{-1}\nabla_{\Gamma_{0}}\hat{f}$ and obtain the second equality. ∎

Remark 3.1**.**

Note that the chain rule for tangential gradients (3.15) holds for any choice of orientation of the unit normals $\nu^{\Gamma_{0}}$ and $\nu^{\Gamma}$ as a result of (3.17).

Let us denote the given extension of the tangential gradient of the surface parametrisation appearing in (3.15) by $B=\left(b_{ij}\right)_{i,j}$,

[TABLE]

and furthermore denote its determinant by $b=\det B$ and the entries of its inverse by $B^{-1}=\left(b^{ij}\right)_{i,j}$. We observe with the orthogonality result $\nabla_{\Gamma_{0}}\phi^{\top}(\nu^{\Gamma}\circ\phi)=0$ which follows from the property that the restriction maps $\nabla_{\Gamma_{0}}\phi:T_{(\cdot)}\Gamma_{0}\rightarrow T_{\phi(\cdot)}\Gamma$, that

[TABLE]

Consequently, we have

[TABLE]

We can now compute the pull-back of the Laplace-Beltrami operator onto the reference surface as follows.

Lemma 3.3** (Laplace-Beltrami operator).**

Given any $f:\Gamma\rightarrow\mathbb{R}$ twice differentiable, the pull-back of the Laplace-Beltrami operator is given by

[TABLE]

Proof.

By the chain rule for tangential gradients (3.15), we can express the Laplace-Beltrami operator as

[TABLE]

The last step follows from the observations (3.19) and (3.20). We continue by proving that the remaining terms vanish. Recalling Jacobi’s formula for the derivative of a determinant $\underline{D}_{j}^{\Gamma_{0}}\text{det}B=\text{det}B\,\text{trace}\left(B^{-1}\underline{D}_{j}^{\Gamma_{0}}B\right)$ and computing the derivative of the inverse matrix $\underline{D}_{j}^{\Gamma_{0}}B^{-1}=-B^{-1}\underline{D}_{j}^{\Gamma_{0}}BB^{-1}$ gives

[TABLE]

It therefore follows after relabelling indices that

[TABLE]

Differentiating $b_{mj}:=\underline{D}_{j}^{\Gamma_{0}}\phi_{m}+(\nu_{m}^{\Gamma}\circ\phi)\nu_{j}^{\Gamma_{0}}$ yields

[TABLE]

By the symmetry of the Weingarten map $\underline{D}_{l}^{\Gamma_{0}}\nu_{j}^{\Gamma_{0}}=\underline{D}_{j}^{\Gamma_{0}}\nu_{l}^{\Gamma_{0}}$, we see that the last two terms cancel. We next interchange tangential derivatives [14, Lemma 2.6]

[TABLE]

to obtain

[TABLE]

Substituting into (3.22), we arrive at the following expression for the remaining terms

[TABLE]

Examining the first term, we have

[TABLE]

Since $G_{\Gamma_{0}}=\nabla_{\Gamma_{0}}\phi^{\top}\nabla_{\Gamma_{0}}\phi+\nu^{\Gamma_{0}}\otimes\nu^{\Gamma_{0}}$ and thus $G_{\Gamma_{0}}^{-1}\nu^{\Gamma_{0}}=\nu^{\Gamma_{0}}$, the first term vanishes. For the second and third term, we observe that

[TABLE]

Therefore as a consequence of the orthogonality results $\nabla_{\Gamma_{0}}\phi^{\top}(\nu^{\Gamma}\circ\phi)=0$ and $\nabla_{\Gamma_{0}}(\nu^{\Gamma}\circ\phi)^{\top}(\nu^{\Gamma}\circ\phi)=0$ which can be seen by

[TABLE]

we conclude $\textup{I}+\textup{II}=0$. ∎

We next compute the specific form of the coefficients appearing in the pull-back of the tangential gradient (3.15) and the Laplace-Beltrami operator (3.21), for the particular case of a graph-like parametrisation over the reference surface.

Lemma 3.4** (Graphical case).**

*Assuming that the parametrisation of the hypersurface $\Gamma$ has the particular graph-like representation described in (3.9) for a given height function $h:\Gamma_{0}\rightarrow\mathbb{R}$, then the inverse and determinant of the tensor $G_{\Gamma_{0}}$ defined in (3.14) simplify to give *

[TABLE]

Here $A:=\left(I+h\mathcal{H}^{\Gamma_{0}}\right)^{-1}$ and $\{\kappa_{j}^{\Gamma_{0}}\}_{j}$ denotes the non-zero eigenvalues of the extended Weingarten map $\mathcal{H}^{\Gamma_{0}}$.

Proof.

Differentiating the given surface parametrisation $\phi(x)=x+h(x)\nu^{\Gamma_{0}}(x)$, we obtain

[TABLE]

Expanding the tensor $G_{\Gamma_{0}}=\nabla_{\Gamma_{0}}\phi^{\top}\nabla_{\Gamma_{0}}\phi+\nu^{\Gamma_{0}}\otimes\nu^{\Gamma_{0}}$ and cancelling the orthogonal terms with the tensor product identity $(a\otimes b)(c\otimes d)=(b\cdot c)a\otimes d$, yields

[TABLE]

Taking the inverse with the identity $\left(I+a\otimes b\right)^{-1}=I-\frac{a\otimes b}{1+a\cdot b}$ we obtain (3.23). For (3.24), we take the determinant and apply $\det(I+a\otimes b)=1+a\cdot b$, which leads to

[TABLE]

Since $A^{-1}=I+h\mathcal{H}^{\Gamma_{0}}$ has eigenvalues $1$ and $\{1+h\kappa_{j}^{\Gamma_{0}}\}_{j=1}^{n}$, we deduce $\det(A^{-1})=\prod_{j=1}^{n}\left(1+h\kappa_{j}^{\Gamma_{0}}\right)$ and thus obtain the stated result for $\sqrt{g_{\Gamma_{0}}}=\sqrt{detG_{\Gamma_{0}}}$. ∎

3.4. The unit normal and extended Weingarten map

We continue by computing expressions for the pull-back onto the reference surface $\Gamma_{0}$, of the unit normal $\nu^{\Gamma}$ and extended Weingarten map $\mathcal{H}^{\Gamma}$ for a general parametrised hypersurface $\Gamma$ as given in (3.8). To obtain an expression for the unit normal, we smoothly extend the given surface parametrisation $\phi:\Gamma_{0}\rightarrow\Gamma$ to a $C^{1}-$diffeomorphic mapping $\bar{\phi}:U\rightarrow V$ between some open sets $U$ and $V$ containing $\Gamma_{0}$ and $\Gamma$ respectively. The existence of such a mapping is gauranteed by the Whitney extension theorem [31]. We now have a level-set description of $\Gamma$

[TABLE]

consequently leading to the following expression for the unit normal vector field due to (3.2).

Lemma 3.5** (Unit normal).**

The pull-back of the unit normal vector field $\nu^{\Gamma}$ of the parametrised surface $\Gamma$ described in (3.8) onto the reference surface $\Gamma_{0}$, is given by

[TABLE]

Note that (3.25) can be shown to be independent of the extension chosen. As an example of a possible extension of the given surface parametrisation, we now consider the case of a graph-like surface.

Corollary 3.1** (Graphical case).**

Assuming that the hypersurface $\Gamma$ has the particular graph-like form described in (3.9), then the pull-back of the unit normal vector field $\nu^{\Gamma}$ is given by

[TABLE]

Here the orientation has been chosen to coincide with the reference surface $\Gamma_{0}$ when the height function is identically zero. Recall that $A:=(I+h\mathcal{H}^{\Gamma_{0}})^{-1}$.

Proof.

We extend the given surface parametrisation $\phi:\Gamma_{0}\rightarrow\Gamma$ defined by

[TABLE]

to a thin tubular neighbourhood $U=\{x\in\mathbb{R}^{n+1}\,|\,|d^{\Gamma_{0}}(x)|<\delta\}$ around $\Gamma_{0}$ of width $\delta>0$ as follows

[TABLE]

For $\delta>0$ sufficiently small, its image $V=\bar{\phi}(U)$ is contained within the neighbourhood in which the Fermi coordinates $(a^{\Gamma_{0}}(x),d^{\Gamma_{0}}(x))$ are well defined. Consequently, the extension $\bar{\phi}:U\rightarrow V$ which equivalently acts upon the Fermi coordinates as follows

[TABLE]

can be seen to be a bijective mapping. Computing its derivative and evaluating on the reference surface $\Gamma_{0}$, recalling that $\nabla d^{\Gamma_{0}}=\nu^{\Gamma_{0}}$ and $\nabla a^{\Gamma_{0}}=\mathcal{P}_{\Gamma_{0}}$ on $\Gamma_{0}$ by (3.6) and (3.7), we obtain

[TABLE]

Hence taking the inverse with the identity $(I+a\otimes b)^{-1}=1-\frac{a\otimes b}{1+a\cdot b}$ and recalling that $A\nu^{\Gamma_{0}}=\nu^{\Gamma_{0}}$, we deduce

[TABLE]

and thus obtain the stated result. Note that $\nu^{\Gamma_{0}}-A\nabla_{\Gamma_{0}}h\neq 0$ since the matrix $A=(I+h\mathcal{H}^{\Gamma_{0}})^{-1}$ maps $A:T\Gamma_{0}\rightarrow T\Gamma_{0}$. ∎

We next compute the pull-back of the extended Weingarten map $\mathcal{H}^{\Gamma}$ for a general parametrised surface $\Gamma$. Since the restriction of the derivative of the surface parametrisation maps $\nabla_{\Gamma_{0}}\phi(\cdot):T_{(\cdot)}\Gamma_{0}\rightarrow T_{\phi(\cdot)}\Gamma$, we consequently have

[TABLE]

for all $j=1,...,n+1$. Differentiating, we obtain

[TABLE]

Next, we define $L:\Gamma_{0}\rightarrow\mathbb{R}^{(n+1)\times(n+1)}$ by

[TABLE]

and rewrite (3.27) as

[TABLE]

It therefore follows from an application of the chain rule and the symmetry of the extended Weingarten map that

[TABLE]

We can then extend the tangential derivative $\nabla_{\Gamma_{0}}\phi$ to an invertible matrix as previously discussed in (6.17), to obtain the following result.

Lemma 3.6** (Extended Weingarten map).**

Let the orientation of the parametrised hypersurface $\Gamma$ described in (3.8), be fixed by a choice of a unit normal vector field $\nu^{\Gamma}$. Then the pull-back of the extended Weingarten map is given by

[TABLE]

Note that the matrix $L(x)$ given in (3.28) is symmetric even though the tangential derivatives do not necessarily commute, as by interchanging the derivatives we obtain

[TABLE]

and since $\underline{D}_{m}^{\Gamma_{0}}\phi\cdot\left(\nu^{\Gamma}\circ\phi\right)=0$, for all $m=1,...,n+1$, we see that the last two terms vanish.

3.5. The normal derivative at the boundary

We conclude this section by computing the pull-back of the normal derivative at the boundary for functions defined over the parametrised bulk-surface described in (3.10).

Lemma 3.7** (Normal derivative).**

Given any $u:\bar{D}\rightarrow\mathbb{R}$ sufficiently smooth, the pull-back of its normal derivative is given by

[TABLE]

where $G=\nabla\phi^{\top}\nabla\phi$ and $g=det(G)$ denoting its determinant.

Proof.

Differentiating $u=\hat{u}\circ\phi^{-1}$ and substituting in the expression (3.25) for the pull-back of the unit normal $\nu^{\Gamma}$, where the orientation has been chosen to be in the outer direction to the domain $D$ gives

[TABLE]

We next observe with the decomposition $\nabla\phi=\nabla_{\Gamma_{0}}\phi+\frac{\partial\phi}{\partial\nu_{\Gamma_{0}}}\otimes\nu^{\Gamma_{0}}$ and the orthogonality result $\nabla_{\Gamma_{0}}\phi^{\top}(\nu^{\Gamma}\circ\phi)=0$ that

[TABLE]

Since $\phi$ maps the boundary $\Gamma_{0}$ onto $\Gamma$, it follows that $\frac{\partial\phi}{\partial\nu_{\Gamma_{0}}}\cdot(\nu^{\Gamma}\circ\phi)>0$ and thus

[TABLE]

We now continue by showing that the normal component of $\frac{\partial\phi}{\partial\nu_{\Gamma_{0}}}$ can be expressed as the ratio between the bulk $\sqrt{g}$ and the surface area element $\sqrt{g_{\Gamma_{0}}}$. This will be achieved in the context of exterior algebras.

Let $\tau_{1},...,\tau_{n}$ represent an orthonormal basis of the tangent space $T\Gamma_{0}$ and thus $\{\tau_{1},...,\tau_{n},\nu^{\Gamma_{0}}\}$ forms a basis of $\mathbb{R}^{n+1}$. The determinant of linear map corresponding to $\nabla\phi$ evaluated on the boundary $\Gamma_{0}$ can be expressed in the notation of exterior algebras as follows

[TABLE]

Hence it follows

[TABLE]

We thus obtain the stated result with the following observations

[TABLE]

∎

4. First applications of the domain mapping method to complex random geometries involving random surfaces

We will now consider two model elliptic problems posed on complex random domains involving random surfaces. In particular, the first problem will be posed on a sufficiently smooth random surface and the second on a random bulk-surface. In both cases, the complete random domain mapping will be assumed to be known. Furthermore, we will assume that the computational domain was chosen to coincide with the expected domain, and thus will assume in both cases that $\mathbb{E}[\phi]=0$. We will now employ the domain mapping method, and reformulate both equations onto their corresponding expected domain and prove well-posedness as well as a regularity result.

4.1. An elliptic equation on a random surface

Let $\Gamma(\omega)$ represent a random, compact $C^{2}-$hypersurface in $\mathbb{R}^{n+1}$ prescribed by

[TABLE]

for a given random field $\phi\in L^{\infty}(\Omega;C^{2}(\Gamma_{0};\mathbb{R}^{n+1}))$ defined over a fixed, compact $C^{2}-$hypersurface $\Gamma_{0}\subset\mathbb{R}^{n+1}$. We will assume that the random domain mapping $\phi(\omega,\cdot):\Gamma_{0}\rightarrow\Gamma(\omega)$ is a $C^{2}-$ diffeomorphism for almost every $\omega$ and furthermore satisfies the uniform bounds

[TABLE]

for some constant $C>0$ independent of $\omega$. We consider the following model elliptic equation on the random surface

[TABLE]

for a given random field $f(\omega,\cdot):\Gamma(\omega)\rightarrow\mathbb{R}.$ Our goal is to analyse the mean solution defined by

[TABLE]

Reformulating (4.3) onto the expected domain with the calculation of the Laplace-Beltrami operator provided in Lemma 3.2 yields

[TABLE]

where the random coefficient is given by

[TABLE]

with $g_{\Gamma_{0}}(\omega)=\det G_{\Gamma_{0}}(\omega)$. Multiplying through by surface area element $\sqrt{g_{\Gamma_{0}}(\omega)}$ and integrating by parts, we arrive at the following mean-weak formulation on the fixed deterministic domain $\Gamma_{0}$.

Problem 4.1** (Mean-weak formulation).**

Given $\hat{f}\in L^{2}(\Omega;L^{2}(\Gamma_{0}))$, find $\hat{u}\in L^{2}(\Omega;H^{1}(\Gamma_{0}))$ such that

[TABLE]

for every $\hat{\varphi}\in L^{2}(\Omega;H^{1}(\Gamma_{0}))$. Here, we have set $\mathcal{D}_{\Gamma_{0}}(\omega)=\sqrt{g_{\Gamma_{0}}(\omega)}G_{\Gamma_{0}}^{-1}(\omega).$

We denote the associated bilinear form $a(\cdot,\cdot):L^{2}(\Omega;H^{1}(\Gamma_{0}))\times L^{2}(\Omega;H^{1}(\Gamma_{0}))\rightarrow\mathbb{R}$ and linear functional $l(\cdot):L^{2}(\Omega;L^{2}(\Gamma_{0}))\rightarrow\mathbb{R}$ by

[TABLE]

Thus the mean-weak formulation can be written more succiently as

[TABLE]

Proposition 4.1**.**

Under the uniformity assumptions (4.2) on the random domain mapping, there exists constants $C_{D_{\Gamma_{0}}},C_{g_{\Gamma_{0}}}>0$ such that the singular values $\sigma_{i}$ of $\mathcal{D}_{\Gamma_{0}}$ and the surface area element $\sqrt{g_{\Gamma_{0}}}$ are bounded above and below by

[TABLE]

for all $x\in\Gamma_{0}\text{ and a.e. }\omega$.

Proof.

We can rewrite $G_{\Gamma_{0}}$ using the orthogonality $\nabla_{\Gamma_{0}}\phi^{\top}(\nu^{\Gamma}\circ\phi)=0,$ as follows

[TABLE]

Examining each term separately, we see that the inverse is given by

[TABLE]

Hence it follows

[TABLE]

Therefore with (4.2), we have uniform bounds above and below on the singular values of $G_{\Gamma_{0}}(\omega)$ and hence obtain the estimates (4.10) and (4.11). ∎

A direct consequence of the above uniform bounds on the random coefficients is the existence and uniqueness of a solution to (4.6) gauranteed by the Lax-Milgram theorem.

Theorem 4.1**.**

Given any $\hat{f}\in L^{2}(\Omega;L^{2}(\Gamma_{0}))$, there exists a unique solution $\hat{u}$ to the mean-weak formulation (4.6) that satisfies the energy estimate

[TABLE]

Proof.

The stability estimate (4.12) follows from the coercivity of $a(\cdot,\cdot)$. ∎

By considering the original surface equation (4.3) on $\Gamma(\omega)\in C^{2}$, we would expect from standard elliptic surface regularity results that for given $f(\omega)\in L^{2}(\Gamma(\omega))$, the pathwise solution belongs to $u(\omega)\in H^{2}(\Gamma(\omega))$ and therefore $\hat{u}(\omega)\in H^{2}(\Gamma_{0})$ for a.e. $\omega$. However since the $H^{2}$ a-priori estimate on $u(\omega)$ will naturally depend on the geometry of the realisation $\Gamma(\omega)$, it is not immediately clear whether the solution to the mean-weak formulation belongs to $\hat{u}\in L^{2}(\Omega;H^{2}(\Gamma_{0}))$. We will therefore continue by explicitly treating all arising constants and their dependency on the geometry of the random domain.

Theorem 4.2** (Regularity).**

Given any $\hat{f}\in L^{2}(\Omega;L^{2}(\Gamma_{0}))$, the solution to (4.6) belongs to $\hat{u}\in L^{2}(\Omega;H^{2}(\Gamma_{0}))$ and furthermore satisfies the following estimate

[TABLE]

Proof.

Let us consider the push-forward $u=\hat{u}\circ\phi^{-1}$ of realisations of the weak solution onto $\Gamma(\omega)$ for almost every $\omega$, which as a result of the tensor structure $L^{2}(\Omega;H^{1}(\Gamma_{0}))\cong L^{2}(\Omega)\otimes H^{1}(\Gamma_{0})$ is a pathwise weak solution of

[TABLE]

with $f=\hat{f}\circ\phi^{-1}.$ Since for almost every $\omega\in\Omega$, $\Gamma(\omega)$ is $C^{2}$ and $f(\omega)\in L^{2}(\Gamma(\omega))$, it follows that $u(\omega)\in H^{2}(\Gamma(\omega))$ and therefore $\hat{u}(\omega)\in H^{2}(\Gamma_{0})$. For the a-priori estimate (4.13), it was shown in [14] through a series of integration by parts and interchanging of tangential derivatives that the $H^{2}$ semi-norm satisfies

[TABLE]

with

[TABLE]

Here $H^{\Gamma(\omega)}=trace\left(\mathcal{H}^{\Gamma(\omega)}\right)$ is the mean-curvature. Hence with the uniform bounds (4.2) on the random domain mapping and the previously calculated expression (3.29) for the Weingarten map, we obtain an upper bound on the constant $c(\omega)$ independent of $\omega$. Thus, with the PDE (4.14) pointwise we have the bound

[TABLE]

We can now pull-back onto the expected domain, applying the norm equivalence of the pull-back transformation

[TABLE]

where the constants are independent of $\omega$ due to bounds (4.2), and the stability estimate (4.12) to obtain

[TABLE]

and thus the stated result. ∎

4.2. A coupled elliptic system on a random bulk-surface

For the second problem, we consider a coupled elliptic system on a random bulk-surface motivated by the deterministic case analysed in [12]. More precisely, the geometric setting is as follows. We let $\{\Gamma(\omega)\}$ denote a family of random, compact $C^{2}-$hypersurfaces in $\mathbb{R}^{n+1}$ enclosing open domains $D(\omega)$ and will denote the outer unit normal by $\nu^{\Gamma(\omega)}.$ The family of random domains will be prescribed by the mapping

[TABLE]

where the reference surface $\Gamma_{0}\subset\mathbb{R}^{n+1}$ will also be a compact $C^{2}-$hypersurface with open interior $D_{0}$. We will assume that the domain mapping is a $C^{2}-$diffeomorphism for $a.e.$ $\omega\in\Omega$ and additionally satisfies

[TABLE]

for a constant $C>0$ independent of $\omega$. The proposed coupled elliptic system on the random bulk-surface reads as follows

[TABLE]

Here $\alpha,\beta>0$ are given positive constants and $f(\omega,\cdot):D(\omega)\rightarrow\mathbb{R}$ and $f_{\Gamma}(\omega,\cdot):\Gamma(\omega)\rightarrow\mathbb{R}$ are prescribed random fields. As with our previous problem, our quantity of interest is the mean solution, that is the pair $\left(\mathbb{E}[u],\mathbb{E}[v]\right)$ defined by

[TABLE]

Let us continue by reformulating the system (4.17) onto the expected domain $\overline{D_{0}}$ with our previously calculated expressions for the Laplace-Beltrami operator (3.21) and the normal derivative (3.30) giving

[TABLE]

Here the random coefficients are

[TABLE]

with $g(\omega)=detG(\omega)$, $g_{\Gamma_{0}}(\omega)=detG_{\Gamma_{0}}(\omega)$. For convenience, we have set $\hat{f}_{\Gamma_{0}}=f_{\Gamma}\circ\phi$. To derive a mean-weak formulation, we follow the variational approach presented in [12]. We begin by multiplying through the bulk equation (4.18a) by the area element $\sqrt{g}$ and integrating by parts which gives

[TABLE]

Similarly, for the surface equation (4.18c) we integrate by parts recalling that the hypersurface $\Gamma_{0}$ is without boundary, to obtain

[TABLE]

Taking the weighted sum and substituting in the reformulated Robin boundary condition (4.18b), we arrive at the following mean-weak formulation:

Problem 4.2** (Mean-weak formulation).**

Given any $\hat{f}\in L^{2}(\Omega;L^{2}(D_{0}))$ and $\hat{f}_{\Gamma_{0}}\in L^{2}(\Omega;L^{2}(\Gamma_{0}))$, find $\hat{u}\in L^{2}(\omega;H^{1}(D_{0}))$ and $\hat{v}\in L^{2}(\Omega;H^{1}(\Gamma_{0}))$ such that

[TABLE]

for every $\hat{\varphi}\in L^{2}(\Omega;H^{1}(D_{0}))$ and $\hat{\xi}\in L^{2}(\Omega;H^{1}(\Gamma_{0})).$ Here we set $\mathcal{D}_{\Gamma_{0}}(\omega)=\sqrt{g_{\Gamma_{0}}(\omega)}G_{\Gamma_{0}}^{-1}(\omega)$.

We denote the associated bilinear form and linear functional stated above by

[TABLE]

where we have set $H=L^{2}(D_{0})\times L^{2}(\Gamma_{0})$ and $V=H^{1}(D_{0})\times H^{1}(\Gamma_{0})$ to be Hilbert spaces equipped with respective inner products

[TABLE]

The mean-weak formulation thus reads as follows

[TABLE]

The following uniform bounds on the random bulk coefficients follow immediately from the assumption (4.16) on the random domain mapping. Furthermore, the derived bounds on the surface coefficients presented in Proposition 4.1 also hold since the tangential derivatives of the surface parametrisation and its inverse are also uniformly bounded as a consequence of (4.16).

Proposition 4.2** (Uniform bounds).**

There exists constants $C_{g},C_{D}>0$ such that the bulk area element $\sqrt{g(\omega)}$ and the singular values $\sigma_{i}$ of $D(\omega)$ are uniformly bounded for all $x\in D_{0}$ and a.e. $\omega$ by

[TABLE]

Theorem 4.3**.**

Given any $(\hat{f},\hat{f}_{\Gamma_{0}})\in H$, there exist a unique solution $(\hat{u},\hat{v})\in L^{2}(\Omega;V)$ to (4.22) which satisfies the energy estimate

[TABLE]

Proof.

With our uniform bounds (4.23), (4.10) on the random bulk and surface coefficients, we can now proceed in verifying all the conditions of the Lax-Milgram theorem are satisified. For a coercivity estimate, we argue

[TABLE]

For the continuity of the bilinear form $a(\cdot,\cdot)$, we apply the Cauchy-Schwarz inequality with the boundedness of the trace operator $\|f\|_{L^{2}(\Gamma_{0})}\leq c_{T}\|f\|_{H^{1}(D_{0})}$ as follows

[TABLE]

Thus we have the existence and uniqueness of a solution to (4.22). The estimate (4.25) then follows from coercivity of $a(\cdot,\cdot).$ ∎

Theorem 4.4** (Regularity).**

Given any $\hat{f}\in L^{2}(\Omega;L^{2}(D_{0}))$ and $\hat{f}_{\Gamma_{0}}\in L^{2}(\Omega;L^{2}(\Gamma_{0}))$, the mean-weak solution $(\hat{u},\hat{v})$ to (4.22) satisfies

[TABLE]

Furthermore, we have

[TABLE]

where the constant $C>0$ depends only the geometry of the reference domain $\overline{D_{0}}$ and the uniform bound (4.16) on the random domain mapping.

Proof.

Observe that for $a.e.$ $\omega\in\Omega$, the solution $(\hat{u},\hat{v})$ satisfies for every $\hat{\varphi}\in H^{1}(D_{0})$ and $\hat{\xi}\in H^{1}(\Gamma_{0})$,

[TABLE]

Setting $\hat{\varphi}=0$ gives

[TABLE]

Hence we see that $\hat{v}(\omega)$ is the pathwise weak solution to the elliptic surface equation

[TABLE]

It therefore follows form the surface regularity result given in Theorem 4.2 since $\hat{u}(\omega)\in L^{2}(\Gamma_{0})$, that $\hat{v}(\omega)\in H^{2}(\Gamma_{0})$ for a.e. $\omega$ and furthermore

[TABLE]

where the constant $C>0$ is independent of $\omega$. To obtain higher regularity of the bulk quantity, we set $\hat{\xi}=0$ yielding

[TABLE]

This is precisely the weak formulation of the following elliptic boundary value problem subject to the reformulated Robin boundary condition

[TABLE]

Since the coefficients are sufficiently regular, more precisely

[TABLE]

and the boundary is sufficiently smooth $\Gamma_{0}\in C^{2}$, we can apply standard regularity results [22] to deduce $\hat{u}(\omega)\in H^{2}(D_{0})$ for a.e. $\omega$ with the estimate

[TABLE]

Here the constant $C>0$ is independent of $\omega$ since all the coefficients are uniformly bounded and furthermore, $\mathcal{D}(\omega)$ is uniformly elliptic in $\omega$. Combining (4.28) and (4.29) with the stability estimate (4.25) and boundedness of the trace operator leads to

[TABLE]

and hence the stated result. ∎

5. An abstract numerical analysis of elliptic equations on random curved domains

We continue by considering in an abstract setting, the mean-weak formulation of general elliptic equations on random curved domains after being transformed onto the expected domain via the given stochastic domain mapping. Working in this abstract framework, we will present and analyse a finite element discretisation coupled with the Monte-Carlo method to approximate our quantity of interest, the mean solution. As the expected domain is assumed to be curved, the proposed finite element method will involve perturbations of the variational set up corresponding to the approximation of the domain. An optimal error bound in the energy norm for our non-conforming approach is derived with the help of the first lemma of Strang with suitable assumptions on the finite element space approximation and arising consistency error. Furthermore, an $L^{2}(\Omega;L^{2})$-type estimate is proved by a standard duality argument.

5.1. Abstract mean-weak formulation

Let $V$ and $H$ denote separable Hilbert spaces for which the embedding $V\hookrightarrow H$ is dense and continuous. We assume that we are in the setting where we have a sample dependent bilinear form $\tilde{a}(\omega;\cdot,\cdot):V\times V\rightarrow\mathbb{R}$ and linear functional $\tilde{l}(\omega;\cdot):H\rightarrow\mathbb{R}$ corresponding to the path-wise weak formulation

[TABLE]

of the elliptic equation after being reformulated onto the expected domain. For convenience, we will omit the pull-back notation for functions $\hat{u}$ since all the subsequent analysis will be considered on the expected domain. The mean-weak formulation will thus in general read as follows:

Problem 5.1** (Mean-weak formulation).**

Find $u\in L^{2}(\Omega;V)$ such that for every $\varphi\in L^{2}(\Omega;V)$ we have

[TABLE]

We denote the associated bilinear form $a(\cdot,\cdot):L^{2}(\Omega;V)\times L^{2}(\Omega;V)\rightarrow\mathbb{R}$ and linear functional $l(\cdot):L^{2}(\Omega;H)\rightarrow\mathbb{R}$ by

[TABLE]

and shall assume all the requirements of the Lax-Milgram theorem are satisfied thus ensuring the existence and uniqueness of the solution.

5.2. Abstract formulation of the finite element discretisation

For a given $h\in(0,h_{0})$, let $\mathcal{V}_{h}$ be a finite dimensional space that will represent a finite element space and let $V_{h}$ and $H_{h}$ denote the space $\mathcal{V}_{h}$ endowed with respective norms $\|\cdot\|_{V_{h}}$ and $\|\cdot\|_{H_{h}}$. We assume that $V_{h}$ and $H_{h}$ are Hilbert spaces and furthermore that $V_{h}\hookrightarrow H_{h}$ is uniformly embedded, that is

[TABLE]

for a constant $c>0$ independent of $h.$ In practice, the spaces $V_{h}$ and $H_{h}$ will represent equivalent Hilbert spaces to the continuous solution spaces $V$ and $H$ but posed over a discrete approximation of the curved domain, with $h$ denoting the discretisation parameter. We introduce the sample-dependent bilinear form and linear functional

[TABLE]

that are perturbations approximating their continuous counterparts and will assume $\tilde{a}_{h}(\omega:\cdot,\cdot)$ is uniformly $V_{h}$-elliptic and bounded and additionally $\tilde{l}_{h}(\omega;\cdot)$ is uniformly bounded. More precisely, there exists constants $c_{1},c_{2},c_{3}>0$ independent of $\omega$ and $h$ such that

[TABLE]

The finite element approximation of the mean-weak formulation (5.1) for a given a finite dimensional subspace $\mathcal{V}_{h}\subset V_{h}$ will then take the following form:

Problem 5.2** (Semi-discrete problem).**

Find $U_{h}\in L^{2}(\Omega;\mathcal{V}_{h})$ such that

[TABLE]

for all $\phi_{h}\in L^{2}(\Omega;\mathcal{V}_{h}).$

By our uniform assumptions of the bilinear form $\tilde{a}(\omega;\cdot,\cdot)$ and the linear functional $\tilde{l}(\omega;\cdot)$, we deduce the existence and uniqueness of a solution to the semi-discrete problem.

Theorem 5.1**.**

There exists a unique solution $U_{h}\in L^{2}(\Omega;V_{h})$ to the semi-discrete problem (5.5) that satisfies

[TABLE]

with the constant $C>0$ independent of $h\in(0,h_{0}).$

Observe that if we let $\{\chi_{j}\}_{j=1}^{N}$ be a basis of $\mathcal{V}_{h}$ and express $U_{h},\phi_{h}\in L^{2}(\Omega;\mathcal{V}_{h})\cong L^{2}(\Omega)\otimes\mathcal{V}_{h}$ in the form

[TABLE]

where $U(\omega)=(U_{1}(\omega),...,U_{N}(\omega))^{\top}\in L^{2}(\Omega)^{N}$ and $\Phi(\omega)=(\phi_{1}(\omega),...,\phi_{N}(\omega))^{\top}\in L^{2}(\Omega)^{N}$, then (5.5) can be rewritten as

[TABLE]

Here the random stiffness matrix $S(\omega)=(S_{ij}(\omega))_{i,j=1,...,N}$ and load vector $F(\omega)=(F_{j}(\omega))_{j=1,...,N}$ are given by $S_{ij}(\omega)=\tilde{a}_{h}(\omega;\chi_{j},\chi_{i}),$ $F_{j}(\omega)=\tilde{l}_{h}(\omega;\chi_{j}).$ Since $\phi_{j}(\omega)\in L^{2}(\Omega)$ are arbitrary, we deduce that the semi-discrete problem is equivalent to finding $U\in L^{2}(\Omega;\mathbb{R}^{N})$ which satisfies

[TABLE]

5.3. Assumptions on the finite element approximation and the continuous equations

We now state all the necessary assumptions that will be required in deriving an error estimate for the semi-discrete solution. In order to compare our semi-discrete solution with the continuous solution, we first need to assume the existence of a lifting map.

Assumption 5.1** (Lifting map).**

There exists a linear mapping $\Lambda_{h}:\mathcal{V}_{h}\rightarrow V$ for which there exists constants $c_{1},c_{2}>0$ independent of $h\in(0,h_{0})$ such that for all $\chi_{h}\in\mathcal{V}_{h}$

[TABLE]

We denote the lifted finite dimensional space by $V_{h}^{l}:=\Lambda_{h}\mathcal{V}_{h}$. Next, we introduce the Hilbert space $Z_{0}\hookrightarrow V$ which shall represent a space consisting of functions of higher regularity for which we assume we have the following interpolation estimate.

Assumption 5.2** (Approximation of finite element space).**

There exists a well-defined interpolation operator $I_{h}:Z_{0}\rightarrow V_{h}^{l}$ for which there exists $c>0$ such that

[TABLE]

Naturally, the lifting map and interpolation operator can be extended to random functions in a pathwise sense

[TABLE]

and the previous estimates (L1),(L2), (I1) hold for their respective norms $\|\cdot\|_{L^{2}(\Omega;H)}$ and $\|\cdot\|_{L^{2}(\Omega;V)}$. We continue by imposing bounds on the consistency error arising from the pertubation of the variational form. For this, we will assume the existence of an inverse lifting map $\Lambda_{h}:L^{2}(\Omega;Z_{0})\rightarrow L^{2}(\Omega;V_{h})$ and will denote inverse lift of a function $w$ by $w^{-l}$.

Assumption 5.3** (Consistency error).**

Given any $W_{h},\phi_{h}\in L^{2}(\Omega;\mathcal{V}_{h})$ with corresponding lifts $w_{h},\chi_{h}\in L^{2}(\Omega;V_{h}^{l})$, we have the bounds

[TABLE]

Furthermore, for any $w,\varphi\in L^{2}(\Omega;Z_{0})$ with inverse lifts $w^{-l},\varphi^{-l}$ we have

[TABLE]

Our final assumption will be on the regularity of an associated dual problem that will enable us to derive an $L^{2}(\Omega;H)$ error estimate using the standard Aubin-Nitsche trick. The associated dual problem reads as follows:

Problem 5.3** (Dual problem).**

For a given $g\in L^{2}(\Omega;H)$, find $w(g)\in L^{2}(\Omega;V)$ such that

[TABLE]

Here $(\cdot,\cdot)_{L^{2}(\Omega;H)}$ denotes the inner product on the Hilbert space $L^{2}(\Omega;H)$.

Assumption 5.4** (Regularity of dual problem).**

The solution $w(g)$ to the dual problem belongs to space $L^{2}(\Omega;Z_{0})$ and furthermore satisfies

[TABLE]

for a constant $c>0$ independent of both $g$ and $h\in(0,h_{0})$.

5.4. Error estimates for the semi-discrete solution

Recall that the abstract finite element space $\mathcal{V}_{h}$ is not necessarily contained in the Hilbert space $V$. However, with the assumed existence of a lifting map

[TABLE]

we can lift the discrete bilinear form $a_{h}(\cdot,\cdot)$ and the linear functional $l_{h}(\cdot)$ onto the space $L^{2}(\Omega;V_{h}^{l})$ by the following relations for $w_{h}=\Lambda_{h}W_{h},\varphi_{h}=\Lambda_{h}\phi_{h}\in L^{2}(\Omega;V_{h}^{l})$

[TABLE]

thus inducing a third variational problem equivalent to (5.5).

Problem 5.4** (Lifted semi-discrete problem).**

Find $u_{h}\in L^{2}(\Omega;V_{h}^{l})$ such that for every $\varphi_{h}\in L^{2}(\Omega;V_{h}^{l})$ we have

[TABLE]

Since $L^{2}(\Omega;V_{h}^{l})$ is contained in the solution space $L^{2}(\Omega;V)$, the lifted semi-discrete problem fits into the abstract non-conforming finite element setting considered in the first lemma of Strang [30]. We will now present these results in the context of our random Hilbert space setting.

Lemma 5.1** (First lemma of Strang).**

Let $u_{h}$ denote the solution to the lifted semi-discrete problem (5.11) and assume that the bilinear form $a_{h}^{l}(\cdot,\cdot)$ is uniformly $L^{2}(\Omega;V_{h}^{l})$-elliptic, i.e. for some $\alpha>0$

[TABLE]

for all $\varphi\in L^{2}(\Omega;V_{h}^{l})$ and $h\in(0,h_{0})$. Then there exists a constant $C>0$ independent of $h$ such that

[TABLE]

Theorem 5.2** (Error estimates).**

Let $u$ denote the solution of the continuous problem (5.1) and assume that it is sufficiently regular $u\in L^{2}(\Omega;Z_{0})$ and let $U_{h}$ be the discrete solution of (5.5) with lift $u_{h}=\Lambda_{h}U_{h}.$ Then with the assumptions listed in section 5.3 satisfied, there exists a constant $c>0$ such that for all $h\in(0,h_{0})$ we have the error estimate

[TABLE]

Proof.

It follows from the uniform ellipticity assumption (5.2) on the bilinear form $a_{h}(\cdot,\cdot)$ and the norm equivalence of the lifting map, that for any $\varphi_{h}=\Lambda_{h}\phi_{h}\in L^{2}(\Omega;V_{h}^{l})$ we have

[TABLE]

Therefore the bilinear form $a_{h}^{l}(\cdot,\cdot)$ is uniformly coercive and thus we can apply the first lemma of Strang. Substituting $\varphi_{h}=I_{h}u$ into the estimate (5.12) and inserting the consistency bounds (P1), (P2) gives

[TABLE]

Hence with the interpolation estimate (I1) applied to $u\in L^{2}(\Omega;Z_{0})$ we obtain

[TABLE]

For the $L^{2}(\Omega;H)-$estimate, we use a standard duality argument. Given $g\in L^{2}(\Omega;H)$ and an arbitrary $w_{h}\in L^{2}(\Omega;V_{h}^{l})$ we have

[TABLE]

Choosing $w_{h}=I_{h}w(g)$ and applying the interpolation estimate (I1) to the solution of the dual problem which is assumed (R1) to be sufficiently regular $w(g)\in L^{2}(\Omega;Z_{0})$ gives

[TABLE]

We bound the consistency error in the second term with (P2) giving

[TABLE]

To obtain a bound of order $h^{2}$ for the third term, we begin by rewriting it as follows

[TABLE]

Now we are able to apply the estimate (P3) to the last term since both $u,w(g)\in L^{2}(\Omega;Z_{0})$ and can then follow a similar argument as to the previous cases for the first two terms which leads to

[TABLE]

Combining the results gives the stated result

[TABLE]

∎

We conclude our abstract error analysis by combining our finite element discretisation with the Monte-Carlo method to estimate our quantity of interest, the mean solution $E[u]$. Recall, that for an arbitrary Hilbert space $\mathcal{H}$, the Monte-Carlo estimator of the expectation of a random variable $Y\in L^{2}(\Omega;\mathcal{H})$ is a $\mathcal{H}$-valued random variable $E_{M}[Y]:\otimes_{i=1}^{M}\Omega\rightarrow\mathcal{H}$ defined by

[TABLE]

where $M\in\mathbb{N}$ is the chosen number of samples taken and $\hat{Y}_{i}$ are independent identically distributed copies of the random variable $Y$. Furthermore, we have the following well-known convergence result, see [24].

Lemma 5.2** (Monte-Carlo convergence rate).**

For a given $M\in\mathbb{N}$ and a $\mathcal{H}$-valued random variable $Y\in L^{2}(\Omega;\mathcal{H})$, the Monte-Carlo estimator satisfies the convergence rate

[TABLE]

Therefore, if we consider the error between the mean solution $\mathbb{E}[u]$ and our discrete approximation $\mathbb{E}[u_{h}]$ in the $L^{2}(\Omega^{M};H)$ norm, and decompose it into the error arising from the finite element discretisation and the statistical error for the Monte-Carlo approximation, we obtain the following bound

[TABLE]

A similar argument in the $L^{2}(\Omega;V)$ leads to the following convergence rates.

Theorem 5.3**.**

Let all the conditions from Theorem 5.2 be satisfied. Then we have the following error estimates

[TABLE]

6. Discretisation of the reformulated elliptic PDEs on their expected domains

In this section, we apply the results from the abstract theory to two finite element discretisation schemes for the reformulations of the two model elliptic equations. In each case, we will verify that all the listed assumptions in abstract setting are satisfied hence giving the stated convergence rate.

6.1. The elliptic equation on a random surface

To discretise the reformulation of the elliptic equation

[TABLE]

on the expected domain, we propose a semi-discrete scheme using linear Lagrangian surface finite elements [14]. Our computational domain $\Gamma_{h}$ approximating the smooth expected hypersurface $\Gamma_{0}$ will be a polyhedral surface

[TABLE]

consisting of finitely many non-degenerate triangles whose vertices are taken to lie on the surface $\Gamma_{0}$ and have the maximum diameter bounded above by $h>0$. The triangulation will be assumed to be shape regular and quasi-uniform, in the sense that the in-ball radius of each element is uniformly bounded below by $ch$, for some constant $c>0$. In order to lift functions between the continuous and discrete surface, we shall assume that the projective mapping $a:\Gamma_{h}\rightarrow\Gamma_{0}$ decribed in (3.5) is bijective and define the lift and inverse lift of functions $f$ and $g$ given over $\Gamma_{h}$ and $\Gamma_{0}$ respectively by

[TABLE]

where $x(a)$ denotes the inverse of the projection mapping $a$. We introduce the linear finite element space on $\Gamma_{h}$

[TABLE]

and define the lifted finite element space by

[TABLE]

The finite element discretisation of the mean-weak formulation reads as follows.

Problem 6.1** (Semi-discrete scheme).**

Find $U_{h}\in L^{2}(\Omega;S_{h})$ such that

[TABLE]

for every $\phi_{h}\in L^{2}(\Omega;S_{h}).$

In the context of the abstract framework, the finite dimensional space $\mathcal{V}_{h}$ is taken to be the finite element space $S_{h}$ and the Hilbert spaces $V_{h},H_{h}$ are given by $H^{1}(\Gamma_{h})$ and $L^{2}(\Gamma_{h})$. Furthermore, the abstract sample-dependent discrete bilinear form $\tilde{a}_{h}(\omega;\cdot,\cdot):H^{1}(\Gamma_{h})\times H^{1}(\Gamma_{h})\rightarrow\mathbb{R}$ and linear functional $\tilde{l}(\omega;\cdot):L^{2}(\Gamma_{h})\rightarrow\mathbb{R}$ are given by

[TABLE]

With the uniform bounds on the random coefficients (4.10), (4.11), we deduce that $\tilde{a}_{h}(\omega:\cdot,\cdot)$ is uniformly $L^{2}(\Omega;H^{1}(\Gamma_{0}))$-elliptic and bounded, and additionally $\tilde{l}(\omega;\cdot)$ is uniformly bounded as presumed in (5.2 - 5.4), and hence obtain existence and uniqueness of a semi-discrete solution to (6.4). We continue by checking the stated assumptions in the abstract error analysis. In particular, we begin with the norm equivalence (L1),(L2) of the lifting map $\Lambda_{h}:\mathcal{V}_{h}\rightarrow V$ given by $\Lambda_{h}\chi_{h}=\chi_{h}^{l}$. A proof of these estimates can be found in [14, Lemma 4.2].

Lemma 6.1** (Equivalence in norms of lifts).**

There exists constants $c_{1},c_{2}>0$ independent of $h$ such that for any $\chi_{h}\in S_{h}$ with lift $\chi_{h}^{l}\in S_{h}^{l}$ we have

[TABLE]

For the interpolation assumption (I1), we set the Hilbert space $Z_{0}$ consisting of functions of higher regularity to be $H^{2}(\Gamma_{0})$. It follows from the Sobolev embedding that $H^{2}(\Gamma_{0})\subset C^{0}(\Gamma_{0})$ for $n\leq 3$ and therefore we can introduce the interpolation operator $I_{h}:H^{2}(\Gamma_{0})\rightarrow S_{h}^{l}$ defined by

[TABLE]

where $\hat{I}_{h}:C^{0}(\Gamma_{h})\rightarrow S_{h}$ denotes the standard Lagrangian interpolatant defined element-wise on $\Gamma_{h}$. The following estimate was proved in [14, Lemma 4.3].

Lemma 6.2** (Interpolation estimate).**

Given any $\eta\in H^{2}(\Gamma_{0})$, there exists a constant $c>0$ independent of $h$ such that

[TABLE]

To derive the assumed bounds (P1),(P2) and (P3) on the approximation of the discrete bilinear forms, we first need a preliminary result on the order of approximation of the geometry, see [14, Lemma 4.1].

Lemma 6.3** (Geometric error bounds).**

Let $\delta_{h}^{\Gamma_{0}}$ denote the surface element corresponding to the transformation from $\Gamma_{0}$ to $\Gamma_{h}$ under the lifting map $d\sigma(a(x))=\delta_{h}(x)d\sigma_{h}(x)$ and define

[TABLE]

where $\mathcal{P}_{h}:=I-\nu_{h}\otimes\nu_{h}$ is the projection operator mapping onto the tangent space of the discrete surface $\Gamma_{h}$ defined element-wise. Then we have the estimates

[TABLE]

We can now bound the consistency error as follows.

Lemma 6.4** (Consistency error).**

*Given any $(W_{h},\phi_{h})\in L^{2}(\Omega;S_{h})\times L^{2}(\Omega;S_{h})$ with lifts

$(w_{h},\varphi_{h})\in L^{2}(\Omega;S_{h}^{l})\times L^{2}(\Omega;S_{h}^{l})$, we have*

[TABLE]

Proof.

Lifting the discrete integral in the linear functional $l_{h}(\cdot)$ onto the smooth surface $\Gamma_{0}$ with the projective mapping $a(\cdot)$ leads to

[TABLE]

Hence with the uniform bound (4.11) on the random coefficient $\sqrt{g_{\Gamma_{0}}(\omega)}$ and the order $h^{2}$ approximation of the geometric pertubation (6.9), we obtain the estimate (6.11). For (6.12), we begin by applying the chain rule to lift $W_{h}(\omega,x)=w_{h}(\omega,a(x))$

[TABLE]

Suppressing the parameter $x$, we deduce

[TABLE]

Therefore, we can express the pertubation error in the approximation of the bilinear form $a(\cdot,\cdot)$ by

[TABLE]

and hence with the uniform bounds (4.10), (4.11) on the random coefficients and the geometric estimates (6.9), (6.10) we obtain (6.12). ∎

For the regularity assumption (R1) on the associated dual problem

[TABLE]

which due the symmetry of $\mathcal{D}_{\Gamma_{0}}$ and thus of $a(\cdot,\cdot)$, is precisely the mean-weak formulation, we have the results presented in Theorem 4.2.

6.2. The coupled elliptic system

We next apply the results from the abstract framework to the second model problem of the coupled elliptic system

[TABLE]

on a random bulk-surface. Our proposed finite element discretisation of the system reformulated on the expected domain and the subsequent analysis will be based on the approach presented in [12]. For the computational domain, we approximate the open bulk $D_{0}\subset\mathbb{R}^{n+1}$ by a polyhedral domain

[TABLE]

consisting of closed $(n+1)-$simplices with maximum diameter uniformly bounded above by positive constant $h>0$ and will assume that the triangulation $\mathcal{T}_{h}$ is quasi-uniform. We denote the induced discrete surface $\Gamma_{h}=\partial D_{h}$ and the associated triangulation by

[TABLE]

and impose the same assumptions on $\mathcal{T}_{h}$ as were listed in the previous example. A piece-wise diffeomorphic mapping $G_{h}:D_{h}\rightarrow D_{0}$ from the discrete bulk to the continuous can be constructed by fixing the interior simplices (simplices with at most one vertex on the boundary $\Gamma_{0}$) and using the projective mapping $a^{\Gamma_{0}}(\cdot)$ to define a diffeomorphism $\Lambda_{h,k}:K\rightarrow K^{e}$ between the boundary simplices $K$ (simplices with at least two vertices on $\Gamma_{0}$) and the exact curved simplices $K^{e}$,

[TABLE]

Details on the precise form of $\Lambda_{h,K}$ can be found in [12]. We are therefore able to define lifts and inverse lifts of functions on the bulk domain by

[TABLE]

Note that, the diffeomorphism $\Lambda_{h,K}$ is chosen such that the mapping $G_{h}$ coincides with the projective mapping

[TABLE]

on the boundary of the discrete bulk and hence the bulk lift agrees with the surface lifting map described in (6.1) on $\partial D_{h}$. For convenience, we will denote the sub-triangulation consisting of all boundary simplices by

[TABLE]

and define the corresponding sets

[TABLE]

where the lifting maps $G_{h},G_{h}^{-1}$ differ from the identity mapping. We introduce the linear finite element spaces on the discrete bulk and discrete surface by

[TABLE]

and denote the corresponding lifted finite element spaces by

[TABLE]

An important feature of our finite element spaces is that the trace of a function $\phi_{h}\in V_{h}$ belongs to $S_{h}$ and similarly the trace of $\varphi_{h}\in V_{h}^{l}$ belongs to $S_{h}^{l}$ as a result of (6.16). The finite element discretisation of the mean-weak formulation then reads as follows.

Problem 6.2** (Semi-discrete problem).**

Find a pair $(U_{h},V_{h})\in L^{2}(\Omega;V_{h}\times S_{h})$ such that

[TABLE]

for every $(\phi_{h},\zeta_{h})\in L^{2}(\Omega;V_{h}\times S_{h}).$

Here the abstract finite dimensional space is $\mathcal{V}_{h}=V_{h}\times S_{h}$ and the Hilbert spaces $V_{h},H_{h}$ are given by $H^{1}(D_{0})\times H^{1}(\Gamma_{0})$ and $L^{2}(D_{0})\times L^{2}(\Gamma_{0})$ respectively. We denote the associated bilinear form and linear functional

[TABLE]

to be the respective left hand side and right hand side of the semi-discrete variational problem 6.2. By the uniform bounds on the random coefficients (4.23), (4.10), we deduce the existence and uniqueness of a semi-discrete solution using a similar argument to the continuous problem. We proceed in a similar manner and check that the assumptions of the abstract analysis are satisfied. The norm equivalence (L1), (L2) of the lifting mapping which in this setting $\Lambda_{h}:V_{h}\times S_{h}\rightarrow V_{h}^{l}\times S_{h}^{l}$ is given component-wise by

[TABLE]

follows from the estimates on the surface lifting map given Lemma 6.1 in combination with the following bulk lifting norm equivalence derived in [12, Proposition 4.9].

Lemma 6.5** (Bulk lift estimates).**

There exists constants $c_{1},c_{2}>0$ independent of $h$, such that for any $\phi_{h}:D_{h}\rightarrow\mathbb{R}$ with lift $\varphi_{h}=\phi_{h}^{l}:D_{0}\rightarrow\mathbb{R}$ we have

[TABLE]

For the interpolation assumption (I1), we set the abstract function space $Z_{0}=H^{2}(D_{0})\times H^{2}(\Gamma_{0})$ and define the interpolation operator component-wise

[TABLE]

with $\tilde{I}_{h}$ denoting the standard Lagrangian intepolation operator and have the following estimate .

Lemma 6.6** (Interpolation estimate).**

There exists a well-defined interpolation operator

[TABLE]

such that for any $(\eta,\xi)\in H^{2}(D_{0})\times H^{2}(\Gamma_{0})$ we have

[TABLE]

The next step will entail bounding the consistency error arising from the geometric approximation of the domain. Estimates for the surface pertubation have previously been given in Lemma 6.3. For the bulk approximation, we recall that the lifting mapping $G_{h}:D_{h}\rightarrow D_{0}$ is defined to be the identity on interior simplices and a $C^{1}-$diffeomorphism for simplices near the boundary. Therefore the corresponding bulk error will be comprised of two parts; the first part will be related to the smallness of the neighbourhood around $\Gamma_{0}$ in which the lifted boundary simplices lie in and the second part is the associated geometric error of the boundary simplices approximating the corresponding exact curved simplex. We begin with the latter and state geometric bulk estimates on the diffeomorphic mapping $G_{h}$, for which a proof of the bounds (6.24) and (6.25) can be found in [12, Proposition 4.7].

Lemma 6.7** (Geometric bulk estimates).**

Let $\delta_{h}^{D_{0}}=|det(\nabla G_{h})|$ be the volume element corresponding to the transformation $G_{h}:D_{h}\rightarrow D_{0}$ and set

[TABLE]

Then we have the following estimates for a constant $c>0$ independent of $\omega$,

[TABLE]

Proof.

The estimate (6.26) follows from the observation

[TABLE]

and the uniform bounds (4.23) on the random coefficient $\mathcal{D}(\omega)$. ∎

To obtain a bound on the open neighbourhood containing the boundary simplices, we have the subsequent narrow band inequality [12, Lemma 4.10].

Lemma 6.8** (Narrow band trace inequality).**

Given any $\delta<\delta_{\Gamma_{0}}$, let $\mathcal{N}_{\delta}$ be a narrow band in the interior domain $D_{0}$ around the boundary $\Gamma_{0}$ defined by

[TABLE]

Then for any $\eta\in H^{1}(D_{0})$ we have

[TABLE]

The consistency error can now be bounded as follows.

Lemma 6.9** (Consistency error).**

Assume $f\in L^{2}(\Omega;H^{1}(D_{0}))$. Then for any $\phi_{h},W_{h}\in L^{2}(\Omega;V_{h})$ and $\zeta_{h},X_{h}\in L^{2}(\Omega;S_{h})$ with corresponding lifts $\varphi_{h},w_{h}$ and $\xi_{h},\chi_{h}$ we have

[TABLE]

Furthermore, for any $\varphi,w\in L^{2}(\Omega;H^{2}(D_{0}))$ and $\xi,\chi\in L^{2}(\Omega;H^{2}(\Gamma_{0}))$ with inverse lifts $\varphi^{-l},w^{-l}$ and $\xi^{-l},\chi^{-l}$ we have

[TABLE]

Proof.

For the estimate (6.28), we begin by lifting the discrete integrals in $l_{h}(\cdot)$ onto their respective continuous counterparts recalling that the set of all boundary simplices $B_{h}$ is the region in which the diffeomorphic mapping $G_{h}$ differs from the identity and thus where $\delta_{h}^{D_{0}}=det(\nabla G_{h})\neq 1$,

[TABLE]

Substituting the geometric bulk and surface estimates (6.25), (6.9) with the uniform bounds on the random coefficients (4.23), (4.10) leads to

[TABLE]

To obtain a bound of order $h^{2}$ on the bulk term, we will now apply the narrow trace band inequality. We choose $\delta>0$ such that $0<h<\delta<ch$ for some constant $c>0$, thus giving

[TABLE]

With a similar estimate on the test function $\varphi_{h}$, we obtain (6.28). For (6.29) and (6.30), we apply the chain rule to the lifts $\varphi_{h}(\omega,G_{h}(x))=\phi_{h}(\omega,x)$ and $w_{h}(\omega,G_{h}(x))=W_{h}(\omega,x)$ to deduce

[TABLE]

We can therefore express the perturbation error in our approximation of $a(\cdot,\cdot)$ as follows

[TABLE]

Here we have again used the fact that the diffeomorphic mapping $G_{h}$ is the identity on interior simplices and consequently $\delta_{h}^{D_{0}}=1$ and $R_{h}^{D_{0}}=I$ on $D_{h}\setminus B_{h}.$ We now apply the geometric estimates and bounds on the random coefficients to obtain

[TABLE]

For the last term, we observe by the boundedness of the trace operator $\|f\|_{L^{2}(\Gamma_{0})}\leq c_{T}\|f\|_{H^{1}(D_{0})}$ that

[TABLE]

Examining the bulk term, we see that we are unable to apply the narrow band inequality Lemma 6.8, to the derivative of $\varphi_{h}(\omega)$ and $w_{h}(\omega)$ since the functions only belong to the space $V_{h}\subset H^{1}(D_{0})$, resulting in the bound of order $h$ given in (6.29). However, considering sufficiently regular functions $\varphi,w\in L^{2}(\Omega;H^{2}(D_{0}))$, we are able to employ Lemma 6.8 attaining the estimate of order $h^{2}$ given in (6.30). ∎

The regularity assumption (R1) on the associated dual problem follows again from the symmetry of the bilinear for $a(\cdot,\cdot)$ and the previously derived regularity result given in Theorem 4.4. Hence all the assumptions of the abstract theory are satisfied and we have the stated convergence rate given in Theorem 5.3.

7. Numerical results

In this section, we numerically verify the stated convergence rates of the two proposed finite element discretisations of the reformulated model elliptic problems. In both cases, the numerical scheme has been implemented in DUNE [3, 9].

7.1. Random Surface

As a model for the random surface $\Gamma(\omega)$, we consider a graph-like representation over the unit sphere $\Gamma_{0}=S^{2}$

[TABLE]

where the prescribed height function $h(\omega,\cdot):\Gamma_{0}\rightarrow\mathbb{R}$, will take the form of a truncated spherical harmonic expansion

[TABLE]

with independent, uniformly distributed random coefficients $\lambda_{l,m}\sim U(-1,1)$. Here $\epsilon_{tol}>0$ is a parameter controlling the maximum deviation of the fluctuating surface which in practice will be set to $\epsilon_{tol}=0.1$ and $Y_{l}^{m}$ denotes the spherical harmonic function of degree $l$ and order $m$, which correspond to the eigenvalues of the Laplace-Beltrami operator. For further details on exact form of the spherical harmonics, we refer the reader to [15, 1]. Realisations of the random surface for different samples are given below in Figure 3.

To numerical verify the convergence rate, we set the exact pull-back solution to be given by

[TABLE]

with $\nu_{1},\nu_{2}\sim U(-1,1)$ and $\sigma_{tol}>0$ a constant controlling the largest deviation of pathwise solution. This in turn determines the random data $\hat{f}$ given in the reformulated elliptic equation (4.4). We observe the following errors for the approximation $\mathbb{E}[\hat{u}]-E_{M}[\hat{u}_{h}]$ in $L^{2}(\Omega^{M};L^{2}(\Gamma_{0}))$ and $L^{2}(\Omega^{M};H^{1}(\Gamma_{0}))$ and thus the stated convergence results.

7.2. Random bulk-surface

For the coupled-elliptic system on a random bulk-surface, we adopt a similar approach to the random surface numerical example and prescribe the curved boundary to the random bulk $D(\omega)$ which for simplicity is taken to lie in $\mathbb{R}^{2}$, as a graph

[TABLE]

over the unit circle. Here the random height function will given by a truncated Fourier series

[TABLE]

with independent, uniformly distributed random coefficients $\lambda_{n},\hat{\lambda}_{n}\sim U(-1,1)$. We extend the given boundary process in the normal direction into the interior with a sufficiently smooth blending function to form the stochastic domain mapping

[TABLE]

Here the precise form of the chosen blending function $L_{\delta}(\cdot):\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is given by

[TABLE]

Realisations of the image of the reference domain mappped under the random domain mapping (7.4) are provided in Figure 4.

We set the pull-back of the path-wise bulk solution to be given by

[TABLE]

with uniformly distributed random coefficient $\lambda\sim U(-1,1)$ and $\epsilon_{tol}=0.1$. This determines the pull-back of the path-wise surface solution $\hat{v}$ by the reformulated Robin boundary condition

[TABLE]

from which the data $f$ and $\hat{f}_{\Gamma_{0}}$ can then be computed. Note that in practice, the expectation $\mathbb{E}[\hat{v}]$ and its surface derivative are approximated with Monte-Carlo sampling to sufficiently high accuracy. We observe the following errors and experimental order of convergence for the approximations of the bulk $\mathbb{E}[\hat{u}]-E_{M}[\hat{u}_{h}]$ and the surface $\mathbb{E}[\hat{v}]-E_{M}[\hat{v}_{h}]$ mean solutions.