Method of Green's potentials for elliptic PDEs in domains with random apertures

TL;DR

This paper introduces a boundary integral approach using Green's potentials for elliptic PDEs in domains with random apertures, reducing computational complexity and enabling efficient uncertainty quantification.

Contribution

The paper develops a novel boundary integral method with Green's kernels for elliptic PDEs in domains with random apertures, improving efficiency over traditional discretization methods.

Findings

The method reduces problem dimension by focusing on random apertures.

Multilevel Monte Carlo achieves optimal $ ext{O}( ext{epsilon}^{-2})$ complexity.

Numerical results validate the effectiveness of the proposed approach.

Abstract

Problems with topological uncertainties appear in many fields ranging from nano-device engineering to the design of bridges. In many of such problems, a part of the domains boundaries is subjected to random perturbations making inefficient conventional schemes that rely on discretization of the whole domain. In this paper, we study elliptic PDEs in domains with boundaries comprised of a deterministic part and random apertures, and apply the method of modified potentials with Green's kernels defined on the deterministic part of the domain. This approach allows to reduce the dimension of the original differential problem by reformulating it as a boundary integral equation posed on the random apertures only. The multilevel Monte Carlo method is then applied to this modified integral equation and its optimal $\epsilon^{-2}$ asymptotical complexity is shown. Finally, we provide the qualitative analysis of the proposed technique and support it with numerical results.

Tables

Table 1. Table 1 : Convergence with analytical Green’s kernel.

$N_2$	${\Vert \mu -{\mu }_h\Vert }_{L^{\infty }}$	rate	${\Vert u-u_h\Vert }_{L^{\infty }}$	rate	${\Vert u-u_h\Vert }_{H^1}$	rate
8	5.074e-1	-	3.323e-2	-	9.462e-2	-
20	7.008e-2	1.87	4.506e-3	2.36	2.992e-2	1.56
68	5.215e-3	2.06	1.820e-4	3.33	2.058e-3	2.64
260	3.533e-4	2.01	4.898e-6	2.34	1.312e-5	3.87
1028	2.258e-5	2.00	3.133e-7	2.00	8.176e-7	2.00
4100	1.428e-6	1.99	1.970e-8	2.00	5.140e-8	2.00

Table 2. Table 2 : Convergence with numerical Green’s kernel.

$N_1$	$N_2$	${\Vert \mu -{\mu }_h\Vert }_{L^{\infty }}$	rate	${\Vert u-u_h\Vert }_{L^{\infty }}$	rate	${\Vert u-u_h\Vert }_{H^1}$	rate
48	8	5.071e-1	-	3.323e-2	-	1.048e-1	-
104	20	7.004e-2	1.87	4.507e-3	2.36	3.460e-2	1.53
328	68	5.212e-3	2.07	1.820e-4	3.33	2.462e-3	2.64
1232	260	3.530e-4	2.01	4.901e-6	2.34	1.335e-5	3.98
4848	1028	2.257e-5	2.00	3.135e-7	2.00	8.319e-7	1.99
19336	4100	1.422e-6	1.99	1.971e-8	2.00	5.230e-8	2.00

Table 3. Table 3 : Empirical rates α 𝛼 \displaystyle\alpha of the weak convergence, β 𝛽 \displaystyle\beta of the decay of the variance and ρ 𝜌 \displaystyle\rho of the growth of the cost for boundary integral equations with different kernels.

	Analytical GF	Numerical GF	Schur complement
$\alpha$	$1.88$	$2.02$	$1.89$
$\beta$	$4.33$	$4.31$	$4.42$
$\rho$	$1.58$	$1.97$	$2.52$

Equations

$$-{\nabla }^{2}u(x,\omega )$$

$${\alpha }_{1}u(x,\omega )+{\beta }_1\frac{\partial u(x,\omega )}{\partial n}$$

$${\alpha }_{2}u(x,\omega )+{\beta }_2\frac{\partial u(x,\omega )}{\partial n}$$

$$\displaystyle\displaystyle L^{p}\Big{(}\Omega;$$

$$\displaystyle\displaystyle:=\left\{u:\Omega\to H^{1}(D)\;\Big{|}\begin{array}[]{l}u\text{ is strongly measurable and }\left\|u\right\|_{{L^{p}(\Omega;H^{1}(D))}}<\infty\end{array}\right\}$$

$$\displaystyle\displaystyle\left\|u\right\|_{L^{p}(\Omega;H^{1}(D))}^{p}=\left\{\begin{array}[]{l l}\displaystyle{\int_{\Omega}\left\|u(\cdot,\omega)\right\|_{H^{1}(D)}^{p}d\mathbb{P}(\omega)}&\qquad\text{if }0<p<\infty,\\[10.00002pt] ess\sup_{\omega\in\Omega}\left\|u(\cdot,\omega)\right\|_{H^{1}(D)}&\qquad\text{if }p=\infty.\end{array}\right.$$

$$u(x,\omega )=u_1(x)+u_2(x,\omega ),$$

$$-{\nabla }^2{u}_1(x)$$

$${\alpha }_1{u}_1(x)+{\beta }_1\frac{\partial u_1(x)}{\partial n}$$

$$-{\nabla }^2{u}_2(x,\omega )$$

$${\alpha }_1{u}_2(x,\omega )+{\beta }_1\frac{\partial u_2(x,\omega )}{\partial n}$$

$${\alpha }_2{u}_2(x,\omega )+{\beta }_2\frac{\partial u_2(x,\omega )}{\partial n}$$

$$\varphi (x)=b_2(x)-{\alpha }_2{u}_1(x)-{\beta }_2\frac{\partial u_1(x)}{\partial n}.$$

$$-{\nabla }^2{G}_1(x,\xi )$$

$${\alpha }_1{G}_1(x,\xi )+{\beta }_1\frac{\partial G_1(x,\xi )}{\partial n_x}$$

$$u_2(x,\omega )$$

$$=\int {\text{displaylimits}}_{\partial D_2(\omega )}{G}_1(x,y(\omega )){\nu }_2(y(\omega ))dl(y(\omega ))=\int {\text{displaylimits}}_0^1{G}_1(x,{\eta }_{\omega }(t))\mu ({\eta }_{\omega }(t))dt,$$

$$-\frac{{\beta }_2}{2}\frac{\mu ({\eta }_{\omega }(s))}{\left|{\eta }_{\omega }^{\prime }(s)\right|}$$

$$+{\alpha }_2{\int }_0^1{G}_1({\eta }_{\omega }(s),{\eta }_{\omega }(t))\mu ({\eta }_{\omega }(t))dt=\varphi ({\eta }_{\omega }(s)),$$

$$-\frac{1}{2}\frac{\mu ({\eta }_{\omega }(s))}{\left|{\eta }_{\omega }^{\prime }(s)\right|}+{\int }_0^1\frac{\partial G_1({\eta }_{\omega }(s),{\eta }_{\omega }(t))}{\partial n_{{\eta }_{\omega }(s)}}\mu ({\eta }_{\omega }(t))dt=\varphi ({\eta }_{\omega }(s)),s\in [0;1].$$

$${\int }_0^1{G}_1({\eta }_{\omega }(s),{\eta }_{\omega }(t))\mu ({\eta }_{\omega }(t))dt=\varphi ({\eta }_{\omega }(s)),s\in [0;1],$$

$$(A\mu )(s)={\int }_0^{1}K({\eta }_{\omega }(s),{\eta }_{\omega }(t))\mu ({\eta }_{\omega }(t))dt$$

$$(A\mu )(s)\approx (A_h\mu )(s)=h\sum _{k=0}^{N-1}K({\eta }_{\omega }(s),{\eta }_{\omega }(kh))\mu ({\eta }_{\omega }(kh)),s\in [0;1].$$

$$-\frac{1}{2}\frac{\mu ({\eta }_{\omega }(k^{\prime }h))}{\left|{\eta }_{\omega }^{\prime }(k^{\prime }h)\right|}+h\sum _{k=0}^{N-1}\frac{\partial G_1({\eta }_{\omega }(k^{\prime }h),{\eta }_{\omega }(kh))}{\partial n_{{\eta }_{\omega }(s)}}\mu ({\eta }_{\omega }(kh))=\varphi ({\eta }_{\omega }(k^{\prime }h)),$$

$$k^{\prime }=0,...,N-1.$$

$$(v,w{)}_h=Q_h(vw),$$

$$Q_{h}g=h\sum _{k=0}^{N-1}\sum _{j=1}^J{w}_{j}g((k+{\zeta }_j)h),0<{\zeta }_1<{\zeta }_2<...<{\zeta }_J<1$$

$$\sum _{j=1}^J{w}_j=1,w_j>0,\text{for }1\le j\le J.$$

$$(A_h{\mu }_h,\chi {)}_h=(\varphi ,\chi {)}_h,\forall \chi \in S_h.$$

$$v_k(s)=\{\begin{array}{ll}1-|s-kh|/h,& \text{if }|x-kh|\le h,\\ 0,& \text{otherwise}.\end{array}$$

$$\sum _{k=0}^{N-1}{a}_{l,k}{\mu }_h({\eta }_{\omega }(kh))=(\varphi ,v_l{)}_h,l=0,...,N-1,$$

$$\displaystyle\displaystyle a_{l,k}=h^{2}\sum_{k^{\prime}=0}^{N-1}\sum_{j=1}^{J}w_{j}G_{1}\Big{(}\eta_{\omega}\big{(}(k^{\prime}+\zeta_{j})h\big{)},\eta_{\omega}(kh)\Big{)}{v_{l}}\big{(}(k^{\prime}+\zeta_{j})h\big{)}.$$

$$J$$

$${\zeta }_1$$

$$w_1$$

Full text

††journal: Journal of Scientific Computing

1 Introduction

Scientific and technological development is often linked to the increase in the requirements to the accuracy of mathematical models and numerical methods. As an example, consideration of uncertainties in model inputs and parameters has been attracting a lot of attention of a research community in recent years and a large database of methods for the numerical treatment of boundary-value problems with random coefficients and random input data has been collected. Problems with topological uncertainties have been also studied and their importance has been recognized in many applied fields ranging from nano-device engineering and analysis of micro electromechanical systems (MEMS) Agarwal2007 ; Zhu2004 to the design of bridges Canuto2009 . Other applications include flows over rough surfaces Tartakovsky2006a ; Zayernouri2013 , surface imaging Tsong2005 , corrosion or wear of surfaces, homogenization of random heterogeneous media Savvas2014 and even modelling of blood flow Park2012 .

Existing numerical methods for PDEs in random domains differ by the way of approximating the spatial and random components of the solution. For example, the random solution of a boundary value problem under the small noise assumption can be often represented as a sum of a deterministic component corresponding to the fixed nominal boundary and a small random perturbation which can be quantified using the methods of the “shape calculus” Honda2005 ; Harbrecht2008 .

Alternatively, the original problem in a random domain can be transformed to the problem with random coefficients posed on a deterministic reference domain by means of the random change of variables estimated from a series of auxiliary PDEs Tartakovsky2006b . In conjunction with the stochastic Galerkin method, this approach was considered in Tartakovsky2006a ; Kundu2014 ; Harbrecht2014 while the stochastic collocation approximation in random space was applied in Castrillon2014 ; Castrillon2015 . An equivalent Lagrangian approach was also proposed in Agarwal2007 where the mapping to the reference domain was combined with the stochastic spectral boundary element approximation.

Similarly to the domain mapping method, the random displacement field can be applied directly to the mesh-based representation of the geometry producing the new mesh with random coordinates of nodes but the same fixed connectivity. The main advantage of the mesh-based formulation is that the structure of the underlying linear system remains unchanged enabling reusability of the existing deterministic solvers. This idea was proposed in Mohan2011 in combination with the polynomial chaos approximation in random space and later was studied in Harbrecht2014 in the context of the Quasi-Monte Carlo method applied to the random interface problem.

In fictitious domain methods, the original problem is reformulated on a larger deterministic domain containing all realizations of the random boundary. The enclosing domain can be chosen arbitrarily allowing for simple discretizations which do not have to conform with the random boundaries at the cost of adding new variables to enforce the true boundary conditions. For example, the authors of Nouy2007 ; Nouy2008 enriched the finite element approximation spaces with the suitably constructed functions which allow for the explicit representation of a solution in terms of the random variables describing the geometry. In Savvas2014 ; Lang2013 ; Nouy2010 , this method was successfully applied to a number of problems with stochastic material interfaces. The authors of Canuto2007 satisfied the boundary conditions by introducing the Lagrange multiplier which transformed the original elliptic equation into the larger saddle-point problem. However, the information on the random geometry in the resulting linear system was encoded only in the part of the matrix coupling the primal variable and the boundary supported Lagrange multiplier. Such localization property is a certain advantage of this approach over the domain mapping methods which propagate the boundary uncertainty to the whole domain. Additionally, the method requires no assumptions on the size of the random displacements which is a major limitation of the perturbation techniques.

The aim of this paper is to construct an efficient and accurate numerical method with good localization properties in the sense described above. Our motivation for such formulation is driven by the problems with only certain (often relatively small) part of the boundary being subjected to the random perturbations. In particular, the emphasis of this effort is on elliptic equations in arbitrary (otherwise deterministic) domains with random apertures. The fully discrete formulations of conventional solvers for such problems can lack efficiency due to the necessity in the discretization of the whole physical domain. In this regard, the semi-analitical approximations hold a vast potential. Here we propose to adapt the method of Green’s potentials for elliptic equations to the case of domains with random boundaries. Green’s potentials are the layer potentials with the modified kernels given by the suitably constructed Green’s functions. We will show that the proposed method allows to formulate the original boundary value problem in terms of the integral equations on the random part of the boundary only leading to the potentially significant computational savings. The first step towards the practical application of this approach was done in Melnikov1977 where the so-called method of “modified potentials” was introduced. Later it was successfully applied to both stationary and time-dependent deterministic problems Melnikov2014 ; Melnikov1996 . It is worth noting that the importance of Green’s functions has been already recognized in various areas of the uncertainty quantification Collins2015 ; BarajasSolano2013 ; Neuman1993 ; Neuman1996 ; Manolis1996 . Here we apply the Multilevel Monte Carlo (MLMC) method Giles2008 for the statistical approximation. However, we note that any method of collocation type can be trivially adopted to the proposed numerical technique.

The paper is organized as follows. In section 2, we formulate the problem and introduce the equations and necessary analytical tools. In section 3, we discuss the discretization scheme for the given equations. The complexity analysis of the proposed scheme is given in section 4. Finally, the numerical examples in section 5 are provided in support of the obtained analytical results.

2 Problem setting

Let $\displaystyle(\Omega,\mathcal{F},\mathbb{P})$ be a complete probability space with a set of outcomes $\displaystyle\Omega$, a sigma algebra of events $\displaystyle\mathcal{F}$ and a probability measure $\displaystyle\mathbb{P}$ defined on it. For each outcome $\displaystyle\omega\in\Omega$, define $\displaystyle D(\omega)$ to be the realization of a random domain with a boundary comprised of deterministic and random parts $\displaystyle\partial D(\omega):=\partial D_{1}\cup\partial D_{2}(\omega)$. $\displaystyle\partial D_{1}$ is assumed to be Lipschitz and $\displaystyle\partial D_{2}(\omega)$ is assumed to be analytic. We are concerned with the solutions of the following boundary value problem

[TABLE]

We assume that, for each $\displaystyle\omega$, the solution $\displaystyle u(x,\omega)$ to the above problem exists, is unique and belongs to $\displaystyle H^{1}(D)$, the space of square integrable functions with square integrable first derivatives. Additionally, we require $\displaystyle u({x},\omega)$ to be a Bochner integrable function with values in $\displaystyle H^{1}(D)$, i.e., $\displaystyle u({x},\omega)\in L^{p}(\Omega;H^{1}(D))$, the function space given by

[TABLE]

with the corresponding norm

[TABLE]

For compactness, we will use $\displaystyle L^{p}(\Omega)$ instead of $\displaystyle L^{p}(\Omega;H^{1}(D))$ in later discussions.

Denote by $\displaystyle D_{1}$ the reference deterministic domain containing all realizations of the random boundary $\displaystyle\partial D_{2}(\omega)$. This definition is similar to that used in the fictitious domain methods. However, we explicitly require that the deterministic part of the boundary $\displaystyle\partial D_{1}$ is also the boundary of $\displaystyle D_{1}$. This definition is depicted graphically in Figure 1. Due to linearity of the operators in (2), one can represent the solution $\displaystyle u(x,\omega)$ as a superposition of two functions

[TABLE]

where $\displaystyle u_{1}(x)$ is the deterministic component which satisfies the boundary value problem in the reference domain $\displaystyle D_{1}$

[TABLE]

and the random component $\displaystyle u_{2}(x,\omega)$ can be determined from the following homogeneous boundary value problem

[TABLE]

with the boundary condition on $\displaystyle\partial D_{2}(\omega)$ defined by the trace of $\displaystyle u_{1}(x)$ as follows

[TABLE]

2.1 Method of Green’s potentials

The particularly simple structure of the problem in (2) is very well suited for the construction of efficient solvers. One of such solvers, namely the method of Green’s potentials, is proposed in this section.

Firstly, define the Green’s function corresponding to the boundary value problem (6) as a solution of the complementary problem

[TABLE]

where $\displaystyle\delta(\xi)$ is the Dirac measure of unit mass at point $\displaystyle\xi$. With this function at hand, the solution $\displaystyle u_{2}(x,\omega)$ of the problem in (2) allows the representation in the form of the single-layer potential Jaswon1963 ; Symm1963

[TABLE]

where $\displaystyle\mu(\eta_{\omega}(t))=\nu_{2}(\eta_{\omega}(t))\left|\eta_{\omega}^{\prime}(t)\right|$ and $\displaystyle\eta_{\omega}(t)$ defines a smooth parameterization of the random boundary curve. Note that in the above expression, the first integral over the deterministic part of the boundary vanishes because of the Green’s function $\displaystyle G_{1}(x,\xi)$.

Using (8) and the jump conditions of the derivative of the single-layer potential on the boundary, the unknown density $\displaystyle\mu(\eta_{\omega}(t))$ of the potential in (2.1) can be obtained from the following boundary integral equation

[TABLE]

which is a Fredholm equation of the second kind. The pure Neumann problem ($\displaystyle\alpha_{2}=0$, $\displaystyle\beta_{2}=1$) yields the similar equation

[TABLE]

It is well known that equations of the second kind are well-posed. On the other side, Dirichlet boundary conditions ($\displaystyle\alpha_{2}=1$, $\displaystyle\beta_{2}=0$) convert (12) into the Fredholm equation of the first kind

[TABLE]

which is intrinsically ill-posed and must be treated with a special care. It has been established that equations of the first kind with logarithmically singular kernels admit unique solutions when the conformal radius of the boundary is not equal to one Yan1988 ; Cheng1993 . Therefore, we assume that all boundaries satisfy this condition.

Note that the traditional approach in solving the Dirichlet problems with methods of potential is based on representing the solution in the form of the double-layer potential which results in equations of the second kind. The use of the single-layer potentials, however, has the advantage of satisfying the governing equation on the deterministic boundary. As a result, when the length of the boundary $\displaystyle\partial D_{1}$ is large, the method of Green’s potentials can lead to significant computational savings compared to the traditional approaches relying on the discretization of the whole boundary.

Of course, the efficiency of the method of Green’s potentials relies on the availability of Green’s functions for the specific geometries of the domain. Unfortunately, analytical expressions for the Green’s functions are known only for very simple domains and the construction of approximations adds an additional level of complexity to the proposed scheme. However, in the case of uncertain domains, the deterministic complementary problem (9)-(10) has to be solved only once and the value of the Green’s function at any field point is then readily available through the simple matrix-vector product which can be done very efficiently. The implementation aspects of this approach are discussed in succeeding sections.

3 Discretization scheme

3.1 Spatial discretization

The boundary integral equation (13) has been extensively studied in the literature as the classical equation of potential theory Atkinson1991 ; Atkinson1997 ; GRAHAM1985 ; Jeon1997 ; Sloan2000 ; Cheng1993 ; Cheng1995 . It is an equation of the second kind with a continuous kernel and any approximation technique is expected to work well. For the sake of completeness, we will use the classical Nyström method. We start with the boundary integral operator

[TABLE]

with a kernel $\displaystyle K$ and approximate it with the trapezoidal rule on the uniform grid with the step $\displaystyle h~{}=~{}1/N$ for some integer $\displaystyle N$

[TABLE]

After collocating the resulting discrete operator at the quadrature nodes this yields the following numerical scheme for the integral equation in (13)

[TABLE]

The order of this scheme is determined by the order of the trapezoidal quadrature rule in (15), i.e., it is at least $\displaystyle O(h^{2})$ for sufficiently smooth $\displaystyle\mu$ Atkinson1997 .

The direct extension of the scheme in (16) to the equations of the first kind is not possible due to the singularity of the kernel. In this case, we will use the quadrature technique proposed in Sloan1992 for the first kind Fredholm equations with logarithmic kernels on closed curves. It is a fully discrete method of qualocation type based on the composite quadrature rule, i.e., both the integral operator and the Galerkin projection are approximated with suitable quadratures. In particular, we use the trapezoidal approximation of the integral operator in (15) and then project it on the test space $\displaystyle S_{h}$ of 1-periodic smoothest splines of order $\displaystyle r$ with the discrete inner product

[TABLE]

where

[TABLE]

and

[TABLE]

The problem now can be formulated as follows: find $\displaystyle\mu_{h}$ such that

[TABLE]

For $\displaystyle r=2$, the basis $\displaystyle(v_{0},...,v_{N-1})$ of $\displaystyle S_{h}$ is represented by the classical hat functions

[TABLE]

Given the basis, one can write the discrete formulation of the problem: find $\displaystyle\mu_{h}$ such that

[TABLE]

where

[TABLE]

It was shown in Sloan1992 that the following unique choices of the quadrature rules are optimal in terms of stability of the approximation

[TABLE]

The choice of nodes in (3.1’) gives the $\displaystyle O(h^{3})$ order of uniform convergence. The scheme with nodes in (3.1) has $\displaystyle O(h^{2})$ accuracy but the linear system in (18) has the simpler form

[TABLE]

The discrete systems in (16) and (20) have similar form and we use them to solve equations of the second and first kind respectively. The combined linear system for the problems with both Dirichlet and Neumann boundary conditions can be written as

[TABLE]

where the elements $\displaystyle\mathbf{G}_{k^{\prime},k}$ of the matrix $\displaystyle\mathbf{G}$ are given by

[TABLE]

$\displaystyle\mathcal{I}_{k^{\prime}-k}$ is the indicator function and $\displaystyle\mathcal{D}$, $\displaystyle\mathcal{N}$ denote the collections of nodes with Dirichlet and Neumann boundary conditions respectively. The quadrature and collocation points for the case of the Neumann and Dirichlet boundary conditions are depicted in Figure 2.

After the density of the potential is determined from (21), one can calculate the solution $\displaystyle u_{2}(x,\omega)$ at any field point by evaluating the integral with the trapezoidal quadrature rule (15) as

[TABLE]

3.2 Evaluation of Green’s functions for arbitrary domains

Recall the definition of the Green’s function for the Laplace operator as a solution of the following boundary value problem

[TABLE]

It was mentioned previously that the proposed numerical technique relies heavily on the ability to evaluate Green’s functions for the domains of arbitrary shapes. We outline here several methods which allow to solve (9)-(10) in a computationally attractive way.

3.2.1 Analytical Green’s functions.

In certain cases, it is possible to solve (9)-(10) analytically. The fundamental solution of the 2-D Laplace equation is one of such examples of an exceptional importance. It satisfies the equation (9) in the entire space and has the form

[TABLE]

The Green’s function to (9)-(10) can then be written as a sum of (22) and a “corrector” function aiming to satisfy the boundary condition in (10)

[TABLE]

Defined in this way, the corrector function is called the regular component of the Green’s function. It solves the following complementary problem

[TABLE]

For simple geometries possessing certain symmetry properties, the problem in (24)-(25) admits the closed form solution which can be constructed via the method of images. Two classical examples of such Green’s functions are given below.

Dirichlet problem in the upper half-plane $\displaystyle D_{1}(x_{1},x_{2})=\{x_{1},x_{2}\geq 0\}$

[TABLE]

Dirichlet problem in the disk $\displaystyle D_{1}(r,\varphi)=\{0\leq r<a,0\leq\varphi\leq 2\pi\}$

[TABLE]

More examples can be found, for instance, in Duffy2001 ; Melnikov2011 . Additionally, the infinite product representation of Green’s functions arising from applying the method of images was discussed in Melnikov2012 ; Reshniak2013 .

3.2.2 Direct numerical approximation of Green’s functions.

By its definition, Green’s function of the boundary value problem is the inverse of the corresponding differential operator. Numerical methods can be viewed as implicitly approximating such inverse operators. For instance, it was shown in Tottenham1970 that the finite element solution $\displaystyle u_{h}(x)$ of the boundary value problem with homogeneous boundary conditions has the form

[TABLE]

where $\displaystyle G_{h}(x,\xi)$ is the FE-Green’s function, i.e., the projection of the exact Green’s function on the finite element space $\displaystyle V_{h}$. Solving for $\displaystyle G_{h}(x,\xi)$ yields

Theorem 1** (Hartmann2013 ).**

Let $\displaystyle\mathbf{K}$ be the stiffness matrix of the linear system arising after the finite element disretization. The FE-Green’s function has the form

[TABLE]

where $\displaystyle\mathbf{v}(x)=(v_{1}(x),...,v_{M}(x))$ are the basis functions of the FE-space $\displaystyle V_{h}$.

For the spectral finite element method, the formula (26) simplifies to

[TABLE]

where $\displaystyle(v_{1},...,v_{M})\in V_{h}$ are the $\displaystyle M$ leading eigenfunctions of the differential operator and $\displaystyle\mathbf{\Lambda}$ is the diagonal matrix with the corresponding eigenvalues. One can immediately recognize in (27) the truncation of the classical eigenfunction representation of the Green’s function which is guaranteed to exist by the Mercer’s theorem Fasshauer2016 . Spectral representations of the Green’s functions can be found, for instance, in Duffy2001 . As an example, we provide the Green’s function for the

Dirichlet problem in the rectangle $\displaystyle D_{1}(x_{1},x_{2})=\{0\leq x_{1}\leq a,0\leq x_{2}\leq b\}$

[TABLE]

3.2.3 Numerical approximation of the regular part of Green’s functions.

Green’s functions generally do not belong to the function spaces approximated by the $\displaystyle span(v_{1}(x),...,v_{M}(x))$ resulting in the very slow convergence of the representations in (26)-(27). For instance, solutions to the Poisson equation are usually constructed in $\displaystyle H^{1}(D_{1})$ but the solution to the problem in (9)-(10) is not in $\displaystyle H^{1}(D_{1})$ since the delta function $\displaystyle\delta(\xi)\notin H^{-1}$ for $\displaystyle d\geq 2$, where $\displaystyle d$ is the physical dimension of the problem. Analogously, analytical expansions like (28) do not have uniform error estimates which seriously limits their immediate practical utilization.

In certain cases, one can obtain uniformly convergent spectral representations by partial summation of the series leading to the explicit extraction of the singularity. Justification of this approach with practical examples can be found in Melnikov1998 . For instance, the series in (28) is transformed to the following form

[TABLE]

where $\displaystyle\nu=n\pi/b$, $\displaystyle z=x_{1}+ix_{2}$, $\displaystyle z_{1}=(x_{1}+a)+ix_{2}$, $\displaystyle z_{2}=(x_{1}-a)+ix_{2}$, $\displaystyle\zeta=\xi_{1}+i\xi_{2}$, $\displaystyle\zeta_{1}=(\xi_{1}+a)+i\xi_{2}$, $\displaystyle\zeta_{2}=(\xi_{1}-a)+i\xi_{2}$, $\displaystyle\zeta_{1}^{*}=(\xi_{1}+a)-i\xi_{2}$, $\displaystyle\zeta_{2}^{*}=(\xi_{1}-a)-i\xi_{2}$, $\displaystyle E(z)=\left|e^{\pi z/b}-1\right|$ and

[TABLE]

The remainder term $\displaystyle R_{M}(x,\xi)$ of the $\displaystyle M$-term truncation of the expansion in (29) has the upper bound

[TABLE]

which reveals the extremely high rate of convergence.

Similarly, the corrector function $\displaystyle\psi(x,\xi)$ in (23) is harmonic everywhere in $\displaystyle D_{1}$ and thus can be efficiently approximated with any conventional numerical method. For instance, the finite element approximation has the form

[TABLE]

where $\displaystyle\mathbf{K}$ is the same stiffness matrix as in (26) and $\displaystyle\boldsymbol{\mathfrak{g}}(\xi)$ encodes the negative trace of the fundamental solution on the boundary.

One can also construct the regular part of the Green’s function in the form of the single layer potential

[TABLE]

where $\displaystyle\mathfrak{G}(x,y)$ is the fundamental solution of the differential operator, $\displaystyle\mu^{\psi}(\eta_{1}(t),\xi)=\nu_{1}(\eta_{1}(t),\xi)\left|\eta_{1}^{\prime}(t)\right|$ and $\displaystyle\eta_{1}(t)$ defines a parameterization of the boundary $\displaystyle\partial D_{1}$. It is natural to build the approximate solution of the above equation with the same method used for the approximation of the original integral equation, e.g., with the scheme given in section 3.1. For this, note that similarly to (30), the density $\displaystyle\mu^{\psi}(\eta_{1}(t),\xi)$ encodes the negative trace of the fundamental solution $\displaystyle\mathfrak{G}(\cdot,\xi)$ on the boundary $\displaystyle\partial D_{1}$.

Denote by $\displaystyle\boldsymbol{\mathfrak{G}}$ the marix of the linear system in (21) for the problem formulated over the whole boundary $\displaystyle D=D_{1}\cup D_{2}$ and using fundamental solution $\displaystyle\mathfrak{G}(x,\xi)$ instead of the Green’s function $\displaystyle G_{1}(x,\xi)$. By taking \displaystyle\boldsymbol{\mu}^{\psi}(\xi)=\big{(}\mu^{\psi}_{0}(\xi),...,\mu^{\psi}_{N-1}(\xi)\big{)} with $\displaystyle\mu^{\psi}_{k}(\xi)=\mu^{\psi}(\eta_{1}(kh),\xi)$, we get

[TABLE]

where $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$ is a $\displaystyle N_{1}\times N_{1}$ submatrix of $\displaystyle\boldsymbol{\mathfrak{G}}$ with both $\displaystyle\eta_{\omega}(k^{\prime}h)\in D_{1}$ and $\displaystyle\eta_{\omega}(kh)\in D_{1}$. The matrix $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$ needs to be inverted only once and the approximation of the regular part of the Green’s function at any point is then readily available as

[TABLE]

It is easy to verify that in this case the matrix $\displaystyle\mathbf{G}\in\mathbb{R}^{N_{2}\times N_{2}}$ of the original linear system in (21) takes the form of the Schur complement of $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$

[TABLE]

where by analogy with $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$, $\displaystyle\boldsymbol{\mathfrak{G}}_{22}$ is a $\displaystyle N_{2}\times N_{2}$ submatrix of $\displaystyle\boldsymbol{\mathfrak{G}}$ with $\displaystyle\eta_{\omega}(k^{\prime}h)\in D_{2}$ and $\displaystyle\eta_{\omega}(kh)\in D_{2}$, and similarly for $\displaystyle\boldsymbol{\mathfrak{G}}_{21}\in\mathbb{R}^{N_{2}\times N_{1}}$ and $\displaystyle\boldsymbol{\mathfrak{G}}_{12}\in\mathbb{R}^{N_{1}\times N_{2}}$. We will use this approach in the subsequent sections. Two examples of the approximate Green’s function obtained in this way are given in Figure 3.

3.3 Statistical discretization

Let $\displaystyle u(x,\omega)$ be the solution of the random partial differential equation and denote by $\displaystyle f(\omega):=f(u(x,\omega))$ some functional of $\displaystyle u(x,\omega)$. The expected value of $\displaystyle f(\omega)$ can be approximated by the Monte-Carlo (MC) estimator of the form

[TABLE]

where the deterministic functions $\displaystyle f^{m}$ represent i.i.d. realizations of $\displaystyle f(\omega)$.

Monte Carlo method is a purely statistical technique which ignores any information about regularity of the functions in the physical space. The multilevel Monte Carlo method can exploit some of this information by relating the sampling error of the estimator to the convergence properties of the spatial discretization Giles2008 . For instance, consider a hierarchical family of nested discretizations with step sizes

[TABLE]

where $\displaystyle q\in\mathbb{N}\setminus 1$ is a refinement parameter. Denote by $\displaystyle u_{l}(x,\omega)$ the approximation of $\displaystyle u(x,\omega)$ at the level $\displaystyle l$ and let $\displaystyle f_{l}(\omega)=f(u_{l}({x},\omega))$. Then the solution at the finest discretization level $\displaystyle L$ is given as the telescoping sum

[TABLE]

Taking advantage of the linearity of the expectation, we obtain

[TABLE]

By setting $\displaystyle\Delta_{l}(\omega)=f_{l}(\omega)-f_{l-1}(\omega),\;l=0,...,L,$ the above expression yields the multilevel Monte-Carlo estimator

[TABLE]

where $\displaystyle\Delta_{0}(\omega)=f_{0}(\omega)$ and $\displaystyle M_{l}$ is the number of random samples at the level $\displaystyle l$.

4 Complexity analysis

4.1 Asymptotical complexity of the MLMC estimator

Consider the cost function of the MLMC estimator

[TABLE]

where $\displaystyle M_{l}$ is the number of samples at the level $\displaystyle l$ and $\displaystyle C_{l}$ is the cost of generating the single realization of $\displaystyle f_{l}(\omega)$. The following theorem shows that the above cost can be optimized by appropriately balancing the errors of the Monte Carlo estimators at different levels.

Theorem 2** (Giles2008 ).**

If there exist independent estimators $\displaystyle\mathbb{E_{MC}}\left[{\Delta_{l}}\right]$ based on $\displaystyle M_{l}$ Monte Carlo samples, each with expected cost $\displaystyle C_{l}$ and variance $\displaystyle V_{l}$, and positive constants $\displaystyle\alpha$, $\displaystyle\beta$, $\displaystyle\rho$ such that $\displaystyle\min(\beta,\rho)\leq 2\alpha$ and

\displaystyle\Big{|}\mathbb{E}\left[{f_{l}(\omega)-f\big{(}u(x,\omega\big{)}}\right]\Big{|}\lesssim h_{l}^{\alpha}, 2. 2.

$\displaystyle V_{l}\lesssim h_{l}^{\beta}$, 3. 3.

$\displaystyle C_{l}\lesssim h_{l}^{-\rho}$,

then for any $\displaystyle\epsilon<e^{-1}$ there are values $\displaystyle L$ and $\displaystyle M_{l}$ for which the multilevel estimator (33) has a mean-square-error with bound

[TABLE]

with a computational complexity $\displaystyle C_{\mathbb{ML}}$ with bound

[TABLE]

Corollary. According to Theorem 2, the optimal $\displaystyle\epsilon^{-2}$ complexity of the MLMC estimator is achieved when $\displaystyle\rho<\beta$, i.e., when the decay of the variance $\displaystyle V_{l}$ of the level corrections is faster than the growth of the computational cost $\displaystyle C_{l}$. For Lipschitz continuous $\displaystyle f$ and given the convergence rate $\displaystyle\alpha$ of the spatial discretization scheme, the variance $\displaystyle V_{l}$ decays with order $\displaystyle\beta=2\alpha$. In practice, $\displaystyle\rho\in(2,3]$ for different methods yielding $\displaystyle\alpha>\frac{\rho}{2}\in(1,1.5]$ as an optimality condition which is easily satisfied with most practical schemes including the one utilized in this paper. However, the complexity estimate in (34) has a silent constant which determines the actual relative performace of different algorithms with the same asymptotic complexity. In the following sections, we provide the complexity analysis for the proposed numerical scheme and use it to justify our choice of the method of Green’s potentials.

4.2 Error component analysis

Consider the following decomposition of the total error of the MLMC estimator

[TABLE]

where $\displaystyle u$ and $\displaystyle\tilde{u}_{L}$ are the exact and the approximate values of the potential (2.1). The error components $\displaystyle I$ and $\displaystyle II$ correspond to the spatial approximation error and the sampling error. To achieve the desired accuracy $\displaystyle\epsilon$, it is sufficient to balance the total error between these two components in the following way

[TABLE]

The error analysis for each of the components is provided below.

4.2.1 Spatial discretization error.

Jensen’s inequality gives

[TABLE]

and the estimate of the error component $\displaystyle I$ can be derived from the convergence properties of the spatial discretization scheme. The order of such scheme is in turn controlled by the convergence properties of the boundary integral equation solver and the numerical error of the single-layer potential evaluation.

For instance, it was shown in Sloan1992 that the qualocation-type discretization (20) of the first-kind integral equation in (14) admits the following estimate for the approximation error of $\displaystyle\mu\in H^{t}([0,1])$

[TABLE]

provided that $\displaystyle s>\frac{1}{2}$, $\displaystyle s+\frac{1}{2}<t\leq s+\alpha$ and the right hand side of the integral equation is continuous and 1-periodic. In the case of optimal regularity, i.e., for $\displaystyle\mu\in H^{t}([0,1])$ with $\displaystyle t>\alpha+1/2$, the following error bound is valid

[TABLE]

due to the embedding of $\displaystyle H^{s}$ ($\displaystyle s>1/2$) in $\displaystyle C_{p}$, the space of 1-periodic continuous functions. The order of convergence is $\displaystyle\alpha=2$ for the scheme in (20) and hence the method is $\displaystyle O(h_{l}^{2})$ accurate.

The error of the Nyström method in (16) for the integral equations of the second kind is given by the error of the corresponding numerical quadrature rule. Since the trapezoidal rule is exteremely well-behaved for smooth periodic functions, the order of the Nyström-trapezoidal method is also at least $\displaystyle O(h_{l}^{2})$ Atkinson1997 .

To determine the overall order of the spatial discretization at the level $\displaystyle l$, consider the error in approximating the single-layer potential (2.1)

[TABLE]

where the first term on the right side gives the error of the numerical scheme with the exact Green’s function and the second term is the error due to the approximation of the Green’s kernel itself. We get for the first component that

[TABLE]

where $\displaystyle R(x)$ is the error of the trapezoidal rule in (15). For a periodic function $\displaystyle\mu\in H^{t}$, $\displaystyle R(x)$ is also $\displaystyle O(h^{t})$ (Sloan1992, , Lemma A1), and since $\displaystyle t>\alpha$ in (39), one gets

[TABLE]

Similarly, the Green’s kernel is a harmonic and thus analytic function at any internal point of the domain which gives the estimate for the error in the approximation of the derivatives of the potential

[TABLE]

where $\displaystyle i=(i_{1},i_{2})$ is a multi-index and $\displaystyle f^{(i)}=\displaystyle{\frac{\partial^{|i|}f}{\partial x_{1}^{i_{1}}\partial x_{2}^{i_{2}}}}$.

Now consider the approximation of the Green’s kernel. From (23) and (31), we have

[TABLE]

By analogy with (40), there holds the error estimate

[TABLE]

By combining (39), (40) and (42), one gets the error bound

[TABLE]

The estimate for the first error component follows trivially from (41) as

[TABLE]

Hence, the condition in (37) on the spatial discretization error can be satisfied by choosing the number of levels according to

[TABLE]

Remark. We used the embedding argument in (39), the analyticity of the Green’s kernel away from the boundary and the high accuracy of the trapezoidal rule for smooth periodic functions to find the errors (40)-(41) of the PDE solution $\displaystyle u(x)$ and its derivatives $\displaystyle u^{(i)}(x)$. For Galerkin methods, similar estimates can be determined directly from the energy errors of the form (38) using the Aubin-Nitsche-duality argument sauter2010boundary ; steinbach2007numerical . Any such method to solve boundary integral equations can be trivially applied instead of the one utilized in this paper which might be of interest for the problems in higher dimensions or in regions with piecewise smooth boundaries. In this case, the results of the MLMC theory remain valid but with possibly different values for the rates $\displaystyle\alpha$, $\displaystyle\beta$, $\displaystyle\rho$ in conditions $\displaystyle(1)-(3)$ of Theorem 2. For planar regions with analytic contours, the high accuracy of the quadrature based methods combined with the trivial construction of their matrices motivated our choice of the solvers in (16)-(20).

4.2.2 Sampling error.

From the definition of the norm in (4), we obtain

[TABLE]

By virtue of the independence of i.i.d. samples $\displaystyle u^{m}$, we have \displaystyle\mathbb{C}ov\Big{(}u^{(i),m},u^{(i),m^{\prime}}\Big{)}~{}=~{}0 and

[TABLE]

where $\displaystyle\overline{V}(u)=\mathbb{E}\left[{\left\|u-\mathbb{E}\left[{u}\right]\right\|_{H^{1}(D)}^{2}}\right]$.

Due to the independence of the MC estimators at each level, one gets the error of the MLMC estimator as follows

[TABLE]

where $\displaystyle\overline{V}_{l}=\overline{V}(\Delta_{l})$ and $\displaystyle\Delta_{l}=u_{l}-u_{{l-1}}$.

Thus, the sampling error can be estimated as

[TABLE]

The above condition can be satisfied by taking

[TABLE]

where the coefficients $\displaystyle a_{l}$ are the weights assigning certain part of the sampling error to each level. It was shown in Giles2008 that the optimal cost in Theorem 2 is achieved with the following choice

[TABLE]

Note that, in view of (43), one has $\displaystyle\overline{V}_{l}=O(h_{l}^{2\alpha})$, i.e., $\displaystyle\beta=2\alpha\geq 4$ in Theorem 2.

4.2.3 Complexity of the method of Green’s potentials with different kernels

In this section, we study and compare the complexity of the method of Green’s potentials for three different choices of the kernel function: fundamental solution, analytical Green’s kernel and approximate Green’s kernel.

4.2.4 Fundamental solution

Here we provide the complexity analysis of the scheme in section 3.1 with the boundary integral equations formulated on the whole boundary $\displaystyle\partial D=\partial D_{1}~{}\cup~{}\partial D_{2}$. Consider the overall cost of the MLMC estimator with $\displaystyle M_{l}$ as in (46)

[TABLE]

The costs $\displaystyle C_{l}$ of the BIE solver consist of the two main components

the cost $\displaystyle C^{a}_{l}$ of assembling the matrix of the linear system, 2. 2.

the cost $\displaystyle C^{s}_{l}$ of solving the linear system.

With the appropriate enumeration of the degrees of freedom, one can write the linear system in (21) as

[TABLE]

where the block $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$ involves only points on the fixed boundary $\displaystyle\partial D_{1}$ and the remaining blocks depend on realizations of the boundary $\displaystyle\partial D_{2}(\omega)$. Using the Schur complement, the above system can be reduced to the simpler one

[TABLE]

Since $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$ is fixed, it should be inverted only once using, e.g., LU decomposition. Hence the cost of evaluating $\displaystyle\boldsymbol{\mathfrak{G}}_{11}^{-1}$ can be neglected.

Let $\displaystyle N_{1,l}$ and $\displaystyle N_{2,l}$ be the numbers of degrees of freedom corresponding to the fixed and random parts of the boundary $\displaystyle\partial D_{1}$ and $\displaystyle\partial D_{2}$ at the level $\displaystyle l\in[0,L]$. Then the total number of degrees of freedom is $\displaystyle N_{l}=N_{1,l}+N_{2,l}$. The construction of matrices $\displaystyle\boldsymbol{\mathfrak{G}}_{22}$, $\displaystyle\boldsymbol{\mathfrak{G}}_{21}$ and $\displaystyle\boldsymbol{\mathfrak{G}}_{12}$ requires $\displaystyle O(N_{l}^{2}-N_{1,l}^{2})$ evaluation operations. Evaluation of the Schur complement $\displaystyle\boldsymbol{\mathfrak{G}}_{22}-\boldsymbol{\mathfrak{G}}_{21}\boldsymbol{\mathfrak{G}}_{11}^{-1}\boldsymbol{\mathfrak{G}}_{12}$ requires $\displaystyle O(N_{1,l}N_{2,l}N_{l}+N_{2,l}^{2})$ arithmetic operations assuming standard matrix multiplication algorithm and that $\displaystyle\boldsymbol{\mathfrak{G}}_{11}^{-1}$ is given. Evaluation of the right hand side requires additional $\displaystyle O(N_{1,l}N_{2,l}+N_{2,l})$ arithmetic operations since $\displaystyle\boldsymbol{\mathfrak{G}}_{11}^{-1}f_{1}$ is fixed. Hence the cost $\displaystyle C_{l}^{a}$ of assembling the matrix involves $\displaystyle O(N_{1,l}N_{2,l}N_{l}+N_{2,l}^{2}+N_{1,l}N_{2,l}+N_{2,l})$ arithmetic operations and $\displaystyle O(N_{l}^{2}-N_{1,l}^{2})$ evaluation operations.

For direct linear solvers, the cost of solving the resulting linear system is equal to $\displaystyle O(N_{2,l}^{\gamma})$ for some $\displaystyle\gamma\in(2,3]$; we assume that $\displaystyle\gamma=3$. Evaluation of $\displaystyle\mu_{1}$ requires only $\displaystyle O(N_{1,l}N_{2,l}+N_{1,l})$ arithmetic operations since $\displaystyle\boldsymbol{\mathfrak{G}}_{11}^{-1}f_{1}$ and $\displaystyle\boldsymbol{\mathfrak{G}}_{11}^{-1}\boldsymbol{\mathfrak{G}}_{12}(\omega)$ are available from the first equation. Therefore, the overall cost of solving the linear system can be estimated as $\displaystyle C_{l}^{s}=O(N_{2,l}^{3}+N_{1,l}N_{2,l}+N_{1,l})$.

Omitting evaluation operations, one may conclude that the total cost at each level $\displaystyle l\in[0,L]$ behaves asymptotically as

[TABLE]

since $\displaystyle N_{l}\simeq|\partial D|h_{l}^{-1}$.

Taking into account (49) and since $\displaystyle\rho=3$, $\displaystyle\beta=2\alpha\geq 4$, the $\displaystyle\epsilon$-cost of the MLMC estimator follows from Theorem 2 as

[TABLE]

4.2.5 Analytical Green’s kernel.

For analytical Green’s kernel, the cost $\displaystyle C_{l}^{a}$ of assembling the matrix involves only $\displaystyle O(N_{2,l}^{2})$ evaluation operations and the cost of solving the linear system is $\displaystyle C_{l}^{s}=O(N_{2,l}^{3})$. Omitting evaluation operations, the total asymptotical cost at each level $\displaystyle l\in[0,L]$ is just

[TABLE]

Hence $\displaystyle\rho=3$ and the $\displaystyle\epsilon$-cost of the MLMC estimator follows from Theorem 2 as

[TABLE]

4.2.6 Approximate Green’s kernel.

Both (32) and (47) have the Schur complement of $\displaystyle\boldsymbol{\mathfrak{G}}_{11}$ as their matrices. Hence, the cost $\displaystyle C_{l}^{a}$ of assembling the matrix involves the same $\displaystyle O(N_{1,l}N_{2,l}N_{l}+N_{2,l}^{2})$ arithmetic operations. However, evaluation of the right hand side does not require additional operations. The cost of solving the linear system is again $\displaystyle C_{l}^{s}=O(N_{2,l}^{3})$, and the total cost at each level has the following asymptotical behavior

[TABLE]

The overall asymptotical complexity of the MLMC estimator is thus given by

[TABLE]

Corollary. It is clear from (50), (52) and (54) that any choice of the kernel function results in the optimal asymptotical complexity of the MLMC estimator. However, the direct comparison of (49), (51) and (53) reveals that the classical method requires extra $\displaystyle O(N_{1,l}N_{2,l}N_{l}+N_{2,l}^{2}+N_{1,l}N_{2,l}+N_{l})$ and $\displaystyle O(N_{1,l}N_{2,l}+N_{l})$ operations at each level when compared to the methods with analytical and numerical Green’s kernels respectively. In the next section, we provide the numerical assessment of the impact of the choice of the Green’s kernel on the performance of the proposed algorithm.

5 Numerical results

Example 1.

In the first example, we test the accuracy of the proposed discretization scheme for the fixed deterministic boundary. Consider the problem

[TABLE]

where $\displaystyle D$ is the square domain with a single aperture. Figure 4 shows the discretization of the boundary for the case of analytical and numerical Green’s kernels. For BIE with exact kernel, only the boundary $\displaystyle\partial D_{2}$ of the aperture has to be discretized while for BIE with approximate kernel, it is necessary to discretize both boundaries. $\displaystyle G_{1}(x,\xi)$ is the Green’s function (29) for the square bounded by $\displaystyle\partial D_{1}$ with the source located at the center of the aperture. This choice of the boundary condition on $\displaystyle\partial D_{2}$ suggests $\displaystyle G_{1}(x,\xi)$ as the analytical solution of the above problem.

Convergence properties of the numerical scheme (18) with the quadrature nodes in (3.1) are given in Tables 1 and 2. The reference density of the potential $\displaystyle\mu$ was evaluated with the higher order scheme using the quadrature nodes in (3.1’). The apparent rates of convergence are defined as $\displaystyle\alpha_{h}=\displaystyle{\frac{\log(e_{l}/e_{l-1})}{\log(h_{l}/h_{l-1})}}$. It is seen that the errors have the order of convergence $\displaystyle\alpha=2$ in all norms as predicted by analysis. Obviously, the errors in Table 2 have larger values but the difference is not large and the order is not reduced.

Results in Tables 1 and 2 are also presented graphically in Figure 5 which illustrates the supremum norm of the error along the isocontours of the boundary. Superiority of the analytical Green’s kernel is obvious near the deterministic boundary $\displaystyle\partial D_{1}$ but both approaches show good results far from the boundaries.

Figure 6 illustrates the costs $\displaystyle C^{a}$ of assembling the matrix and $\displaystyle C^{s}$ of solving the resulting system. As predicted, the cost $\displaystyle C^{s}$ is asymptotically dominant for BIE with analytical Green’s kernel while the cost $\displaystyle C^{a}$ dominates in the case of numerical Green’s function. Figure 7 depicts the overall empirical computational costs of different methods along with the speedup of the proposed schemes over the standard BIE method with fundamental solution as kernel. It is seen that methods with both analytical and numerical Green’s kernels are superior to standard approach and perform exceptionally well for small ratios of $\displaystyle\frac{|\partial D_{2}|}{|\partial D_{1}|}$.

Example 2.

For the second example, consider the problem

[TABLE]

where $\displaystyle D$ and $\displaystyle\partial D_{1}$ are the same as in Example 1 and $\displaystyle\partial D_{2}$ is the aperture with the following parametrization

[TABLE]

The radius of the aperture is defined as

[TABLE]

with the mean radius $\displaystyle\overline{R}(t)$ and the random coefficients $\displaystyle a_{n}(\omega)=\mathcal{U}(-\sqrt{3},\sqrt{3})$, $\displaystyle b_{n}(\omega)=\mathcal{U}(-\sqrt{3},\sqrt{3})$. Coefficient $\displaystyle\sigma_{r}$ controls intensity of the random perturbation. Coordinates of the center are also random variables

[TABLE]

where $\displaystyle(\overline{x}_{c},\overline{y}_{c})$ is the mean location of the center and the coefficients $\displaystyle\sigma_{x}$, $\displaystyle\sigma_{y}$ control deviation from the mean. For this particular example, we set $\displaystyle(\overline{x}_{c},\overline{y}_{c})=(0.3,0.4)$, $\displaystyle\sigma_{x}=\sigma_{y}=0.05$, $\displaystyle\overline{R}=0.15$, $\displaystyle\sigma_{r}=0.01$ and $\displaystyle s=10$.

The forcing term $\displaystyle f(x)$ is chosen such that the analytical solution in the deterministic domain $\displaystyle D_{1}$ without the aperture is given by

[TABLE]

Figure 8 illustrates five different realizations of the random geometry and isolines of the solutions corresponding to two particular realizations.

Consider the functional

[TABLE]

where $\displaystyle C$ is the contour parallel to the aperture with the offset $\displaystyle d=0.01$. A single realization of the geometry and the corresponding contour $\displaystyle C$ is depicted on Figure 9 along with 200 realizations of the value of the functional.

The observed rates $\displaystyle\alpha$, $\displaystyle\beta$, $\displaystyle\rho$ in Theorem 2 for the proposed numerical scheme are given in Table 3. It is seen that, due to the second order of convergence of the spatial approximation, the variance decay is much faster than the growth of the cost for all types of kernels. Therefore, the theoretical $\displaystyle\epsilon$-complexity of the MLMC estimator is proportional to $\displaystyle\epsilon^{-2}$ and Figure 10 confirms this prediction. One can also see that the proposed schemes with both analytical and numerical Green’s kernels perform better than the standard method in section 4.2.4 as expected.

6 Conclusion

In this paper, we considered the problem of approximating solutions to elliptic PDEs in domains comprised of deterministic and random boundaries. We used the numerical scheme based on boundary integral representation of solutions to such PDEs and proposed to use Green’s functions as the kernels of the corresponding potentials. The proposed numerical scheme can be applied to problems in arbitrary domains with random apertures and does not require the explicit knowledge of analytical Green’s functions which can be considered as a major contribution of this work. We showed that when a large number of repetitive solutions is required, as is the case of Monte Carlo simulations, the proposed scheme can lead to significant reduction of computational complexity compared to standard BIE techniques.

Data availability

The authors declare that all data supporting the findings of this study are available within the article.

Acknowledgements

This material is based upon work supported in part by: the U.S. Department of Energy, Office of Science, Early Career Research Program under award number ERKJ314; U.S. Department of Energy, Office of Advanced Scientific Computing Research under award numbers ERKJ331 and ERKJ345; the National Science Foundation, Division of Mathematical Sciences, Computational Mathematics program under contract number DMS1620280; and by the Laboratory Directed Research and Development program at the Oak Ridge National Laboratory, which is operated by UT-Battelle, LLC., for the U.S. Department of Energy under contract DE-AC05-00OR22725.