Abstract
To model subsurface flow in uncertain heterogeneous or fractured media an elliptic equation with a discontinuous stochastic diffusion coefficient—also called random field—may be used. In case of a one-dimensional parameter space, Lévy processes allow for jumps and display great flexibility in the distributions used. However, in various situations (e.g. microstructure modeling), a one-dimensional parameter space is not sufficient. Classical extensions of Lévy processes on two parameter dimensions suffer from the fact that they do not allow for spatial discontinuities [see for example Barth and Stein (Stoch Part Differ Equ Anal Comput 6(2):286–334, 2018)]. In this paper a new subordination approach is employed [see also Barth and Merkle (Subordinated gaussian random fields. ArXiv e-prints, arXiv:2012.06353 [math.PR], 2020)] to generate Lévy-type discontinuous random fields on a two-dimensional spatial parameter domain. Existence and uniqueness of a (pathwise) solution to a general elliptic partial differential equation is proved and an approximation theory for the diffusion coefficient and the corresponding solution provided. Further, numerical examples using a Monte Carlo approach on a Finite Element discretization validate our theoretical results.
Similar content being viewed by others
1 Introduction
Over the last decade partial differential equations with stochastic operators\ data\ domain became a widely studied object. This branch of research is oftentimes called uncertainty quantification. Especially for problems where data is sparse or measurement errors are unavoidable, like subsurface flow problems, the theory provides an approach to quantify this uncertainty. There are two main approaches to discretize the uncertain problem: intrusive and non-intrusive methods. The former require the solution of a high dimensional partial differential equation, where the dimensionality depends on the smoothness of the random field or process (see for example [7, 22, 30] and the references therein). The latter consist of (essentially) sampling methods and require multiple solutions of a low dimensional problem (see, among others [1, 9, 11, 12, 29, 34]). Up to date mainly Gaussian random fields were used to model the diffusivity in an elliptic equation (as a model for a subsurface flow problem). Gaussian random fields have the advantage that they may be used in both approaches and that they are stochastically very well understood objects. A great disadvantage is however, that the distributions underlying the field are Gaussian and therefore the model lacks flexibility, in the sense that the field cannot have pointwise marginal distributions with heavy-tails. Further, Gaussian random fields that are commonly used in applications, like Matérn fields, have \({\mathbb {P}}\)-almost surely spatial continuous paths. In view of these limitations, some extensions of Gaussian models for the diffusivity in elliptic PDEs have been explored in the literature: The paper [29] presents a multilevel Monte Carlo method for elliptic equations with jump-diffusion coefficients for two-phase random media, where the diffusion coefficient assumes two different values on random domains. Discontinuous models for the diffusion coefficient in elliptic PDEs have also been used in the context of Bayesian inversion, where geometric priors or level-set priors are used as a model for the diffusion, see e.g. [27, 28]. In the recent article [20], the diffusion coefficient of an elliptic PDE is modeled as a transformation of smoothed Lévy noise fields. The authors investigate the existence, integrability and approximation of the corresponding stochastic PDE solution. The presented model yields a high distributional flexibility of the diffusion coefficient, however, the considered fields are continuous in space.
In this paper we propose a two-dimensional subordinated Gaussian random field as stochastic diffusion coefficient in an elliptic equation. The subordinated Gaussian random field is a type of a (discontinuous) Lévy field. Different subordinators display unique patterns in the discontinuities and have varied marginal distributions (see [8]). Naturally the spatial regularity of a subordinated Gaussian random field depends on the subordinator. We prove existence and uniqueness of a solution to the elliptic equation in a pathwise sense and provide different discretization schemes.
We structured the rest of the paper as follows: in Sect. 2 we introduce a general pathwise existence and uniqueness result for a stochastic elliptic equation under mild assumptions on the coefficient. These assumptions accommodate the subordinated Gaussian random fields we introduce in Sect. 3. In Sect. 4 we approximate the specific diffusion coefficient which is used in this paper and show convergence of the elliptic equation with the approximated coefficient to the unapproximated solution in Sect. 5. Section 6 provides spatial approximation methods and in Sect. 7 numerical examples are presented.
2 The stochastic elliptic problem
In this section we introduce the framework of the general stochastic elliptic boundary value problem which allows for discontinuous diffusion coefficients. For the general setting and pathwise existence theory we follow Barth and Stein [11]. In the following, let \((\varOmega ,{\mathcal {F}},{\mathbb {P}})\) be a complete probability space. Let \((B,\Vert \cdot \Vert _B)\) be a Banach space and Z be a B valued random variable, i.e. a measurable function \(Z:\varOmega \rightarrow B\). The space \(L^p(\varOmega ;B)\) contains all strongly measurable B-valued random variables with \(\Vert Z\Vert _{L^p(\varOmega ;B)}<+\infty \), for \(p\in [1,+\infty ]\), where the norm is defined by
2.1 Problem formulation
Let \({\mathcal {D}}\subset {\mathbb {R}}^d\) for \(d\in {\mathbb {N}}\) be a bounded, connected Lipschitz domain. We consider the equation
where \(a:\varOmega \times {\mathcal {D}}\rightarrow {\mathbb {R}}_+\) is a stochastic (jump-diffusion) coefficient and \(f:\varOmega \times {\mathcal {D}}\rightarrow {\mathbb {R}}\) is a (measurable) random source function. Further, we impose the following boundary conditions
where we assume to have a decomposition \(\partial {\mathcal {D}}=\varGamma _1\overset{\cdot }{\cup }\varGamma _2\) with two \((d-1)\)-dimensional manifolds \(\varGamma _1,~\varGamma _2\) such that the exterior normal derivative \(\overrightarrow{n}\cdot \nabla u\) on \(\varGamma _2\) is well-defined for every \(u\in C^1(\overline{{\mathcal {D}}})\). Here, \(\overrightarrow{n}\) is the outward unit normal vector to \(\varGamma _2\) and \(g:\varOmega \times \varGamma _2\rightarrow {\mathbb {R}}\) a measurable function. Note that we just reduce the theoretical analysis to the case of homogeneous Dirichlet boundary conditions to simplify notation. The extension to the inhomogeneous case is straightforward.
We now state assumptions under which the elliptic boundary value problem has a unique solution.
Assumption 2.1
Let \(H:=L^2({\mathcal {D}})\). We assume that
-
i.
for any fixed \({\underline{x}}\in {\mathcal {D}}\) the mapping \(\omega \mapsto a(\omega ,{\underline{x}})\) is measurable, i.e. \(a(\cdot ,{\underline{x}})\) is a (real-valued) random variable,
-
ii.
for any fixed \(\omega \in \varOmega \) the mapping \(a(\omega ,\cdot )\) is \({\mathcal {B}}({\mathcal {D}})-{\mathcal {B}}({\mathbb {R}}_+)\)-measurable and it holds \(a_{-}(\omega ):=\underset{{\underline{x}}\in {\mathcal {D}}}{{\text {ess}}\,\inf }\,a(\omega ,{\underline{x}})>0\) and \(a_+(\omega ):=\underset{{\underline{x}}\in {\mathcal {D}}}{{\text {ess}}\,\sup }\,a(\omega ,{\underline{x}})<+\infty \),
-
iii.
\(\frac{1}{a_{-}}\in L^p(\varOmega ;{\mathbb {R}})\), \(f\in L^q(\varOmega ;H)\) and \(g\in L^q(\varOmega ;L^2(\varGamma _2))\) for some \(p,q\in [1,+\infty ]\) such that \(r:=(\frac{1}{p} + \frac{1}{q})^{-1}\ge 1\).
Remark 2.2
Note that Assumption 2.1 implies that the real-valued mappings \(a_-,a_+:\varOmega \rightarrow {\mathbb {R}}\) are measurable. This can be seen as follows: For fixed \(p\ge 1\) consider the mapping
which is well-defined by Assumption 2.1. It follows from the definition of the Lebesgue integral and Assumption 2.1(i) that the mapping \(\omega \mapsto I_p(\omega )\) is \({\mathcal {F}}-{\mathcal {B}}({\mathbb {R}})\) measurable. For a fixed \(\omega \in \varOmega \), by the embedding theorem for \(L^p\) spaces (see [2, Theorem 2.14]), we get
Since this holds for all \(\omega \in \varOmega \) we obtain by Aliprantis and Border [4, Lemma 4.29] that the mapping
is \({\mathcal {F}}-{\mathcal {B}}({\mathbb {R}})\)-measurable. The measurability of \(\omega \mapsto a_-(\omega )\) follows analogously. Note that we do not treat the random coefficient \(a:\varOmega \times {\mathcal {D}}\rightarrow {\mathbb {R}}\) as a \(L^\infty ({\mathcal {D}})\)-valued random variable, since \(L^\infty ({\mathcal {D}})\) is not separable and therefore the strong measurability of the mapping \(a:\varOmega \rightarrow L^\infty ({\mathcal {D}})\) is only guaranteed in a very restrictive setting. Nevertheless, the measurability of the functions \(a_+\) and \(a_-\) allows taking expectations of these real-valued random variables. In order to avoid confusion about that, we use the notation \({\mathbb {E}}(\mathop {\hbox {ess sup}}\nolimits _{{\underline{x}}\in {\mathcal {D}}} |\cdot |^s)^{1/s}\).
2.2 Weak solution
We denote by \(H^1({\mathcal {D}})\) the Sobolev space on \({\mathcal {D}}\) with the norm \(\Vert v\Vert _{H^1({\mathcal {D}})}:=(\Vert v\Vert _{L^2({\mathcal {D}})}^2 + \Vert \nabla v\Vert _{L^2({\mathcal {D}})}^2)^{1/2}\), for \(v\in H^1({\mathcal {D}})\) (see for example [21, Section 5.2]). Further, we denote by T the trace operator \(T:H^1({\mathcal {D}})\rightarrow H^{\frac{1}{2}}(\partial {\mathcal {D}})\) where \(Tv=v|_{\partial {\mathcal {D}}}\) for \(v\in C^\infty (\overline{{\mathcal {D}}})\) (see [18]). We define the subspace \(V\subset H^1({\mathcal {D}})\) as
with the standard Sobolev norm, i.e. \(\Vert \cdot \Vert _V:=\Vert \cdot \Vert _{H^1({\mathcal {D}})}\). We identify \(H=L^2({\mathcal {D}})\) with its dual space \(H'\) and work on the Gelfand triplet \(V\subset H\simeq H'\subset V'\). Hence, Assumption 2.1 guarantees that \(f(\omega ,\cdot )\in V'\) and \(g(\omega ,\cdot )\in H ^ {-\frac{1}{2}}(\varGamma _2)\) for \({\mathbb {P}}\)-almost every \(\omega \in \varOmega \). We multiply Eq. (2.1) by a test function \(v\in V\), integrate by parts and use the boundary conditions (2.2) and (2.3) to obtain
This leads to the following pathwise weak formulation of the problem: For any \(\omega \in \varOmega \), given \(f(\omega ,\cdot )\in V'\) and \(g(\omega ,\cdot )\in H^{-\frac{1}{2}}(\varGamma _2)\), find \(u(\omega ,\cdot )\in V\) such that
for all \(v\in V\). The function \(u(\omega ,\cdot )\) is then called pathwise weak solution to problem (2.1)–(2.3). Here, the bilinear form \(B_{a(\omega )}\) and the operator \(F_\omega \) are given by
and
for fixed \(\omega \in \varOmega \), where the integrals in \(F_\omega \) are understood as the duality pairings:
for \(v\in V\).
Theorem 2.3
(See [11, Theorem 2.5]) Under Assumption 2.1, there exists a unique pathwise weak solution \(u(\omega ,\cdot )\in V\) to problem (2.4) for very \(\omega \in \varOmega \). Furthermore, \(u\in L^r(\varOmega ;V)\) and
where \(C(a_-,{\mathcal {D}},p)>0\) is a constant depending on \(a_-\), p and the volume of \({\mathcal {D}}\).
In addition to the existence of the solution, the following remark gives a rigorous justification for the measurability of the solution mapping
which maps any \(\omega \in \varOmega \) on the corresponding pathwise weak PDE solution.
Remark 2.4
Let \((v_n,~n\in {\mathbb {N}})\subset V\) be an orthonormal basis of the separable Hilbert space V. For every \(n\in {\mathbb {N}}\) we define the mapping
It is easy to see that this mapping is Carathéodory for any \(n\in {\mathbb {N}}\), i.e. \(J_n(\cdot ,v)\) is \({\mathcal {F}}\)-\({\mathcal {B}}({\mathbb {R}})\)-measurable for any fixed \(v\in V\) and \(J_n(\omega ,\cdot )\) is continuous on V for any fixed \(\omega \in \varOmega \). We define the correspondences
for every \(n\in {\mathbb {N}}\). It follows from Aliprantis and Border [4, Corollary 18.8] that this correspondence has a measurable graph, i.e.
Further, by Assumption 2.1 and the Lax–Milgram lemma (see for example [26, Lemma 6.97] and [11, Theorem 2.5]) we know that for every fixed \(\omega \in \varOmega \), there exists a unique solution \(u(\omega ,\cdot )\in V\) satisfying (2.4) for every \(v\in V\). Therefore, we obtain for the graph of the solution mapping:
This implies for an arbitrary measurable set \({\tilde{V}}\in {\mathcal {B}}(V)\)
and therefore
by the projection theorem (see [4, Theorem 18.25]), which gives the measurability of the solution mapping.
3 Subordinated Gaussian random fields
A random field \(W:\varOmega \times {\mathcal {D}} \rightarrow {\mathbb {R}}\) is called an \({\mathbb {R}}-valued\) Gaussian random field (GRF) if for any tuple \(({\underline{x}}_1,\dots ,{\underline{x}}_n)\subset {\mathcal {D}}\) and any number \(n\in {\mathbb {N}}\) the \({\mathbb {R}}^n\)-valued random variable
is multivariate normally distributed (see [3, Section 1.2]). Here \({\underline{x}}^T\) denotes the transpose of the vector \({\underline{x}}\). We denote by
the associated mean and covariance function. The covariance operator \(Q:L^2({\mathcal {D}})\rightarrow L^2({\mathcal {D}})\) of W is defined by
Further, if \({\mathcal {D}}\subset {\mathbb {R}}^d\) is compact and W is centered, i.e. \(m\equiv 0\), there exists a decreasing sequence \((\lambda _i,~i\in {\mathbb {N}})\) of real eigenvalues of Q with corresponding eigenfunctions \((e_i,~i\in {\mathbb {N}})\subset L^2({\mathcal {D}})\) which form an orthonormal basis of \(L^2({\mathcal {D}})\) (see [3, Section 3.2] and [35, Theorem VI.3.2 and Chapter II.3]).
Example 3.1
One important class of continuous GRFs is given by the Matérn family: for a given smoothness parameter \(\nu > 1/2\), correlation parameter \(r>0\) and variance \(\sigma ^2>0\), the Matérn-\(\nu \) covariance function is given by \(q_M({\underline{x}},{\underline{y}})=\rho _M(\Vert {\underline{x}}-{\underline{y}}\Vert _2)\), for \(({\underline{x}},{\underline{y}})\in {\mathbb {R}}_+^d\times {\mathbb {R}}_+^d\), with
where \(\varGamma (\cdot )\) is the Gamma function and \(K_\nu (\cdot )\) is the modified Bessel function of the second kind (see [23, Section 2.2 and Proposition 1]). Here, \(\Vert {\underline{x}}\Vert _2:=(\sum _{i=1}^d {\underline{x}}_i^2)^\frac{1}{2}\) denotes the Euclidean norm of the vector \({\underline{x}}\in {\mathbb {R}}^d\). A Matérn-\(\nu \) GRF is a centered GRF with covariance function \(q_M\).
3.1 Construction of subordinated GRFs
Let \(D>0\) be fixed and \({\mathcal {I}}=[0,D]\) or \({\mathcal {I}}={\mathbb {R}}_+\) with \({\mathbb {R}}_+:=[0,+\infty )\). A real-valued stochastic process \(l=(l(t),~t\in {\mathcal {I}})\) is said to be a Lévy process on \({\mathcal {I}}\) if \(l(0)=0\) \({\mathbb {P}}\)-a.s., l has independent and stationary increments and l is stochastically continuous (see [5, Section 1.3]). One of the most important properties of Lévy processes is the so called Lévy–Khinchin formula.
Theorem 3.2
(Lévy–Khinchin formula, see [5, Th. 1.3.3]) Let l be a real-valued Lévy process on the interval \({\mathcal {I}}\subseteq {\mathbb {R}}_+\). There exist constants, \(\gamma _l \in {\mathbb {R}}\), \(\sigma _l^2\in {\mathbb {R}}_+\) and a measure \(\nu \) on \(({\mathbb {R}},{\mathcal {B}}({\mathbb {R}}))\) such that the characteristic function \(\phi _{l(t)}\), for \(t\in {\mathcal {I}}\), admits the representation
Motivated by Theorem 3.2 we denote by \((\gamma _l,\sigma _l^2,\nu )\) the characteristic triplet of the Lévy process l. A (Lévy-)subordinator is a Lévy process which is non-decreasing \({\mathbb {P}}\)-a.s. By [5, Theorem 1.3.15] it follows that the Lévy triplet of a Lévy-subordinator always admits the form \((\gamma _l,0,\nu )\) with a measure \(\nu \) on \(({\mathbb {R}},{\mathcal {B}}({\mathbb {R}}))\) satisfying
Remark 3.3
Let \(B=(B(t),~t\ge 0)\) be a standard Brownian motion and \(l=(l(t),~t\ge 0)\) be a Lévy subordinator. The stochastic process defined by
is called subordinated Brownian motion and is again a Lévy process (see [5, Theorem 1.3.25]).
In [8] the authors propose a new approach to extend standard subordinated Lévy processes on a higher dimensional parameter space. Motivated by the rich class of subordinated Brownian motions, the authors construct discontinuous random fields by subordinating a GRF on a d-dimensional parameter domain by d one-dimensional Lévy subordinators. In case of a two-dimensional parameter space the construction is as follows: For a GRF \(W:\varOmega \times {\mathbb {R}}_+^2\rightarrow {\mathbb {R}}\) and two (Lévy-)subordinators \(l_1,~l_2\) on [0, D], with a finite \(D>0\), we define the real-valued random field
Figure 1 shows samples of a GRF with Martérn-1.5 covariance function and the corresponding subordinated field where we used Poisson and Gamma processes as subordinators.
Sample of Matérn-1.5-GRF (left), Poisson-subordinated GRF (middle) and Gamma-subordinated GRF (right)
This construction yields a rich class of discontinuous random fields which also admit a Lévy–Khinchin-type formula. Further, the newly constructed random fields are also interesting for practical reasons, since formulas for the covariance functions can be derived which is very useful for applications, e.g. in statistical fitting. For a theoretical investigation of the constructed random fields we refer to Barth and Merkle [8].
3.2 Subordinated GRFs as diffusion coefficients in elliptic problems
In the following, we define the specific diffusion coefficient that we consider in problem (2.1)–(2.3). In order to allow discontinuities, we incorporate a subordinated GRF in the coefficient additionally to a Gaussian component. The construction of the coefficient is done so that Theorem 2.3 is applicable and, at the same time, the coefficient is as versatile as possible.
Definition 3.4
We consider the domain \({\mathcal {D}}=(0,D)^2\) with \(D<+\infty \).Footnote 1 We define the jump-diffusion coefficient a in problem (2.1)–(2.3) with \(d=2\) as
where
-
\({\overline{a}}:{\mathcal {D}}\rightarrow (0,+\infty )\) is deterministic, continuous and there exist constants \({\overline{a}}_+,{\overline{a}}_->0\) with \({\overline{a}}_-\le {\overline{a}}(x,y)\le {\overline{a}}_+\) for \((x,y)\in {\mathcal {D}}\).
-
\(\Phi _1,~\Phi _2:{\mathbb {R}}\rightarrow [0,+\infty )\) are continuous.
-
\(W_1\) and \(W_2\) are zero-mean GRFs on \({\mathcal {D}}\) respectively on \([0,+\infty )^2\) with \({\mathbb {P}}-a.s.\) continuous paths.
-
\(l_1\) and \(l_2\) are Lévy subordinators on [0, D] with Lévy triplets \((\gamma _1,0,\nu _1)\) and \((\gamma _2,0,\nu _2)\) which are independent of the GRFs \(W_1\) and \(W_2\).
Remark 3.5
The first two assumptions ensure that the diffusion coefficient a is positive over the domain \({\mathcal {D}}\). To show the convergence of the approximated diffusion coefficient in Sect. 5.1 we have to impose independence of the GRFs \(W_1\) and \(W_2\) (see Assumption 5.2 and the proof of Theorem 5.5). This assumption is in the sense natural as also one-dimensional Lévy processes admit an additive decomposition into a continuous part and a pure-jump part which are stochastically independent (Lévy–Itô decomposition, see e.g. [5, Theorem 2.4.11]). For the same reason the assumption that the Lévy subordinators are independent of the GRFs is also natural (see for example [5, Section 1.3.2]).
In order to verify Assumption 2.1(i) and (ii) we need the following Lemma. For a proof we refer to Barth and Merkle [8, Lemma 3.2 and Lemma 3.3]).
Lemma 3.6
For fixed \((x,y)\in {\mathcal {D}}\) the mapping \(\omega \mapsto a(\omega ,x,y)\) is \({\mathcal {F}} -{\mathcal {B}}({\mathbb {R}}_+)\)-measurable. Further, for fixed \(\omega \in \varOmega \), the mapping \((x,y)\mapsto a(\omega ,x,y)\) is \({\mathcal {B}}({\mathcal {D}})-{\mathcal {B}}({\mathbb {R}}_+)\)-measurable.
Definition 3.4 guarantees the existence of a pathwise weak solution to problem (2.1), as we prove in the following theorem.
Theorem 3.7
Let a be as in Definition 3.4 and let \(f\in L^q(\varOmega ;H),~g\in L^q(\varOmega ;L^2(\varGamma _2))\) for some \(q\in [1,+\infty )\). Then there exists a unique pathwise weak solution \(u(\omega ,\cdot )\in V\) to problem (2.1) for \({\mathbb {P}}\)-almost every \(\omega \in \varOmega \). Furthermore, \(u\in L^r(\varOmega ;V)\) for all \(r\in [1,q)\) and
where \(C({\overline{a}}_-,{\mathcal {D}})>0\) is a constant depending only on the indicated parameter and the volume of \({\mathcal {D}}\).
Proof
In order to apply Theorem 2.3 we have to verify Assumption 2.1. We have \(a_-(\omega )=\inf \limits _{x\in {\mathcal {D}}}a(\omega ,x)\ge {\overline{a}}_- \) for every fixed \(\omega \in \varOmega \) by Definition 3.4. Further, \(W_2(\omega )\) is continuous on \(K(\omega ):=[0,l_1(\omega ,D)]\times [0,l_2(\omega ,D)]\) and therefore
For \(1\le r<q\) define \(p:=(\frac{1}{r}-\frac{1}{q})^{-1}>0\). We observe that
\({\mathbb {P}}\)-a.s. and hence \(1/a_-\in L^p(\varOmega ;{\mathbb {R}})\). Therefore, Assumption 2.1 holds with \(r=(1/p + 1/q)^{-1}\) and the assertion follows by Theorem 2.3. \(\square \)
4 Approximation of the diffusion coefficient
To simulate the solution to the elliptic equation we need to define a tractable approximation of the diffusion coefficient. Here, we face a new challenge regarding the GRF \(W_2\) which is subordinated by the Lévy processes \(l_1\) and \(l_2\): due to the fact that the Lévy subordinators in general can attain any value in \([0,+\infty )\) we have to consider (and approximate) the GRF \(W_2\) on the unbounded domain \([0,+\infty )\). In most cases where elliptic PDEs of the form (2.1) have been considered with a random coefficient, the problem is stated on a bounded domain, see e. g. [11, 13, 14, 23, 25]. Many regularity results for GRFs formulated for a bounded parameter space cannot easily be transferred to an unbounded parameter space (see also [3, Chapter 1], especially the discussion on p. 13). Even the Karhunen–Loève expansion of a GRF requires compactness of the domain (see e.g. [3, Section 3.2]). We avoid an unbounded parameter domain by bounding the subordinators from above in Sect. 4.1. Furthermore, to show convergence of the solution in Sect. 5 we need to bound the diffusion coefficient from above by a deterministic upper bound A (see Theorem 5.11 and Remark 5.9). Note that the weaker assumption of \(L^p\) integrability of the pathwise upper bound \(a_+\) instead of a deterministic upper bound A of the diffusion coefficient is not enough in this setting. The reason being that the absence of a deterministic upper bound A would lead to an additional \(\omega \)-dependence in the regularity constants in Assumption 5.7 which would be difficult to handle in the convergence theorem in Sect. 5.2 (see also Remark 5.9). Therefore, we cut off the diffusion coefficient at a deterministic level A (see Subsection 4.2). Subsequently we show that this induces an error in the solution approximation which can be controlled and which vanishes for growing A (see Sect. 5.1 and Theorem 5.1). Summarized, we derive an approximation in three steps: First, we bound the subordinators, and, second, we cut off the diffusion coefficient. Finally, we consider approximations of the GRFs and the subordinators and prove the convergence of this approximation of the diffusion coefficient under suitable assumptions.
4.1 First approximation: bounding the Lévy subordinators
For a fixed \(K\in (0,+\infty )\), we define the cut-function \(\chi _K:[0,+\infty )\rightarrow [0,K]\) as \(\chi _K(z):=\min (z,K)\) for \(z\in [0,+\infty )\). Instead of problem (2.1) we consider the following modified problem
and impose the boundary conditions
Here, the diffusion coefficient is defined by
For functions \(f\in L^q(\varOmega ;H)\) and \(g\in L^q(\varOmega ;L^2(\varGamma _2))\) with \(q\in [1,+\infty )\), there exists a weak solution \(u_K\in L^r(\varOmega ;V)\) to problem (4.1)–(4.3) for \(r\in [1,q)\) (see Theorem 3.7).Footnote 2
Remark 4.1
We note that the influence of this problem modification can be controlled: one may choose \(K>0\) such that
for any \(\varepsilon >0\). In other words, pathwise the modified problem coincides with the original one up to a set of samples, whose probability can be made arbitrarily small.
4.2 Second modification: diffusion cut-off
We consider again the cut function \(\chi _{A}(z) :=\min (z,A)\) for \(z\in [0,+\infty )\) with a fixed positive number \(A>0\) and consider the following problem
where we impose the boundary conditions
The diffusion coefficient \(a_{K,A}\) is defined by
Again, Theorem 3.7 applies in this case and yields the existence of a pathwise weak solution \(u_{K,A}\in L^r(\varOmega ;V)\) for \(r\in [1,q)\) if \(f\in L^q(\varOmega ;H)\) and \(g\in L^q(\varOmega ;L^2(\varGamma _2))\). The error of the modification vanishes for growing A as shown in the end of Sect. 5.1.
4.3 Approximation of GRF and subordinators
In the following, we show how to approximate the modified diffusion coefficient \(a_{K,A}\) using approximations \(W_1^{\varepsilon _W}\approx W_1\), \(W_2^{\varepsilon _W}\approx W_2\) of the GRFs and \(l_1^{\varepsilon _l}\approx l_1\), \(l_2^{\varepsilon _l}\approx l_2\) of the Lévy subordinators. We aim to approximate the GRFs \(W_1\) and \(W_2\) in the diffusion coefficient by sampling on a discrete grid. These approximations may be extended to the whole domain using linear interpolation (see [24, 25]). Certain regularity assumptions on the GRFs then allow for a quantification of the corresponding approximation error (see Lemma 4.4). In [10], the authors considered approximations of general Lévy processes and quantified the approximation error in an \(L^s(\varOmega )\)-sense uniformly in \(x\in [0,D]\). Such an approximation may be obtained by piecewise constant approximations constructed from samples of the process \((l_j(x_i),~i=1,\dots ,N_l)\) on a discrete grid \((x_i,~i=1,\dots ,N_l)\subset [0,D]\) with \(N_l\in {\mathbb {N}}\) grid points, for \(j=1,2\) (see Sect. 7, [10, 33]). Motivated by this, we formulate the following working assumptions on the covariance operators of the Gaussian fields, the subordinators and the data of the elliptic problem.
Assumption 4.2
Let \(W_1\) be a zero-mean GRF on \([0,D]^2\) and \(W_2\) be a zero-mean GRF on \([0,K]^2\). We denote by \(q_1:[0,D]^2\times [0,D]^2\rightarrow {\mathbb {R}}\) and \(q_2:[0,K]^2\times [0,K]^2\rightarrow {\mathbb {R}}\) the covariance functions of these random fields and by \(Q_1,Q_2\) the associated covariance operators defined by
for \(\phi \in L^2([0,z_j]^2)\) with \(z=(D,K)\) and \(j=1,2\). We denote by \((\lambda _i^{(1)},e_i^{(1)},~i\in {\mathbb {N}})\) resp. \((\lambda _i^{(2)},e_i^{(2)},~i\in {\mathbb {N}})\) the eigenpairs associated to the covariance operators \(Q_1\) and \(Q_2\). In particular, \((e_i^{(1)},~i\in {\mathbb {N}})\) resp. \((e_i^{(2)},~i\in {\mathbb {N}})\) are ONBs of \(L^2([0,D]^2)\) resp. \(L^2([0,K]^2)\).
-
i.
We assume that the eigenfunctions are continuously differentiable and there exist positive constants \(\alpha , ~\beta , ~C_e, ~C_\lambda >0\) such that for any \(i\in {\mathbb {N}}\) it holds
$$\begin{aligned}&\Vert e_i^{(1)}\Vert _{L^\infty ([0,D]^2)},~\Vert e_i^{(2)}\Vert _{L^\infty ([0,K]^2)}\le C_e,\\&\Vert \nabla e_i^{(1)}\Vert _{L^\infty ([0,D]^2)},~\Vert \nabla e_i^{(2)}\Vert _{L^\infty ([0,K]^2)}\le C_e i^\alpha ,~ \\&\sum _{i=1}^ \infty (\lambda _i^{(1)}+ \lambda _i^{(2)})i^\beta \le C_\lambda < + \infty . \end{aligned}$$ -
ii.
There exist constants \(\phi ,~\psi , C_{lip}>0\) such that the continuous functions \(\Phi _1,~\Phi _2:{\mathbb {R}}\rightarrow [0,+\infty )\) from Definition 3.4 satisfy
$$\begin{aligned} |\Phi _1'(x)|\le \phi \, \exp (\psi |x|),~ |\Phi _2(x)-\Phi _2(y)|\le C_{lip}\,|x-y| \text { for } x,y\in {\mathbb {R}}. \end{aligned}$$In particular, \(\Phi _1\in C^1({\mathbb {R}})\).
-
iii.
\(f\in L^q(\varOmega ;H)\) and \(g\in L^q(\varOmega ;L^2(\varGamma _2))\) for some \(q\in (1,+\infty ).\)
-
iv.
\({\overline{a}}:{\mathcal {D}}\rightarrow (0,+\infty )\) is deterministic, continuous and there exist constants \({\overline{a}}_+,{\overline{a}}_->0\) with \({\overline{a}}_-\le {\overline{a}}(x,y)\le {\overline{a}}_+\) for \((x,y)\in {\mathcal {D}}\).
-
v.
\(l_1\) and \(l_2\) are Lévy subordinators on [0, D] with Lévy triplets \((\gamma _1,0,\nu _1)\) and \((\gamma _2,0,\nu _2)\) which are independent of the GRFs \(W_1\) and \(W_2\). Further, we assume that we have approximations \(l_1^{(\varepsilon _l)},~l_2^{(\varepsilon _l)}\) of these processes and there exist constants \(C_l>0\) and \(\eta >1\) such that for every \(s\in [1,\eta )\) it holds
$$\begin{aligned} {\mathbb {E}}(|l_j(x)-l_j^{(\varepsilon _l)}(x)|^s)\le C_l\varepsilon _l, \end{aligned}$$for \(\varepsilon _l >0\), \(x\in [0,D]\) and \(j=1,2\).
Remark 4.3
Note that the first assumption on the eigenpairs of the GRFs is natural (see [11, 23]). For example, the case that \(Q_1,~Q_2\) are Matérn covariance operators are included. Assumption 4.2(ii) is necessary to be able to quantify the error of the approximation of the diffusion coefficient. Assumption 4.2(iii) is necessary to ensure the existence of a solution and has already been formulated in Assumption 2.1. The last assumption ensures that we can approximate the Lévy subordinators in an \(L^s\)-sense. This can always be achieved under appropriate assumptions on the tails of the distribution of the subordinators, see [10, Assumption 3.6, Assumption 3.7 and Theorem 3.21].
For any numerical simulation we have to approximate the GRF as well as the subordinating Lévy processes, which results in an additional approximation of the coefficient \(a_{K,A}\) given in Eq. (4.8). In the following we want to quantify the error induced by this approximation.
It follows by an application of the Kolmogorov–Chentsov theorem ([16, Theorem 3.5]) that \(W_1\) and \(W_2\) can be assumed to have Hölder-continuous paths with Hölder exponent \(b \in (0,(2\gamma k-2)/(2k))\) for \(0<\gamma \le \min (1,\beta /(2\alpha ))\) and every \(k\in {\mathbb {N}}\) (see the proof of [11, Lemma 3.5]). Further, it follows by an application of the Sobolev embedding theorem that \(W_1\in L^n(\varOmega ;C^{0,\gamma }([0,D]^2))\) and \(W_2\in L^n(\varOmega ;C^{0,\gamma }([0,K]^2))\), i.e.
for every \(n\in [1,+\infty )\) and \(\gamma < \min (1,\beta /(2\,\alpha ))\) (see [13, Proposition 3.1]).
Next, we prove a bound on the error of the approximated diffusion coefficient, where the GRFs are approximated by a discrete evaluation and (bi-)linear interpolation between these points (see [24, 25]).
Lemma 4.4
We consider the discrete grids \(G_1^{(\varepsilon _W)}=\{(x_i,x_j)|~i,j=0,\dots ,M_{\varepsilon _W}^{(1)}\}\) on \([0,D]^2\) and \(G_2^{(\varepsilon _W)}=\{(y_i,y_j)|~i,j=0,\dots ,M_{\varepsilon _W}^{(2)}\}\) on \([0,K]^2\) where \((x_i,~i=0,\dots ,M_{\varepsilon _W}^{(1)})\) is an equidistant grid on [0, D] with maximum step size \(\varepsilon _W\) and \((y_i,~i=0,\dots ,M_{\varepsilon _W}^{(2)})\) is an equidistant grid on [0, K] with maximum step size \(\varepsilon _W\). Further, let \(W_1^{(\varepsilon _W)}\) and \(W_2^{(\varepsilon _W)}\) be approximations of the GRFs \(W_1,~W_2\) on the discrete grids \(G_1^{(\varepsilon _W)}\) resp. \(G_2^{(\varepsilon _W)}\) which are constructed by point evaluation of the random fields \(W_1\) and \(W_2\) on the grids and linear interpolation between the grid points. Under Assumption 4.2(i) it holds for \(n\in [1,+\infty )\):
for \(\gamma < \min (1,\beta /(2\alpha ))\) where \(\beta \) and \(\alpha \) are the parameters from Assumption 4.2.
Proof
Note that for any fixed \(\omega \in \varOmega \) and \({\mathcal {D}}_{ij}:=[x_i,x_{i+1}]\times [y_j,y_{j+1}]\) with \(i,j\in \{1,\dots ,M_{\varepsilon _W}^{(1)}\}\), it holds
This holds since \(W_1^{(\varepsilon _W)}\) is constructed by (bi-)linear interpolation of the GRF \(W_1\) and the piecewise linear interpolants attain their maximum and minimum at the corners (the Hessian evaluated at the (unique) stationary point of the bilinear basis functions is always indefinite). Therefore, for a fixed \((x,y)\in [x_i,x_{i+1}]\times [y_j,y_{j+1}]\) it follows from the intermediate value theorem that \(W_1^{(\varepsilon _W)}(x,y)=W_1(x',y')\) for appropriate \((x',y')\in [x_i,x_{i+1}]\times [y_j,y_{j+1}]\). Using this observation we estimate
where we used Eq. (4.9) in the last step. Equivalently the error bound for \(W_2\) follows. \(\square \)
Remark 4.5
Note that Lemma 4.4 immediately implies for \(m\in [1,+\infty )\)
4.4 Convergence to the modified diffusion coefficient
Given some approximations \(W_1^{(\varepsilon _W)}\approx W_1,~W_2^{(\varepsilon _W)}\approx W_2\) as in Lemma 4.4 and approximations \(l_1^{(\varepsilon _l)}\approx l_1,~l_2^{(\varepsilon _l)}\approx l_2\) as in Assumption 4.2v as well as some fixed constants \(K,A>0\), we approximate the diffusion coefficient \(a_{K,A}\) in (4.8) by \(a_{K,A}^{(\varepsilon _W,\varepsilon _l)}:\varOmega \times {\mathcal {D}}\rightarrow (0,+\infty )\) with
for \((x,y)\in {\mathcal {D}}\). To prove a convergence result for this approximated coefficient (Theorem 4.8) we need the following two technical lemmas. The second can be proved by the use of Sato [32, Proposition 1.16]. For a detailed proof we refer to Barth and Merkle [8].
Lemma 4.6
For \(n,m\in [1,+\infty )\) with \(n\ge m\) and \(\varphi \in L^n([0,D]^2\times \varOmega ;{\mathbb {R}},\lambda \otimes {\mathbb {P}}))\), where \(\lambda \) denotes the Lebesgue measure on \(({\mathbb {R}}^2,{\mathcal {B}}({\mathbb {R}}^2))\), it holds
Proof
The case \(n=m\) is trivial. For \(n>m\) we use Hölder’s inequality and obtain
\(\square \)
Lemma 4.7
Let \(W:\varOmega \times {\mathbb {R}}_+^d\rightarrow {\mathbb {R}}\) be a \({\mathbb {P}}\)-a.s. continuous random field and let \(Z:\varOmega \rightarrow {\mathbb {R}}_+^d\) be a \({\mathbb {R}}_+^d\)-valued random variable which is independent of the random field W. Further, let \(\varphi :{\mathbb {R}}\rightarrow {\mathbb {R}}\) be a deterministic, continuous function. It holds
where \(\zeta (z):={\mathbb {E}}(\varphi (W(z))\) for deterministic \(z\in {\mathbb {R}}_+^d\).
Theorem 4.8
Let \(W_1^{(\varepsilon _W)}\approx W_1,~W_2^{(\varepsilon _W)}\approx W_2\) be approximations of the GRFs on discrete grids as in Lemma 4.4. Further, let \(1\le t\le s<\eta \) and \(0<\gamma <min(1,\beta /(2\alpha ))\) such that \(s\gamma \ge 2\). Under Assumption 4.2 we get the following error bound for the approximation of the diffusion coefficient:
with a constant C which does not depend on the discretization parameters \(\varepsilon _W\) and \(\varepsilon _l\).
Proof
Since the cut function \(\chi _{A}\) is Lipschitz continuous with Lipschitz constant 1 we calculate
First, we consider \(I_1\) and use Assumption 4.2(ii) and the same calculation as in Remark 4.5 to get
The mean value theorem yields, for fixed \((x,y)\in [0,D]^2\) and an appropriately chosen value \(\xi \in (\min (W_1(x,y),W_1^{(\varepsilon _W)}(x,y)),\max (W_1(x,y),W_1^{(\varepsilon _W)}(x,y)))\),
for \({\mathbb {P}}\)-almost every \(\omega \in \varOmega \). As already mentioned in the proof of Lemma 4.4, for any \(\omega \in \varOmega \) and fixed \({\mathcal {D}}_{ij}:=[x_i,x_{i+1}]\times [y_j,y_{j+1}]\) with \(i,j\in \{1,\dots ,M_{\varepsilon _W}^{(1)}\}\), it holds
Therefore, we obtain the pathwise estimate
Finally, we obtain for any \(n_1,n_2\in [1,+\infty )\) with \(1/n_1 + 1/n_2 = 1\) by Hölder’s inequality
where we used Lemma 4.4 and the fact that \(\Vert \exp (\psi |W_1|)\Vert _{L^{sn_1}(\varOmega ;L^\infty ([0,D]^2))}<\infty \) (see [3, Theorem 2.1.1] and the proof of Theorem 5.5 for more details).
For the second summand we calculate:
We use the same calculation as in Remark 4.5 and the Lipschitz continuity of \(\Phi _2\) to calculate for the summand \(I_4\)
where we used Lemma 4.4. It remains to bound the summand \(I_3\). We estimate using Lemma 4.6
We know by Lemma 4.7 that it holds for \((x,y)\in [0,D]^2\)
where
For \((x,y)=(v,w)\) it holds \(\varphi (x,y,v,w)=0\) and for \((x,y)\ne (v,w)\in [0,K]^2\) we get
by Eq. (4.9). Further, we know from Hölder’s inequality for \(\gamma s\ge 2\) that it holds
and therefore we calculate
where we used the Lipschitz continuity of \(\chi _K\) and Assumption 4.2(v) in the last step. Therefore, we finally obtain
which proves that
\(\square \)
5 Convergence analysis
In this section we derive an error bound for the approximation of the solution. We split the error in two components: the first component is associated with the cut-off of the diffusion coefficient we described in Sect. 4.2. The second error contributor corresponds to the approximation of the GRFs and the Lévy subordinators we considered in Sect. 4.4.
Let \(r\in [1,q)\) with q as in Assumption 4.2(iii) and denote by \(u_K\in L^r(\varOmega ;V)\) the weak solution to problem (4.1)–(4.4). Further, let \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\in L^r(\varOmega ;V)\) be the weak solution to the problem
with boundary conditions
Note that Theorem 3.7 also applies to the elliptic problem with coefficient \(a_{K,A}^{(\varepsilon _W,\varepsilon _l)}\). The aim of this section is to quantify the error of the approximation \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)} \approx u_K\).Footnote 3 By the triangle inequality we obtain
for any suitable norm \(\Vert \cdot \Vert \). Here, \(u_{K,A}\) is the solution to the truncated problem (4.5)–(4.8). We consider the two error contributions \(E_1\) and \(E_2\) separately in this section. The results presented in Sects. 5.1 and 5.2 prove the following theorem.
Theorem 5.1
Under Assumptions 5.2 and 5.7 it holds
for admissible values of \(r\ge 2\) (see Remark 5.6 and Theorem 5.11). As a consequence, we obtain
for any \(\delta >0\). Furthermore, it holds
Therefore, the solution \(u_{K,A}\) to the modified problem (4.5)–(4.7) coincides with the solution \(u_K\) to problem (4.1)–(4.3) up to a set of samples, whose probability can be made arbitrarily small for growing A. With Remark 4.1 we conclude that the solution \(u_{K,A}\) coincides with the solution u to problem (2.1)–(2.3) with diffusion coefficient a defined in Eq. (3.1) up to a set of samples, whose probability can be made arbitratily small for growing thresholds K and A.
5.1 Bound on \(E_1\)
Assumption 5.2
We assume that the GRFs \(W_1\) and \(W_2\) occurring in the diffusion coefficient (3.1) are stochastically independent.
The aim of this subsection is to show that the first error contributor \(E_1\) in Eq. (5.4) vanishes for increasing cut-off threshold A. The strategy consists of two separated steps: in the first step we show the stability of the solution, which means that the value \(E_1\) can be controlled by the quality of the approximation of the diffusion coefficient \(a_{K,A}\approx a_K\). In the second step, we show that the quality of the approximation of the diffusion coefficient can be controlled by the cut-off threshold A. The first step is given by Theorem 5.4. In order to prove it we need the following lemma.
Lemma 5.3
For fixed cut-off levels A and K we consider the solution \(u_K\in L^r(\varOmega ;V)\) and its approximation \(u_{K,A}\in L^r(\varOmega ;V)\) for \(r\in [1,q)\). It holds the pathwise estimate
for \({\mathbb {P}}\)-almost every \(\omega \in \varOmega \). Here, the constant \(C({\overline{a}}_-,{\mathcal {D}})\) only depends on the indicated parameters and we define \(a_{K,+}(\omega ):=\max \{1, \underset{(x,y)\in {\mathcal {D}}}{{\text {ess}}\,\sup }\,a_K(\omega ,x,y)\}<\infty \) for \(\omega \in \varOmega \).
Proof
For a fixed \(\omega \in \varOmega \) we consider the variational problem: find a unique \({\hat{w}}\in V\) such that
for all \(v\in V\). By the Lax–Milgram theorem there exists a unique solution \({\hat{w}}\in V\) with
(see Theorem 3.7 and [11, Theorem 2.5]). Therefore, we obtain by Hölder’s inequality
where we used Young’s inequality in the last step. Finally, we obtain
\(\square \)
Theorem 5.4
Let \(f\in L^q(\varOmega ;H)\) and \(g\in L^q(\varOmega ;L^2(\varGamma _2))\) for some \(q\in [1,+\infty )\). Further, for a given number \(t\in (1,+\infty )\) we define the dual number \(n:=\frac{t}{t-1}\). Then, for any for \(r\in [1,q/n)\), holds
Proof
By a direct calculation we obtain
Since \(u_K\) and \(u_{K,A}\) are weak solutions of problem (4.1)–(4.4) resp. (4.5)–(4.8) it holds
\({\mathbb {P}}\)-almost surely and therefore
We estimate using Hölder’s inequality
and therefore we obtain
Using Lemma 5.3 we obtain the pathwise estimate
Using again Hölder’s inequality we have
By assumption it holds \(nr<q\). Therefore, we can choose a real number \(\rho >1\) such that \(nr\rho <q\). We define the dual number \(\rho ':=\frac{\rho }{\rho -1}\in (1,+\infty )\) and use Hölder’s inequality to calculate
Obviously Theorem 3.7 applies also to problem (4.1)–(4.4). Therefore, since \(nr\rho <q\) by assumption we conclude
where we additionally used the fact that \({\mathbb {E}}(a_{K,+}^{nr\rho '})^{\frac{1}{nr\rho '}}<+\infty \) (see Theorem 5.5). \(\square \)
In other words, finding a bound for the error contribution \(E_1\) in Eq. (5.4) reduces to quantifying the quality of the approximation of the diffusion coefficient \(a_{K,A}\approx a_K\). For readability the proof of the following theorem can be found in Appendix A.
Theorem 5.5
For any \(n\in (1,+\infty )\) it holds
Further, for any \(\delta >0\) there exists a constant \(A=A(\delta ,n)>0\) such that
We close this subsection with the following remark on the convergence of the error contributor \(E_1\) in Eq. (5.4).
Remark 5.6
From Theorem 5.4 together with Theorem 5.5 we obtain
for every \(r<q\).
5.2 Bound on \(E_2\)
The aim is to bound the second term of Eq. (5.4) given by
in an appropriate norm. For technical reasons we have to impose an additional assumption on the solution of the truncated problem. The subsequent remarks discuss situations under which this assumption is fulfilled.
Assumption 5.7
We assume that there exist constants \(j_{reg}>0\) and \(k_{reg}\ge 2\) such that
Since \(u_{K,A}\in H^1({\mathcal {D}})\) we already know that \(\nabla u_{K,A}\in L^2({\mathcal {D}})\). Assumption 5.7 requires a slightly higher integrability over the spatial domain. Since this is an assumption on the regularity of the solution \(u_{K,A}\) we denote the above constant by \(C_{reg}\).
Remark 5.8
Note that Assumption 5.7 is fulfilled if there exists \(\theta \in (0,1)\) such that
with some constant C(f, a) with \({\mathbb {E}}(C(f,a)^{k_{reg}})<\infty \). This is true since for any \(\rho \ge 2\) and an arbitrary function \(\varphi \in H^{1+\theta }({\mathcal {D}})\) the inequality
holds for \(\theta :=1-\frac{2}{\rho }\) (see [17, Theorem 6.7]). Here, the constant \(C=C({\mathcal {D}}, \theta )\) depends only on the indicated parameters. Hence, the condition (5.6) implies Eq. (5.5) with \(j_{reg}=\frac{2\theta }{1-\theta }\).
Remark 5.9
Note that there are several results about higher integrability of the gradient of the solution to an elliptic PDE of the form (4.5)–(4.8). For instance [31] yields that the solution \(u_{K,A}\) has \(H^{1+\delta /(2\pi )}\) regularity with \(\delta =min(1,{\overline{a}}_-,A^{-1})\) under mixed boundary conditions and under the assumption that \(a_{K,A}\) is piecewise constant (see [31, Theorem 7.3]). This corresponds to the case where no Gaussian noise is considered (i.e. \(\Phi _1\equiv 0\)) and \({\overline{a}}\) is constant. Another important result is given in [19]. It follows by Dreyfuss [19, Theorem 1] that under the assumption that there exists \(q>2\) with \(f\in L^q({\mathcal {D}})\) \({\mathbb {P}}-a.s.\) there exists a constant \(C=C({\mathcal {D}},\Vert f\Vert _{L^q({\mathcal {D}})},{\overline{a}}_-,A)\) and a positive number \(\vartheta =\vartheta ({\mathcal {D}},\Vert f\Vert _{L^q({\mathcal {D}})},{\overline{a}}_-,A)>0\) only depending on the indicated parameters, such that:
In particular, if the right hand side f of the problem is deterministic, then \(\vartheta \) and the constant C in (5.7) are deterministic and one immediately obtains
for any \(k_{reg}\ge 1\) and a deterministic, positive constant \(\vartheta >0\).
Next, we show that for a given approximation of the diffusion coefficient, the resulting error contributor \(E_2\) is bounded by the approximation error of the diffusion coefficient. Similar to the corresponding assertion we gave in Sect. 5.1 we need the following Lemma for the proof of this error bound. For a proof we refer to Lemma 5.3.
Lemma 5.10
For fixed cut-off levels A, K and fixed approximation parameters \(\varepsilon _W,\varepsilon _l\) we consider the PDE solutions \(u_{A,K}\in L^r(\varOmega ;V)\) and \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\in L^r(\varOmega ;V)\) for \(r\in [1,q)\). It holds the pathwise estimate
for \({\mathbb {P}}\)-almost every \(\omega \in \varOmega \). Here, constant \(C({\overline{a}}_-,{\mathcal {D}})\) depends only on the indicated parameters and \(a_{K,+}(\omega ):=\max \{1, \underset{(x,y)\in {\mathcal {D}}}{{\text {ess}}\,\sup }\,a_K(\omega ,x,y)\}<\infty \) for \(\omega \in \varOmega \).
Theorem 5.11
Let \(r\ge 2\) and \(b,c\in [1,+\infty ]\) be given such that it holds
with a fixed real number \(\gamma \in (0,min(1,\beta /(2\alpha ))\). Here, the parameters \(\eta , \alpha \) and \(\beta \) are determined by the GRFs \(W_1\), \(W_2\) and the Lévy subordinators \(l_1\), \(l_2\) (see Assumption 4.2).
Let \(m,n\in [1,+\infty ]\) be real numbers such that
and let \(k_{reg}\ge 2\) and \(j_{reg}>0\) be the regularity specifiers given by Assumption 5.7. If it holds that
then the approximated solution \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\) converges to the solution \(u_{K,A}\) of the truncated problem for \(\varepsilon _W,\varepsilon _l\rightarrow 0\) and it holds
Proof
By a direct calculation we obtain the pathwise estimate
Since \(u_{K,A}\) (resp. \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\)) is the weak solution to problem (4.5)–(4.8) [resp. (5.1)–(5.3)] we have
\({\mathbb {P}}\)-a.s. and therefore
Using Hölder’s inequality we calculate
and therefore
Next, we apply Lemma 5.10 to obtain the following estimate.
Hence, it remains to bound the norm \(\Vert (a_{K,A}^{(\varepsilon _W,\varepsilon _l)}-a_{K,A})\nabla u_{K,A}\Vert _{L^2({\mathcal {D}})}\). By Hölder’s inequality we obtain
Applying Hölder’s inequality once more we estimate
By assumption it holds \(rm<k_{reg}\). Hence, we can choose a real number \(\rho >1\) such that \(rm\rho <k_{reg}\). We define the dual number \(\rho ':=\frac{\rho }{1-\rho }\in (1,+\infty )\) and use Hölder’s inequality to obtain
where we again used the fact that \({\mathbb {E}}(a_{K,+}^{nr\rho '})^{\frac{1}{nr\rho '}}<+\infty \) (see Theorem 5.5) together with Assumption 5.7. Finally we obtain the estimate
where we applied Theorem 4.8 in the last estimate. \(\square \)
We close this section with a remark on how to choose the parameters A, \(\varepsilon _W\) and \(\varepsilon _l\) to obtain an approximation error smaller than any given threshold \(\delta >0\).
Remark 5.12
For any given parameter K large enough (see Remark 4.1), we choose a positive number \(A>0\) such that the first error contributor satisfies \(E_1=\Vert u_K-u_{K,A}\Vert _{L^r(\varOmega ;V)}< \delta /2\) (see Theorems 5.4 and 5.5). Afterwards, under the assumptions of Theorem 5.11, we may choose the approximation parameters \(\varepsilon _W\) and \(\varepsilon _l\) small enough, such that the second error contributor satisfies \(E_2 = \Vert u_{K,A}-u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\Vert _{L^r(\varOmega ;V)}<\delta /2\). Hence, we get an overall error smaller than \(\delta \) [see Eq. (5.4)].
6 Pathwise sample-adapted finite element approximation
We want to approximate the solution u to the problem (2.1)–(2.3) with diffusion coefficient a given by Eq. (3.1) using a pathwise Finite Element (FE) approximation of the solution \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\) of problem (5.1)–(5.3) where the approximated diffusion coefficient \(a_{K,A}^{(\varepsilon _W,\varepsilon _l)}\) is given by (4.10). Therefore, for almost all \(\omega \in \varOmega \), we have to find a function \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}(\omega ,\cdot )\in V\) such that it holds
for every \(v\in V\). Here, \(K,A,\varepsilon _W,\varepsilon _l\) are fixed approximation parameters. In order to solve this variational problem numerically we consider a standard Galerkin approach with linear elements and assume \({\mathcal {V}}=(V_\ell ,~\ell \in {\mathbb {N}}_0)\) to be a sequence of finite-dimensional subspaces \(V_\ell \subset V\) with \(\dim (V_\ell )=d_\ell \) for \(\ell \in {\mathbb {N}}_0\). We denote by \((h_\ell ,~\ell \in {\mathbb {N}}_0)\) the corresponding sequence of refinement sizes which is assumed to decrease monotonically to zero for \(\ell \rightarrow \infty \). Let \(\ell \in {\mathbb {N}}_0\) be fixed and denote by \(\{v_1^{(\ell )},\dots ,v_{d_\ell }^{(\ell )}\}\) a basis of \(V_\ell \). The (pathwise) discrete version of (6.1) reads:
Find \(u_{K,A,\ell }^{(\varepsilon _W,\varepsilon _l)}(\omega ,\cdot )\in V_\ell \) such that
We expand the function \(u_{K,A,\ell }^{(\varepsilon _W,\varepsilon _l)}(\omega ,\cdot )\) with respect to the basis \(\{v_1^{(\ell )},\dots ,v_{d_\ell }^{(\ell )}\}\):
where the coefficient vector \(\mathbf{c }=(c_1,\dots ,c_{d_\ell })^T\in {\mathbb {R}}^{d_\ell }\) is determined by the linear equation system
with a stochastic stiffness matrix \(\mathbf{B }(\omega )_{i,j}=B_{a_{K,A}^{(\varepsilon _W,\varepsilon _l)}(\omega )}(v_i^{(\ell )},v_j^{(\ell )})\) and load vector \({\mathbf {F}}(\omega )_i=F_\omega (v_i^{(\ell )})\) for \(i,j=1,\dots ,d_\ell \).
Remark 6.1
Let \(({\mathcal {K}}_\ell ,~\ell \in {\mathbb {N}}_0)\) be a sequence of triangulations on \({\mathcal {D}}\) and denote by \(\theta _\ell >0\) the minimum interior angle of all triangles in \({\mathcal {K}}_\ell \). We assume \(\theta _\ell \ge \theta >0\) for a positive constant \(\theta \) and define the maximum diameter of the triangulation \({\mathcal {K}}_\ell \) by \(h_\ell :=\max \limits _{K\in {\mathcal {K}}_\ell }\, \text {diam}(K),\) for \(\ell \in {\mathbb {N}}_0\). Further, we define the finite dimensional subspaces by \(V_\ell :=\{v\in V~|~v|_K\in {\mathcal {P}}_1,K\in {\mathcal {K}}_\ell \},\) where \({\mathcal {P}}_1\) denotes the space of all polynomials up to degree one. The convergence of the FE method is determined by the regularity of the solution: If we assume that for \({\mathbb {P}}-\)almost all \(\omega \in \varOmega \) it holds \(u_{K,A}^{(\varepsilon _W,\varepsilon _l)}(\omega ,\cdot )\in H^{1+\kappa _a}({\mathcal {D}})\) for some positive number \(\kappa _a>0\), the pathwise discretization error is bounded by Céa’s lemma \({\mathbb {P}}\)-a.s. by
(see [11, Section 4] and [26, Chapter 8]). If the bound \(\Vert u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\Vert _{L^2(\varOmega ;H^{1+\kappa _a}({\mathcal {D}}))}\le C_u=C_u(K,A)\) is finite for the fixed approximation parameters K, A, we immediately obtain
We note that, by construction of our random field, we always obtain an interface geometry with fixed angles and bounded jump height, which have great influence on the solution regularity, see e.g. [31].
For general elliptic jump-diffusion problems, one obtains a discretization error of order \(\kappa _a \in (1/2,1)\). In general, we cannot expect the full order of convergence \(\kappa _a=1\) for discontinuous diffusion coefficients. Without special treatment of the interfaces with respect to the triangulation, one cannot expect a convergence rate which is higher than \(\kappa _a=1/2\) for the deterministic interface problem (see [6, 11]). The convergence of the FE method may be improved by the use of triangulations which are adapted to the discontinuities. This is explained in more detail in the following subsection.
6.1 Sample-adapted triangulations
In [11], the authors suggest sample-adapted triangulations to improve the convergence rate of the FE approximation: Consider a fixed \(\omega \in \varOmega \) and assume that the discontinuities of the diffusion coefficient are described by the partition \({\mathcal {T}}(\omega )=({\mathcal {T}}_i,~i=1,\dots , \tau (\omega ))\) of the domain \({\mathcal {D}}\) where \(\tau (\omega )\) describes the number of elements in the partition. We consider finite-dimensional subspaces \({\hat{V}}_\ell (\omega )\subset V\) with (stochastic) dimension \({\hat{d}}_\ell (\omega )\in {\mathbb {N}}\). We denote by \({\hat{\theta }}_\ell (\omega )\) the minimal interior angle within \({\mathcal {K}}_\ell (\omega )\) and assume the existence of a positive number \(\theta >0\) such that \(\inf \{{\hat{\theta }}_\ell (\omega ) ~|~ \ell \in {\mathbb {N}}_0\}\ge \theta \) for \({\mathbb {P}}\)-almost all \(\omega \in \varOmega \). Assume that \({\mathcal {K}}_\ell (\omega )\) is a triangulation of \({\mathcal {D}}\) which is adjusted to the partition \({\mathcal {T}}(\omega )\) in the sense that for every \(i=1,\dots ,\tau (\omega )\) it holds
for all \(\ell \in {\mathbb {N}}_0\), where \(({\overline{h}}_\ell ,~\ell \in {\mathbb {N}}_0)\) is a deterministic, decreasing sequence of refinement thresholds which converges to zero (see Fig. 2).
A sample of a Poisson-subordinated Matérn-1.5-GRF (left) with corresponding sample-adapted triangulations (right)
Sample-adapted triangulations lead to an improved convergence rate for elliptic problems with discontinuous coefficients (see [11, Section 4.1] and Sect. 7). This is especially true for jump-diffusion coefficients with polygonal jump geometry, which is the case for the subordinated GRF considered in this paper (see Fig. 2, [11, 15]). The question arises whether mean squared convergence rates can be derived analytically (cf. Assumption 7.1). While it is not possible to prove convergence rates for the mean squared error for our diffusion coefficient in general due to the regularity of the stochastic solution, in practice one can (at least) recover the convergence rates of the deterministic jump-diffusion problem in the strong error (see Sect. 7, [11, Section 4.1] and [15, Section 2]). In fact, in the non-adapted case it is possible to get even better convergence rates than expected for some examples.
7 Numerical examples
In the following experiments we work on the domain \({\mathcal {D}}=(0,1)^2\) and use a FE method with hat-function basis as described in Sect. 6. Here, we distinguish between the standard FEM approach with standard triangulations and the sample-adapted FEM approach. The aim of our numerical experiments is to compare the sample-adapted with the non-adapted approach and investigate how different Lévy subordinators and GRFs influence the strong convergence rate. In our first example we use Poisson processes with low intensity to investigate the superiority of the presented sample-adapted triangulation. In the second example we use Poisson subordinators with a significantly higher intensity. Besides Poisson subordinators we also use Gamma processes which have infinite activity.
7.1 Strong error approximation
Remark 6.1 and Sect. 6.1 motivate the following assumption on the approximation error of the FE method for non-adapted and adapted triangulations for the rest of this paper (see [11, Assumption 4.4]).
Assumption 7.1
There exist deterministic constants \({\hat{C}}_{u}, C_{u},{\hat{\kappa }}_a,\kappa _a>0\) such that for any \(\varepsilon _W,\varepsilon _l>0\) and any \(\ell \in {\mathbb {N}}_0\), the Finite Element approximation errors of \({\hat{u}}_{K,A,\ell }^{(\varepsilon _W,\varepsilon _l)}\approx u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\) in the (sample-adapted) subspaces \({\hat{V}}_\ell \), respectively \(u_{K,A,\ell }^{(\varepsilon _W,\varepsilon _l)}\approx u_{K,A}^{(\varepsilon _W,\varepsilon _l)}\) in \(V_\ell \), are bounded by
where the constants \({\hat{C}}_{u},C_{u}\) may depend on a, f, g, K, A but are independent of \({\hat{h}}_\ell ,h_\ell ,{\hat{\kappa }}_a\) and \(\kappa _a\).
In each numerical experiment we choose a problem dependent cut-off level K for the subordinators in (4.4) large enough so that its influence is negligibly (see Remark 4.1). Further, we choose the cut-off level for the diffusion coefficient A in (4.8) large enough such that it has no influence in numerical experiments, \(e.g.~ A=50\) and therefore the error induced by the error contributor \(E_1\) in (5.4) can be neglected in our experiments. We estimate the strong error using a standard Monte Carlo estimator. Assume that a sequence of (sample-adapted) finite-dimensional subspaces \(({\hat{V}}_\ell ,~\ell \in {\mathbb {N}}_0)\subset V\) is given where we use the notation of Sect. 6. For readability we only treat the case of pathwise sample-adapted Finite Element approximations in the rest of the theoretical consideration in this subsection. We would like to point out, however, that similar arguments lead to the corresponding results for standard FE approximations.
Under the assumptions of Theorem 5.11 and Assumption 7.1 we obtain
with a constant \(C=C(C_{reg},D,{\overline{a}}_-,{\hat{C}}_u)\). Therefore, in order to equilibrate all error contributions, we choose the approximation parameters \(\varepsilon _W\) and \(\varepsilon _l\) in the following way:
For readability, we omit the cut-off parameters K and A in the following and use the notation \({\hat{u}}_{\ell ,\varepsilon _W,\varepsilon _l}={\hat{u}}_{K,A,\ell }^{(\varepsilon _W,\varepsilon _l)}\). Choosing the approximation parameters \(\varepsilon _W, \varepsilon _l\) according to (7.2), we can investigate the strong error convergence rate by a Monte Carlo estimation of the left hand side of (7.1): for a fixed natural number \(M\in {\mathbb {N}}\) we approximate
where \((u_{ref}^{(i)},~i=1,\dots , M)\) are i.i.d. realizations of the stochastic reference solution \(u_{ref}\approx u_{K,A}\) and \(({\hat{u}}_{\ell ,\varepsilon _W,\varepsilon _l}^{ (i)}, ~i=1\dots ,M)\) are i.i.d. realizations of the FE approximation \({\hat{u}}_{\ell ,\varepsilon _W,\varepsilon _l}\) of the PDE solution on the FE subspace \({\hat{V}}_\ell \). In all examples we choose the sample number M so that the standard deviation of the MC samples is smaller than 10% of the MC estimator itself.
7.2 PDE parameters
In all of our numerical examples we choose \({\overline{a}}\equiv 1/10\), \(f\equiv 10\), \(\Phi _1=1/100\,\exp (\cdot )\) and \(\Phi _2=5\,|\cdot |\). Further, we impose mixed Dirichlet–Neumann boundary conditions if nothing else is explicitly mentioned. To be precise, we split the domain boundary \(\partial {\mathcal {D}}\) by \(\varGamma _1=\{0,1\}\times [0,1]\) and \(\varGamma _2=(0,1)\times \{0,1\}\) and impose the pathwise mixed Dirichlet–Neumann boundary conditions
for \(\omega \in \varOmega \).
We choose \(W_1\) to be a Matérn-1.5-GRF on \({\mathcal {D}}\) with correlation length \(r_1=0.5\) and different variance parameters \(\sigma _1^2\). Further, we set \(W_2\) to be a Matérn-1.5-GRF on \([0,K]^2\) which is independent of \(W_1\) with different variances \(\sigma _2^2\) and correlation lengths \(r_2\). We use a reference grid with \(800\times 800\) equally spaced points on the domain \({\mathcal {D}}\) for interpolation and prolongation.
7.3 Poisson subordinators
In this section we use Poisson processes to subordinate the GRF \(W_2\) in the diffusion coefficient in (3.1). We consider both, high and low intensity Poisson processes and vary the boundary conditions. Further, using Poisson subordinators allows for a detailed investigation of the approximation error caused by approximating the Lévy subordinators \(l_1\) and \(l_2\) according to Assumption 4.2(v).
7.3.1 The two approximation methods
Using Poisson processes as subordinators allows for two different simulation approaches in the numerical examples: the first approach is an exact and grid-independent simulation of a Poisson process using for example the Method of Exponential Spacings or the Uniform Method (see [33, Section 8.1.2]). On the other hand, one may also work with approximations of the Poisson processes satisfying Assumption 4.2(v).
We recall that a Poisson(\(\lambda \))-process l is a Lévy process with \(l(t)\sim Poiss(t\lambda )\) for \(t>0\), i.e. it admits the discrete density
We sample the values of the Poisson(\(\lambda \))-processes \(l_1\) and \(l_2\) on an equidistant grid \(\{x_i,~i=0,\dots ,N_l\}\) with \(x_0=0\) and \(x_{N_l}=1\) and step size \(|x_{i+1}-x_i|\le \varepsilon _l\le 1\) for all \(i=0,\dots ,N_l-1\). Further, we approximate the stochastic processes by a piecewise constant extension \(l_j^{(\varepsilon _l)}\approx l_j\) of the values on the grid:
for \(j=1,2\). Since the Poisson process has independent, Poisson distributed increments, values of the Poisson process at the discrete points \(\{x_i,~i=0,\dots ,N_l\}\) can be generated by adding independent Poisson distributed random variables. In the following we refer to this approach as the approximation approach to simulate a Poisson process. Note that in this case Assumption 4.2(v) holds with \(\eta =+\infty \). In fact, for any \(s\in [1,+\infty )\) we obtain for \(j=1,2\) and an arbitrary \(x\in [0,1)\) with \(x\in [x_i,x_{i+1})\):
which is independent of the specific \(x\in [0,1)\). Note that this also holds for \(x=D=1\) and therefore
For a Poisson process with parameter \(\lambda \) we obtain
where the series converges by the ratio test.
Since the Poisson process allows for both approaches - approximation and exact simulation of the process - the use of these processes are suitable to investigate the additional error in the approximation of the PDE solution resulting from an approximation of the subordinators. Note that the main difference of the proposed approaches is that the approximation approach is grid dependent and the exact simulation of the Poisson process is not. The differences in the computational costs for the two approaches are insignificant in our numerical examples.
7.3.2 Poisson subordinators: low intensity and mixed boundary conditions
In this example we choose \(l_1\) and \(l_2\) to be Poisson(1)-subordinators. Further, the variance parameter of the GRF \(W_1\) is set to be \(\sigma _1^2=1.5^2\) and the variance and correlation parameters of the GRF \(W_2\) are given by \(\sigma _2^2=0.3^2\) and \(r_2=1\).
For independent Poisson(1)-subordinators \(l_1\) and \(l_2\) we choose \(K=8\) as the cut-off parameter [see (4.4)]. With this choice we obtain
for \(j=1,2\), such that this cut-off has no influence in the numerical example. Note that for Matérn-1.5-GRFs we can expect \(\gamma =1\) in Eq. (4.9) (see e.g. [12, Chapter 5], [23, Proposition 9] and [14, Section 2.3]).
We approximate the GRFs \(W_1\) and \(W_2\) by the circulant embedding method (see [24, 25]) to obtain approximations \(W_1^{(\varepsilon _W)}\approx W_1\) and \(W_2^{(\varepsilon _W)}\approx W_2\) as in Lemma 4.4. Since \(\eta =+\infty \) and \(f\in L^q(\varOmega ;H)\) for every \(q\ge 1\) we choose for any positive \(\delta >0\)
to obtain from Theorem 5.11
where we have to assume that \(j_{reg}\ge 2((1+\delta )/\delta -1)\) and \(k_{reg}\ge 2(1+\delta )/\delta \) for the regularity constants \(j_{reg},k_{reg}\) given in Assumption 5.7. For \(\delta =0.05\) we obtain
Therefore, we get \(\gamma =1\) and \(c=2.01/2\) in the equilibration formula (7.2).
Figure 3 shows two different samples of the diffusion coefficient and the corresponding FE approximations of the PDE solution.
Different samples of the diffusion coefficient with Poisson(1)-subordinators and the corresponding PDE solutions with mixed Dirichlet–Neumann boundary conditions
The FE discretization parameters are given by \(h_\ell = 0.4\cdot 2^{-(\ell -1)}\) for \(l=1,\dots ,7\). We set \(u_{ref}={\hat{u}}_{7,\varepsilon _W,\varepsilon _l}\), where the approximation parameters \(\varepsilon _W\) and \(\varepsilon _l\) are choosen according to (7.2) and we compute \(M=100\) samples to estimate the strong error by the Monte Carlo estimator [see (7.3)]. In this experiment we investigate the strong error convergence rate for the sample-adapted FE approach as well as convergence rate for the non-adapted FE approach (see Sect. 6). In Sect. 7.3.1 we described two approaches to simulate Poisson subordinators. We run this experiment with both approaches: first, we approximate the Poisson process via sampling on an equidistant level-dependent grid and, in a second run of the experiment, we simulate the Poisson subordinators exactly using the Uniform Method described in [33, Section 8.1.2]. The convergence results for the both approaches for this experiment are given in the Fig. 4.
Convergence results for Poisson(1)-subordinators using the approximation approach and the Uniform Method with mixed Dirichlet–Neumann boundary conditions
We see a convergence rate of approximately 0.7 for the standard FEM discretization and full order convergence (\({\hat{\kappa }}_a\approx 1\)) for the sample-adapted approach. On the right hand side of Fig. 4 one sees that the sample-adapted approach is more efficient in terms of computational effort if we consider the error-to-(averaged)DOF-plot. Only on the first level the standard FEM approach seems to be more efficient (pre-asymptotic behaviour). If we compare the results for the approximation method with the Uniform Method (see Sect. 7.3.1), we find that, while the convergence rates are the same, the constant of the error in the sample-adapted approach is slightly smaller for the Uniform Method. This shift is exactly the additional error resulting from an approximation of the subordinators in the approximation approach. We also see that, compared to the approximation approach, on the lower levels the averaged degrees of freedom in the sample-adapted FEM approach is slightly higher if we simulate the Poisson subordinators exactly. This is caused by the fact that in this case we do not approximate the discontinuities of the field which are generated by the Poisson processes. This results in a higher average number of degrees of freedom on the lower levels because discontinuities are more likely close to each other.
7.3.3 Poisson subordinators: low intensity and homogeneous Dirichlet boundary conditions
Next, we consider the elliptic PDE under homogeneous Dirichlet boundary conditions. All other parameters remain as in Sect. 7.3.2. Figure 5 shows samples of the diffusion coefficient and the corresponding FE approximation of the PDE solution.
Different samples of the diffusion coefficient with Poisson(1)-subordinators and the corresponding PDE solutions with homogeneous Dirichlet boundary conditions
We estimate the strong error convergence rate for this problem in the same way as in the previous example using \(M=250\) samples and we use the approximation approach to simulate the Poisson subordinators (see Sect. 7.3.1). Convergence results are given in Fig. 6.
Convergence results for Poisson(1)-subordinators using the approximation approach and homogeneous Dirichlet boundary conditions
As in the experiment with mixed Dirichlet-Neumann boundary conditions we obtain convergence order of \(\kappa _a\approx 0.7\) for the standard FEM approach and full order convergence for the sample-adapted approach. Also in case of homogeneous Dirichlet boundary conditions the sample-adapted FEM is more efficient in terms of the averaged number of degrees of freedom.
7.3.4 Poisson subordinators: high intensity and mixed boundary conditions
In this section we want to consider subordinators with higher intensity, resulting in a higher number of discontinuities in the diffusion coefficient. Therefore, we consider \(l_1\) and \(l_2\) to be Poisson(5) processes. We set \(K=1\) in the diffusion coefficient (4.4) and consider the downscaled process
for \(t\in [0,1]\) and \(j=1,2\). With this choice it is reasonable to expect that this cut-off has no numerical influence since
for \(j=1,2\). The reason we consider the downscaled process is that otherwise we would have to simulate the GRF \(W_2\) on the domain \([0,15]^2\) which is very time consuming. Note that the downscaling of the subordinators has no effect on the jump activity. The variance parameter of the field \(W_1\) is chosen to be \(\sigma _1^2=1\) and the parameters of the GRF \(W_2\) are set to be \(\sigma _2^2=0.3^2\) and \(r_2=0.5\). Figure 7 shows samples of the coefficient and the corresponding pathwise FEM solution.
Different samples of the diffusion coefficient with Poisson(5)-subordinators (top) and the corresponding PDE solutions with mixed Dirichlet–Neumann boundary conditions (bottom)
As in the first experiment, we again run this experiment using both methods described in Sect. 7.3.1: the approximation approach using Poisson-distributed increments and the Uniform Method. We use the discretization steps \(h_\ell = 0.1\cdot 1.7^{-(\ell -1)}\) for \(\ell = 1,\dots ,7\) and \(M=150\) samples.
Convergence results for Poisson(5)-subordinators using the approximation approach and the Uniform Method with mixed Dirichlet–Neumann boundary conditions
In Fig. 8 we see that we get almost full order convergence for the sample-adapted FE method for both approximation approaches of the Poisson processes. Compared to the low-intensity examples with Poisson(1)-subordinators given in Sects. 7.3.2 and 7.3.3, we get a slightly lower convergence rate of approximately 0.55 for the standard FEM approach. This holds for both approximation methods of the Poisson subordinators. Hence, we see that the way how the Poisson-subordinators are simulated seems to have no effect on the convergence rate.
7.3.5 Poisson subordination of a GRF with short correlation length: high intensity and mixed boundary conditions
In our construction of the jump-diffusion coefficient, the jumps are generated by the subordinated GRF. To be precise, the number of spatial jumps is determined by the subordinators and the jump intensities (in terms of the differences in height between the jumps) are essentially determined by the GRF \(W_2\). This fact allows to control the jump intensities of the diffusion coefficient by the correlation parameter of the underlying GRF \(W_2\). In the following experiment we want to investigate the influence of the jump intensities of the diffusion coefficient on the convergence rates.
In Sect. 7.3.4 we subordinated a Matérn-1.5-GRF with pointwise standard deviation \(\sigma _2^2= 0.3^2\) and a correlation length of \(r_2=0.5\). In the following experiment we set the standard deviation of the GRF \(W_2\) to \(\sigma _2^2=0.5^2\) and the correlation length to \(r_2=0.1\) and leave all the other parameters unchanged. Figure 9 compares the resulting GRF with the field \(W_2\) with parameters \(\sigma _2^2=0.3^2\) and \(r_2=0.5\) which we used in Sect. 7.3.4.
Samples of a Matérn-1.5-GRF with \(\sigma _2^2=0.3^2\), \(r_2=0.5\) (left) and with parameters \(\sigma _2^2=0.5^2\), \(r_2=0.1\) (right)
Subordinating the GRF with small correlation length (right plots in Fig. 9) result in higher jump intensities in the diffusion coefficient as the subordination of the GRF with higher correlation length (left plots in Fig. 9). Figure 10 shows samples of the diffusion coefficient and the corresponding PDE solutions where the parameters of \(W_2\) are \(\sigma _2^2=0.5^2\) and \(r_2=0.1\).
Different samples of the diffusion coefficient with Poisson(5)-subordinators and the corresponding PDE solutions with mixed Dirichlet-Neumann boundary conditions and small correlation length \(r_2=0.1\) of the GRF \(W_2\)
As expected, the resulting jump coefficient shows jumps with a higher intensity compared to the jump coefficient in the previous experiment where we used the GRF \(W_2\) with parameters \(\sigma _2^2=0.3^2\) and \(r_2=0.5\) (see Fig. 7).
We estimate the strong error convergence rate using this high-intensity jump coefficient using \(M=200\) samples and approximate the Poisson subordinators by the Uniform Method.
Convergence results for Poisson(5)-subordinators using the Uniform Method with mixed Dirichlet–Neumann boundary conditions and GRF parameters \(\sigma _2^2=0.5^2\) and \(r_2=0.1\)
Figure 11 shows that for the GRF \(W_2\) with small correlation length the convergence rates are reduced for both approaches: the standard FEM approach and the sample-adapted version. We cannot preserve full order convergence in the sample-adapted FEM but observe a convergence rate of approximately 0.75. In the non-adapted approach we obtain a convergence rate of approximately 0.45. Looking at the error-to-(averaged)DOF-plot on the right hand side of Fig. 11 we see that still the sample-adapted approach is by a large margin more efficient in terms of computational effort. This experiment confirms our expectations since the FEM convergence rate has been shown to be strongly influenced by the regularity of the jump-diffusion coefficient (see e.g. [11, 31]).
7.4 Gamma subordinators
In order to also consider Lévy subordinators with infinite activity we take Gamma processes to subordinate the GRF in the remaining numerical examples. We set the standard deviation of the GRF \(W_1\) to be \(\sigma _1^2=1.5^2\) and we choose \(\sigma _2^2=0.3^2\) and \(r_2=1\) for the Matérn-1.5-GRF \(W_2\) and leave the other parameters unchanged. For \(a_G,b_G>0\), a \(Gamma(a_G,b_G)\)-distributed random variable admits the density function
where \(\varGamma (\cdot )\) denotes the Gamma function. A Gamma process \((l(t), t\ge 0)\) has independent Gamma distributed increments. Being precise, \(l(t)-l(s)\sim Gamma(a_G\cdot (t-s),b_G)\) for \(0<s<t\) (see [33, Chapter 8]). The following lemma is essential to approximate the Gamma processes.
Lemma 7.2
Let Z be a Gamma\((a_G,b_G)\) distributed random variable for positive parameters \(a_G,b_G>0\). It holds
for all \(n\in {\mathbb {N}}_0\).
Proof
We calculate
where we used that the integral over the Gamma density equals one in the last step. \(\square \)
In our numerical experiments we choose \(l_j\) to be a Gamma(4, 10)-process for \(j=1,2\). Since increments of a Gamma process are Gamma-distributed random variables it is straightforward to generate values of a Gamma process on grid points \((x_i)_{i=0}^{N_l}\subset [0,1]\) with \(|x_{i+1}-x_i|\le \varepsilon _l\) for \(i=0,\dots ,N_l-1\). We then use the piecewise constant extension of the simulated values \(\{l_j(x_i),~i=0,\dots ,N_l-1,~j=1,2\}\) to approximate the Lévy subordinators:
for \(j=1,2\). Note that in this case Assumption 4.2(v) is fulfilled for any fixed \(\eta <+\infty \). To see that we consider a fixed \(s\in {\mathbb {N}}\) with \(s\le \eta \) and calculate for an arbitrary \(x\in [0,1)\) with \(x\in [x_i,x_{i+1})\):
Figure 12 shows samples of the jump-diffusion coefficient with Gamma(4, 10)-subordinator and corresponding FE solution where we used mixed Dirichlet-Neumann boundary conditions.
Different samples of the diffusion coefficient with Gamma(4, 10)-subordinators and the corresponding PDE solutions with mixed Dirichlet–Neumann boundary conditions
We set the diffusion cut-off to \(K=2\) since in this case we obtain
for \(j=1,2\). The use of infinite-activity Gamma subordinators in the diffusion coefficient does not allow anymore for a sample-adapted approach to solve the PDE problem. Hence, we only use the standard FEM approach to solve the PDE samplewise and estimate the strong error convergence. We use \(M=200\) samples to estimate the strong error on the levels \(\ell = 1,\dots ,5\) where we set the non-adaptive FEM solution \(u_{7,\varepsilon _W,\varepsilon _l}\) on level \(L=7\) to be the reference solution. We choose the FEM discretization steps to be \(h_\ell = 0.4\cdot 2^{-(\ell -1)}\) for \(\ell = 1,\dots ,7\).
Convergence results for Gamma(4, 10)-subordinators with mixed Dirichlet–Neumann boundary conditions
Figure 13 shows a convergence rate of approximately 0.8 for the standard-FEM approach. Since we do not treat the discontinuities in a special way we cannot expect full order convergence. In fact, the given convergence is comparably good since in general we cannot prove a higher convergence order than 0.5 for the standard deterministic FEM approach without special treatment of the discontinuities (see [6, 11]). The convergence rate of approximately 0.8 in this example is based on the comparatively large correlation length of the underlying GRF \(W_2\) (see Sect. 7.3.5).
In Sect. 7.3.5 we investigated the effect of a rougher diffusion coefficient on the convergence rate for Poisson(5)-subordinators. In the following experiment we follow a similar strategy and use a shorter correlation length in the GRF \(W_2\) which is subordinated by Gamma processes. Therefore, we choose the parameters of the Matérn-1.5-GRF \(W_2\) to be \(\sigma _2^2=0.3^2\) and \(r_2=0.05\). Figure 14 shows a comparison of the resulting GRFs \(W_2\) with the different correlation lengths.
Samples of a Matérn-1.5-GRF with correlation length \(r_2=1\) (left) and with correlation length \(r_2=0.05\) (right)
In Fig. 15, the GRF with small correlation length results in higher jumps of the diffusion coefficient and stronger deformations of the corresponding PDE solution compared to the previous example (see Fig. 12).
Different samples of the diffusion coefficient with Gamma(4, 10)-subordinators and the corresponding PDE solutions with mixed Dirichlet–Neumann boundary conditions where the correlation length of \(W_2\) is \(r_2=0.05\)
We estimate the strong error taking \(M=200\) samples where we use the non-adapted FEM solution \(u_{9,\varepsilon _W,\varepsilon _l}\) on level \(L=9\) as reference solution and choose the FEM discretization steps to be \(h_\ell = 0.1\cdot 1.5^{-(\ell -1)}\) for \(\ell = 1,\dots ,9\). Figure 16 shows the convergence on the levels \(\ell =1,\dots ,6\).
Convergence results for Gamma(4, 10)-subordinators with mixed Dirichlet–Neumann boundary conditions where the correlation length of \(W_2\) is \(r_2=0.05\)
We observe a convergence rate of approximately 0.45 which is significantly smaller than the rate of approximately 0.8 we obtained in the example where we used a GRF \(W_2\) with correlation length \(r_2=1\) (see Fig. 13). This again confirms that, for subordinated GRFs, the convergence rate of the FE method is highly dependent on the correlation length of the underlying GRF \(W_2\) and the resulting jump-intensity of the diffusion coefficient.
Availability of data and materials
All data used in our experiments have been produced with MATLAB random number generators and no external datasets have been used. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.
Notes
For simplicity we chose a square domain, recangular ones may be considered in the same way
For simplicity we assume one fixed K for all spatial dimensions. The results in the subsequent sections hold for individual independent values in each spatial dimension as well.
The error of the approximation \(u_{K} \approx u\) may be controlled as in Remark 4.1 but cannot be quantified for the solution.
References
Abdulle, A., Barth, A., Schwab, C.: Multilevel Monte Carlo methods for stochastic elliptic multiscale PDEs. Multiscale Model. Simul. 11(4), 1033–1070 (2013). https://doi.org/10.1137/120894725
Adams, R.A., Fournier, J.J.F.: Sobolev Spaces, Pure and Applied Mathematics (Amsterdam), vol. 140, 2nd edn. Elsevier/Academic Press, Amsterdam (2003)
Adler, R.J., Taylor, J.E.: Random Fields and Geometry. Springer Monographs in Mathematics. Springer, New York (2007)
Aliprantis, C.D., Border, K.C.: Infinite Dimensional Analysis, 3rd edn. Springer, Berlin (2006). A Hitchhiker’s Guide
Applebaum, D.: Lévy Processes and Stochastic Calculus, Cambridge Studies in Advanced Mathematics, vol. 116, 2nd edn. Cambridge University Press, Cambridge (2009). https://doi.org/10.1017/CBO9780511809781
Babuška, I.: The finite element method for elliptic equations with discontinuous coefficients. Computing (Arch. Elektron. Rechnen) 5, 207–213 (1970). https://doi.org/10.1007/bf02248021
Babuška, I., Tempone, R., Zouraris, G.E.: Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM J. Numer. Anal. 42(2), 800–825 (2004). https://doi.org/10.1137/S0036142902418680
Barth, A., Merkle, R.: Subordinated gaussian random fields. ArXiv e-prints, arXiv:2012.06353 [math.PR] (2020)
Barth, A., Schwab, C., Zollinger, N.: Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients. Numer. Math. 119(1), 123–161 (2011). https://doi.org/10.1007/s00211-011-0377-0
Barth, A., Stein, A.: Approximation and simulation of infinite-dimensional Lévy processes. Stoch. Part. Differ. Equ. Anal. Comput. 6(2), 286–334 (2018). https://doi.org/10.1007/s40072-017-0109-2
Barth, A., Stein, A.: A study of elliptic partial differential equations with jump diffusion coefficients. SIAM/ASA J. Uncertain. Quantif. 6(4), 1707–1743 (2018). https://doi.org/10.1137/17M1148888
Barth, A., Stein, A.: A multilevel Monte Carlo algorithm for parabolic advection-diffusion problems with discontinuous coefficients. In: Springer Proceedings in Mathematics and Statistics, pp. 445–466. Springer (2020). https://doi.org/10.1007/978-3-030-43465-6-22
Charrier, J.: Strong and weak error estimates for elliptic partial differential equations with random coefficients. Hyper Articles en Ligne, INRIA (2010). http://hal.inria.fr/inria-00490045/en/
Charrier, J., Scheichl, R., Teckentrup, A.L.: Finite element error analysis of elliptic PDEs with random coefficients and its application to multilevel Monte Carlo methods. SIAM J. Numer. Anal. 51(1), 322–352 (2013). https://doi.org/10.1137/110853054
Chen, Z., Zou, J.: Finite element methods and their convergence for elliptic and parabolic interface problems. Numer. Math. 79(2), 175–202 (1998). https://doi.org/10.1007/s002110050336
Da Prato, G., Zabczyk, J.: Stochastic Equations in Infinite Dimensions, Encyclopedia of Mathematics and its Applications, vol. 152, 2nd edn. Cambridge University Press, Cambridge (2014). https://doi.org/10.1017/CBO9781107295513
Di Nezza, E., Palatucci, G., Valdinoci, E.: Hitchhiker’s guide to the fractional Sobolev spaces. ArXiv e-prints, arXiv:1104.4345v3 [math.FA] (2011)
Ding, Z.: A proof of the trace theorem of Sobolev spaces on Lipschitz domains. Proc. Am. Math. Soc. 124(2), 591–600 (1996). https://doi.org/10.1090/S0002-9939-96-03132-2
Dreyfuss, P.: Higher integrability of the gradient in degenerate elliptic equations. Potential Anal. 26(2), 101–119 (2007). https://doi.org/10.1007/s11118-006-9030-4
Ernst, O.G., Gottschalk, H., Kalmes, T., Kowalewitz, T., Reese, M.: Integrability and approximability of solutions to the stationary diffusion equation with lévy coefficient. ArXiv e-prints, arXiv:2010.14912v2 [math.AP] (2021)
Evans, L.C.: Partial Differential Equations. Graduate Studies in Mathematics, vol. 19, 2nd edn. American Mathematical Society, Providence (2010). https://doi.org/10.1090/gsm/019
Frauenfelder, P., Schwab, C., Todor, R.A.: Finite elements for elliptic problems with stochastic coefficients. Comput. Methods Appl. Mech. Eng. 194(2–5), 205–228 (2005). https://doi.org/10.1016/j.cma.2004.04.008
Graham, I.G., Kuo, F.Y., Nichols, J.A., Scheichl, R., Schwab, C., Sloan, I.H.: Quasi-Monte Carlo finite element methods for elliptic PDEs with lognormal random coefficients. Numer. Math. 131(2), 329–368 (2015). https://doi.org/10.1007/s00211-014-0689-y
Graham, I.G., Kuo, F.Y., Nuyens, D., Scheichl, R., Sloan, I.H.: Analysis of circulant embedding methods for sampling stationary random fields. SIAM J. Numer. Anal. 56(3), 1871–1895 (2018). https://doi.org/10.1137/17M1149730
Graham, I.G., Kuo, F.Y., Nuyens, D., Scheichl, R., Sloan, I.H.: Circulant embedding with QMC: analysis for elliptic PDE with lognormal coefficients. Numer. Math. 140(2), 479–511 (2018). https://doi.org/10.1007/s00211-018-0968-0
Hackbusch, W.: Elliptic Differential Equations, Springer Series in Computational Mathematics, vol. 18, 2nd edn. Springer, Berlin (2017). https://doi.org/10.1007/978-3-662-54961-2. Theory and numerical treatment
Iglesias, M.A., Lin, K., Stuart, A.M.: Well-posed Bayesian geometric inverse problems arising in subsurface flow. Inverse Probl. 30(11), 114001, 39 (2014). https://doi.org/10.1088/0266-5611/30/11/114001
Iglesias, M.A., Lu, Y., Stuart, A.: A Bayesian level set method for geometric inverse problems. Interfaces Free Bound. 18(2), 181–217 (2016). https://doi.org/10.4171/IFB/362
Li, J., Wang, X., Zhang, K.: Multi-level Monte Carlo weak Galerkin method for elliptic equations with stochastic jump coefficients. Appl. Math. Comput. 275, 181–194 (2016). https://doi.org/10.1016/j.amc.2015.11.064
Mugler, A., Starkloff, H.J.: On the convergence of the stochastic Galerkin method for random elliptic partial differential equations. ESAIM Math. Model. Numer. Anal. 47(5), 1237–1263 (2013). https://doi.org/10.1051/m2an/2013066
Petzoldt, M.: Regularity results for Laplace interface problems in two dimensions. Z. Anal. ihre Anwend. 20(2), 431–455 (2001). https://doi.org/10.4171/zaa/1024
Sato, K.i.: Lévy processes and Infinitely Divisible Distributions, Cambridge Studies in Advanced Mathematics, vol. 68. Cambridge University Press, Cambridge (2013). Translated from the 1990 Japanese original, Revised edition of the 1999 English translation
Schoutens, W.: Lévy Processes in Finance: Pricing Financial Derivatives. Wiley Series in Probability and Statistics. Wiley (2003). https://books.google.de/books?id=Be4NqIz7h-kC
Teckentrup, A.L., Scheichl, R., Giles, M.B., Ullmann, E.: Further analysis of multilevel Monte Carlo methods for elliptic PDEs with random coefficients. Numer. Math. 125(3), 569–600 (2013). https://doi.org/10.1007/s00211-013-0546-4
Werner, D.: Funktionalanalysis. Springer-Lehrbuch. Springer, Berlin, Heidelberg (2011). https://books.google.de/books?id=580NcHHrAXsC
Acknowledgements
Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—EXC 2075—390740016. The authors would like to thank two anonymous referees for their remarks, which led to a significant improvement of the manuscript.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Proof of Theorem 5.5
Proof of Theorem 5.5
Theorem A.1
(Theorem 5.5) For any \(\delta >0\) and any \(n\in (1,+\infty )\) there exists a constant \(A=A(\delta ,n)>0\) such that
Proof
Step 1: Tail estimation for the coefficient
By Assumption 4.2(ii), the functions \(\Phi _1,\Phi _2\) from Definition 3.4 fulfill
for \(x,y\in {\mathbb {R}}\). By the mean value theorem, for any \(x\in {\mathbb {R}}\) there exists a real number \(\xi \in {\mathbb {R}}\) with \(|\xi |\le |x|\) such that it holds
for positive constants \({\tilde{\phi }}\) and \({\tilde{\psi }}\) which are independent of \(x\in {\mathbb {R}}\).
Since the GRFs \(W_1\) and \(W_2\) are \({\mathbb {P}}-a.s.\) bounded on \({\mathcal {D}}\) resp. on \([0,K]^2\) it follows from Adler and Taylor [3, Theorem 2.1.1] that \(\mu _1:={\mathbb {E}}\left( \underset{(x,y)\in {\mathcal {D}}}{\sup }\,W_1(x,y)\right) <+\infty \) and
for \(m>\mu _1\) with a finite constant \(\sigma _{\mathcal {D}}^2\) defined by
by Assumption 4.2i. For a given \(\varepsilon \in (0,1)\) we choose the real number A such that it holds
With this choice we obtain the bound
This can be seen by the following calculation
where we used (A.1) in the first step, the estimate (A.2) in the third step and condition (A.3) in the last step.
Obviously, an estimation as in Eq. (A.2) holds for the GRF \(W_2\):
for \(m>\mu _2\) with \(\mu _2:={\mathbb {E}}\left( \underset{(x,y)\in [0,K]^2}{\sup }\,W_2(x,y)\right) <+\infty \) and
By the Lipschitz continuity of \(\Phi _2\) we conclude the existence of a constant \(\phi _2>0\) such that
for \(x\in {\mathbb {R}}\). If we again fix a positive \(\varepsilon \in (0,1)\) and choose the real number A such that
we obtain the following bound
This can be seen by the following calculation:
Step 2: Finite moments of the coefficient
In this step we want to show that for any \(n\in [1,+\infty )\) it holds
We use the definition of the coefficient \(a_K\) in (4.4) and Hölder’s inequality to calculate
Therefore, it remains to show that it holds \(I_1,I_2<+\infty \).
By Fubini’s theorem, for every nonnegative random variable X it holds
if the right hand side exists. We use this fact and Eq. (A.1) to estimate for \(I_1\):
where we split the integral and used Eq. (A.2) in the last step. In a similar way, we use Eq. (A.6) to calculate for the second summand \(I_2\):
where we used Eq. (A.5) in the last step. This proves Eq. (A.8).
Step 3: Estimate for the approximation of the diffusion coefficient.
Now, let \(\delta \in (0,1)\) be arbitrary. Choose \(A=A(\delta )>0\) such that
for \(\varepsilon :=1-\sqrt{1-(\delta ^{2s} /C_{a_K}(2s,K,D))}\).
We estimate using Hölder’s inequality:
For the second factor we estimate using the independence of \(W_1\) and \(W_2\)
where we used Eqs. (A.4) and (A.7) in the fourth step and therefore we obtain
\(\square \)
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Merkle, R., Barth, A. Subordinated Gaussian random fields in elliptic partial differential equations. Stoch PDE: Anal Comp 11, 819–867 (2023). https://doi.org/10.1007/s40072-022-00246-w
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s40072-022-00246-w