Abstract
We study the stochastic p-Laplace system in a bounded domain. We propose two new space–time discretizations based on the approximation of time-averaged values. We establish linear convergence in space and 1/2 convergence in time. Additionally, we provide a sampling algorithm to construct the necessary random input in an efficient way. The theoretical error analysis is complemented by numerical experiments.
Similar content being viewed by others
1 Introduction
Let \({\mathcal {O}} \subset {\mathbb {R}}^n\) be a polygonal domain, \(n \ge 1\), \(N \ge 1\), \(T > 0\) be finite. We are interested in the approximation of the solution process \(u :\Omega \times [0,T] \times {\mathcal {O}} \rightarrow {\mathbb {R}}^N\) to the stochastic p-Laplace system. Given an initial datum \(u_0\) and a stochastic forcing term \(G(u) \,\textrm{d}W\) (for the precise assumptions see Assumption 1), u is determined by the relations
where \(S(\xi ):= \left( \kappa + \left|\xi \right| \right) ^{p-2} \xi \in {\mathbb {R}}^{N \times n}\), \(p \in (1,\infty )\) and \(\kappa \ge 0\). Closely related to S is the nonlinear operator \(V(\xi ):= \left( \kappa + \left|\xi \right| \right) ^\frac{p-2}{2} \xi \in {\mathbb {R}}^{N \times n}\).
1.1 Model
The system (1.1) has many important applications in nature. As one major example, it provides a prototype system towards the modeling of non-Newtonian fluids. More specifically, it is closely related to power law fluids [14, 69]. The case \(p=2\) corresponds to the famous stochastic heat equation and has been studied extensively, both analytically e.g. [12, 47, 48] as well as numerically e.g. [1, 24, 36, 42, 67].
The p-Laplace system arises as a stochastically perturbed gradient flow of the energy \(J : W^{1,p}_0({\mathcal {O}}) \rightarrow [0,\infty )\) defined by
where \(\varphi (t):= \int _0^t (\kappa +s)^{p-2}s \,\textrm{d}s\). In the case \(\kappa = 0\) the energy (1.2) corresponds to the classical p-Dirichlet energy.
1.2 Existence and well-posedness
The existence of analytically weak solutions in the space \( L^2\left( \Omega ; C\left( [0,T]; L^2\left( {\mathcal {O}} \right) \right) \right) \cap L^p\big (\Omega ;L^p\big (0,T;W^{1,p}_0\left( {\mathcal {O}}\right) \big )\big )\) can be established by standard monotonicity arguments [61]. It requires a linear growth assumption on the noise coefficient G and \(L^2({\mathcal {O}})\)-integrable initial data.
First results on the existence of strong solutions to stochastically perturbed gradient flows have been obtained by Gess [43]. It includes the degenerate, \(p \ge 2\), p-Laplace equation.
In the literature different generalization of (1.1) have been considered. Well-posedness for merely \(L^1\)-initial data has been addressed in [68]. More general systems, where p is allowed to depend on \((\omega ,t,x)\) respectively on (t, x), are considered in [73, 74] respectively in [9]. It was further extended by Breit and Gmeineder [17] to electro-rheological fluids. Electro-rheological fluids are modeled by (1.1) complemented by an additional divergence-free constraint and a free pressure variable. The singular case \(p \in [1,2)\) has been analyzed in [60, 44] and [5].
1.3 Regularity
Regularity properties of a function u become particularly important when it comes to the approximability of u within a discrete function class. Prominently, discrete tensor spaces generated via a time stepping scheme and a finite element discretization in space can approximate smooth functions more easily compared to non-smooth functions.
Historically, many authors addressed Hölder and \(C^{1,\alpha }\)-regularity of solutions to the deterministic steady p-Laplace equation, e.g. [26, 34, 56, 64, 70,71,72]. In general \(\alpha \in (0,1)\) is an unknown quantity. While \(C^{1,\alpha }\)-regularity can be used to measure the approximation quality of finite elements, they usually fail to produce optimal results, since in general \(\alpha < 1\).
Sobolev regularity provides an alternative scale of smoothness. In the singular case, \(p \in (1,2)\), \(W^{2,2}\)-regularity has been proven by Liu and Barrett [59]. The degenerate setting is more delicate and one can not expect \(\nabla ^2 u\) to be well-defined on \({\{{\nabla u = 0}\}}\). In fact, due to the nonlinear structure of the equation, solutions have limited regularity even for smooth data. A sharp result in the 2-dimensional setting about limited regularity for the Hölder as well as Sobolev scale has been obtained by Iwaniec and Manfredi [51].
The nonlinear character of the equation naturally introduces the additional quantities \(S(\nabla u)\) and \(V(\nabla u)\). In the scalar case both expressions are robust with respect to \(p \in (1,\infty )\). Here it is possible to prove \(V(\nabla u) \in W^{1,2}\) via a difference quotient technique [13, 58] and \(S(\nabla u) \in W^{1,2}\) by a functional inequality [20]. The result \(V(\nabla u) \in W^{1,2}\) generalizes to the vectorial case. However, the functional inequality fails in the vectorial case at least point-wise for \(p \le 2(2-\sqrt{2})\approx 1.1715\). Therefore, it is unclear whether a regularity result \(S(\nabla u) \in W^{1,2}\) is achievable for small p. Nevertheless, for \(p > 2(2-\sqrt{2})\) it is shown in [2] that \(S(\nabla u) \in W^{1,2}\). Regularity for \(S(\nabla u)\) on the Besov and Triebel-Lizorkin scale in the plane for \(p \ge 2\) has been obtained in [3]. Estimates of \(S(\nabla u)\) in terms of Riesz potentials were derived in [53, 54].
Regularity results for the parabolic p-Laplace system were derived in [41, Theorem 6.2.1] (cf. [57] for \(p \ge 2\)). It was shown that \(u \in C^1_{\textrm{weak}}\big (0,T;L^2({\mathcal {O}})\big ) \cap C_{\textrm{weak}}\big ( 0,T;W^{1,p}_0({\mathcal {O}})\big )\). It was extended in [30] to the nonlinear tensor \(V(\nabla u)\). The authors showed, by formally testing the equation with \(-\Delta u\) respectively \(-\partial _t^2 u\), that
Additionally, \(S(\nabla u) \in L^2(0,T;W^{1,2}({\mathcal {O}}))\) was proven in [22] for either \(p\in (1,\infty )\) and \(N=1\) or \(p \ge 2\) and \(N > 1\).
Within the context of the stochastic p-Laplace system it is possible to prove similar spatial regularity as in (1.3). The formal testing needs to be replaced by a suitable application of Itô’s formula as done by Breit [15]. However, the time regularity is limited due to the presence of the stochastic forcing. In the super-quadratic regime partial time regularity can be recovered by exploiting the strong formulation of the p-Laplace system as done by Wichmann [76]. Overall, for appropriate data assumptions it is possible to verify (see Sect. 2.5)
Here \(\varPhi _2(s) := e^{s^2}-1\) and B denotes a Besov-Orlicz space (see Sect. 2.1).
1.4 Approximation
In the past many authors have studied the numerical approximation of the deterministic counterpart of (1.1), e.g. [6,7,8, 10, 11, 16, 19, 25, 30, 33, 35, 37, 75]. The error of the discrete and analytic solution has been expressed in various different quantities. It turned out that the natural quantity to measure the error is given by
where \(V(\xi ) := \left( \kappa + \left|\xi \right| \right) ^\frac{p-2}{2}\xi \) and \(u_{m,h}\) is an approximation of \(u(t_m)\). In [33] it has been proved that the expression \(\left\Vert V(\nabla u) - V(\nabla u_h)\right\Vert _{L^2({\mathcal {O}})}^2\) is equivalent to the energy error \({\mathcal {J}}(u_h) - {\mathcal {J}}(u)\). Here \(u_h\) is a minimizer of the energy on a subspace \(V_h \subset W^{1,p}_0\).
In the steady case, starting with the seminal work by Barrett and Liu [6] and further improvements in [37] and [29], it has been proved that
This settles the question about optimal convergence for piece-wise linear continuous elements. In fact, the paper [30] deals with the parabolic system and optimal convergence under the regularity assumption (1.3) has been achieved, i.e.,
The error quantity (1.6) relies on the fact that the mapping \(t \mapsto \left\Vert V(\nabla u(t))\right\Vert _{L^2({\mathcal {O}})}\) is continuous. However, if the data is not sufficiently smooth, point-values might not be well-defined. In general, dealing with irregular data and therefore irregular solutions is a delicate task. Different methods have been suggested to recover well-defined error quantities.
In [16] the first and the third author develop a numerical scheme based on time averages to circumvent the usage of point evaluations. A more probabilistic approach has been used in [38]. There the authors replace deterministic evaluation points by random ones.
First results for monotone stochastic equations were derived by Gyöngy and Millet [45, 46]. They developed an abstract discretization theory that also covers the system (1.1) in the superquadratic case. In [45] they gave conditions that guarantee plain convergence of their discretization even for the full degenerate, \(\kappa = 0\), system. Contrary, convergence rates have been obtained under a non-degeneracy assumption, \(\kappa > 0\), in [46]. In both cases, they require restrictive assumptions on the regularity of the solution.
In the stochastic case, due to the limited time regularity (1.4), robust error quantities need to be used. Breit, Hofmanová and Loisel [18] use randomized time steps to construct an algorithm that achieves almost optimal convergence in time and optimal convergence in space, i.e., for all \(\alpha \in (0,1/2)\),
Here \(\mathbb {E}_{{\textbf{t}}}\) denotes the expectation with respect to the random time steps \({\textbf{t}}_m\) centered around deterministic grid-points \(t_m\). They use Hölder and Sobolev-Slobodeckij spaces to measure the time regularity of the solution. This excludes the limiting case \(\alpha = 1/2\).
The classical Euler–Maruyama scheme has been analyzed in [39]. The authors show, for general multi-valued monotone equations, convergence of the algorithm with convergence rate 1/4. In particular, the p-Laplace system can be treated. The authors infer convergence towards the spatially discrete solution \(v_h \in V_h\) of
with rate 1/4 measured in the error
The method has the advantage that no regularity on the abstract solution needs to be assumed. However, only a sub-optimal convergence rate can be obtained.
1.5 Main results
Based on deterministic time averages we propose two algorithms (3.14) and (3.15) that are essentially driven by the update rule, \(\mathbb {P}\)-a.s. for all \(m \ge 2\) and \(\xi _h \in V_h\)
The randomness enters through the averaged increments \(\Delta _m \mathbb {W}:= \langle {W}\rangle _m - \langle {W}\rangle _{m-1}\), where \(\langle {W}\rangle _m\) is the time averaged value of W on the interval \([t_{m-1},t_m]\).
Importantly, we manage to achieve optimal convergence in time with rate 1/2 without assuming any time Hölder regularity on the solution process. Instead, we measure time regularity in terms of Nikolskii spaces. Our main results, Theorem 19 and 25, verify under the condition
the optimal convergence
We want to stress that the regularity assumption (1.8) is even weaker than the provable regularity (1.4).
On the other hand if the solution process u has certain amount of time regularity, e.g. (1.4b), one can recover point evaluations, cf. Lemma 17,
For the implementation of (1.7) one needs to sample according to the distribution of the random variable \(\Delta _m \mathbb {W}\). We show in Corollary 33 that \(\Delta _m \mathbb {W}\) is a Gaussian random variable whose variance is slightly reduced compared to the classical increments \(\Delta _m W:= W(t_m) - W(t_{m-1})\). Ultimately, we provide a sampling algorithm (4.9) that samples not only the marginal distributions but the joint distribution of the random vector \((\Delta _m W, \Delta _m \mathbb {W})_{m=1}^M\).
1.6 Outline
In Sect. 2 we formulate the functional analytic setup, construct the multiplicative forcing term G and recall known regularity results. Section 3 introduces the discrete setup and contains the main results Theorem 19 and Theorem 25. Next, Sect. 4 clarifies the construction of the discrete random input vectors and provides a sampling algorithm. Lastly, Sect. 5 contains numerical experiments.
2 Mathematical setup
This section contains classical definitions and preliminary results. It is structured as follows: Sect. 2.1 introduces the function analytical framework. Section 2.2 presents the construction of the stochastic forcing. Section 2.3 is about the nonlinear operators S and V. The solution concept is fixed in Sect. 2.4. Lastly, Sect. 2.5 collects regularity results.
Let \({\mathcal {O}} \subset {\mathbb {R}}^n\) for \(n \ge 1\) be a bounded Lipschitz domain (further assumptions on \({\mathcal {O}}\) will be needed for the regularity of solutions). For some given \(T>0\) we denote by \(I := [0,T]\) the time interval and write \({\mathcal {O}}_T := I \times {\mathcal {O}}\) for the time space cylinder. Moreover let \(\left( \Omega ,{\mathcal {F}}, ({\mathcal {F}}_t)_{t\in I}, \mathbb {P} \right) \) denote a stochastic basis, i.e. a probability space with a complete and right continuous filtration \(({\mathcal {F}}_t)_{t\in I}\). We write \(f \lesssim g\) for two non-negative quantities f and g if f is bounded by g up to a multiplicative constant. Accordingly we define \( > rsim \) and \(\eqsim \). We denote by c a generic constant which can change its value from line to line.
2.1 Function spaces
As usual, \(L^q(\mathcal {O})\) denotes the Lebesgue space and \(W^{1,q}(\mathcal {O})\) the Sobolev space, where \(1\le q < \infty \). We denote by \(W^{1,q}_0(\mathcal {O})\) the Sobolev spaces with zero boundary values. It is the closure of \(C^\infty _0(\mathcal {O})\) (smooth functions with compact support) in the \(W^{1,q}(\mathcal {O})\)-norm. We denote by \(W^{-1,q'}(\mathcal {O})\) the dual of \(W^{1,q}_0(\mathcal {O})\). We do not distinguish in the notation between vector- and matrix-valued functions.
For a Banach space \(\left( X, \left\Vert \cdot \right\Vert _X \right) \) let \(L^q(I;X)\) be the Bochner space of Bochner-measurable functions \(u: I \rightarrow X\) satisfying \(t \mapsto \left\Vert u(t)\right\Vert _X \in L^q(I)\). Moreover, C(I; X) is the space of continuous functions with respect to the norm-topology. We also use \(C^{\alpha }(I;X)\) for the space of Hölder continuous functions. Given an Orlicz-function \(\varPhi : [0,\infty ] \rightarrow [0,\infty ]\), i.e. a convex function satisfying \( \lim _{t \rightarrow 0} \varPhi (t)/t = 0\) and \(\lim _{t \rightarrow \infty } \varPhi (t)/t = \infty \) we define the Luxemburg-norm
The Orlicz space \(L^\varPhi (I;X)\) is the space of all Bochner-measurable functions with finite Luxemburg-norm. For more details on Orlicz-spaces we refer to [32]. Given \(h \in I\) and \(u :I \rightarrow X\) we define the difference operator \(\tau _h: {\{{u: I \rightarrow X}\}} \rightarrow {\{{u: I\cap I - {\{{h}\}} \rightarrow X}\}} \) via \(\tau _h(u) (s) := u(s+h) - u(s)\). The Besov-Orlicz space \(B^\alpha _{\varPhi ,r}(I;X)\) with differentiability \(\alpha \in (0,1)\), integrability \(\varPhi \) and fine index \(r \in (1,\infty ]\) is defined as the space of Bochner-measurable functions with finite Besov-Orlicz norm \(\left\Vert \cdot \right\Vert _{B^\alpha _{\varPhi ,r}(I;X)}\), where
The case \(r = \infty \) is commonly called Nikolskii-Orlicz space. When \(\varPhi (t) = t^p\) we write \(B^\alpha _{p,r}(I;X)\) and call the space Besov space.
Similarly, given a Banach space \(\left( Y, \left\Vert \cdot \right\Vert _Y \right) \), we define \(L^q(\Omega ;Y)\) as the Bochner space of Bochner-measurable functions \(u: \Omega \rightarrow Y\) satisfying \(\omega \mapsto \left\Vert u(\omega )\right\Vert _Y \in L^q(\Omega )\). The space \(L^q_{\mathcal {F}}(\Omega \times I;X)\) denotes the subspace of X-valued \((\mathcal {F}_t)_{t\in I}\)-progressively measurable processes. Let \(\left( U,\left\Vert \cdot \right\Vert _U \right) \) be a separable Hilbert space. \(L_2(U;L^2(\mathcal {O}))\) denotes the space of Hilbert-Schmidt operators from U to \(L^2(\mathcal {O})\) with the norm \(\left\Vert z\right\Vert _{L_2(U;L^2_x)}^2:= \sum _{j\in \mathbb {N}} \left\Vert z(u_j)\right\Vert _{L^2(\mathcal {O})}^2 \) where \({\{{u_j}\}}_{j \in \mathbb {N}}\) is some orthonormal basis of U. We abbreviate the notation \(L^q_\omega L^q_t L^q_x := L^q(\Omega ;L^q(I;L^q(\mathcal {O}))) \) with obvious modifications for Sobolev, Besov and Hölder spaces. Additionally, we write \(L^{q-} := \bigcap _{r< q} L^r\).
2.2 Stochastic integrals
In order to construct the stochastic forcing term, we impose the following conditions:
Assumption 1
-
(a)
Let \(\left( U, \left\Vert \cdot \right\Vert _U \right) \) be a separable Hilbert space. We assume that W is an U-valued cylindrical Wiener process with respect to \((\mathcal {F}_t)_{t\in I}\). Formally W can be represented as
$$\begin{aligned} W = \sum _{j \in \mathbb {N}} u_j \beta ^j, \end{aligned}$$(2.1)where \({\{{u_j}\}}_{j \in \mathbb {N}}\) is an orthonormal basis of U and \({\{{\beta ^j}\}}_{j\in \mathbb {N}}\) are independent 1-dimensional standard Brownian motions.
-
(b)
Let \(v \in L^2_{\mathcal {F}}(\Omega \times I; L^2_x)\). We assume that \(G(v)(\cdot ): U \rightarrow L_{\mathcal {F}}^2(\Omega \times I; L^2_x)\) is given by
$$\begin{aligned} u \mapsto G(v)(u) := \sum _{j \in \mathbb {N}} g_j(\cdot ,v)( u_j, u)_U , \end{aligned}$$where \({\{{g_j}\}}_{j\in \mathbb {N}} \in C(\mathcal {O} \times {\mathbb {R}}^N; {\mathbb {R}}^N)\) with
-
(i)
(sublinear growth) for all \(x \in \mathcal {O}\) and \(\xi \in {\mathbb {R}}^N\) it holds
$$\begin{aligned} \sum _{j\in \mathbb {N}} \left|g_j(x,\xi )\right|^2 \le c_{\text {growth}}(1+\left|\xi \right|^2), \end{aligned}$$(2.2) -
(ii)
(Lipschitz continuity) for all \(x \in \mathcal {O}\) and \(\xi _1, \xi _2 \in {\mathbb {R}}^N\) it holds
$$\begin{aligned} \sum _{j\in \mathbb {N}} \left|g_j(x,\xi _1) -g_j(x,\xi _2) \right|^2 \le c_{\text {lip}} \left|\xi _1 - \xi _2\right|^2. \end{aligned}$$(2.3)
-
(i)
Now it is a classical result, see for example the book of Prévôt and Röckner [66], that we can construct a corresponding stochastic integral.
Proposition 2
Let Assumption 1 be true. Then the operator \(\mathcal {I}\) defined through
defines a bounded linear operator from \(L^2_{\mathcal {F}}(\Omega \times I;L_2(U; L^2_x))\) to \( L^2_\omega C_t L^2_x \). Moreover,
-
(a)
\(\mathcal {I}(G(v))\) is an \(L^2_x\)-valued martingale with respect to \(\left( \mathcal {F}_t\right) _{t \in I}\),
-
(b)
(Itô isometry) for all \(t \in I\) it holds
$$\begin{aligned} \mathbb {E} \left[ \left\Vert \mathcal {I}(G(v)) (t)\right\Vert _{L^2_x}^2 \right] =\mathbb {E} \left[ \int _0^t \left\Vert G(v)\right\Vert ^2_{L_2(U;L^2_x)} \,\textrm{d}s\right] . \end{aligned}$$(2.5)
2.3 Perturbed gradient flow
Let \(\kappa \ge 0\) and \(p \in (1,\infty )\). For \(\xi \in \mathbb {R}^{N \times n}\) we define
and
where \(\varphi (t):= \int _0^t \left( \kappa + s \right) ^{p-2} s \,\textrm{d}s\). The nonlinear functions S and V are closely related. In particular the following lemmata are of great importance. The proofs can be found in e.g. [31]. For more details we refer to [10, 27, 29].
Lemma 3
(V-coercivity) Let \(\xi _1, \xi _2 \in \mathbb {R}^{N \times n}\). Then it holds
Lemma 4
(generalized Young’s inequality) Let \(\xi _1, \xi _2, \xi _3 \in \mathbb {R}^{N \times n}\) and \(\delta >0\). Then there exists \(c_\delta \ge 1\) such that
Lemma 5
Let \(\xi _1, \xi _2, \xi _3 \in \mathbb {R}^{N \times n}\) and \(\delta >0\). Then there exists \(c_\delta \ge 1\) such that
Remark 6
A continuous, convex and strictly increasing function \(\varphi :[0,\infty ) \rightarrow [0,\infty )\) satisfying
is called an N-function. It is called uniformly convex if \(\varphi \in C^1[0,\infty ) \cap C^2(0,\infty )\) and \(\varphi '(t) \eqsim \varphi ''(t) t\) uniformly in \(t > 0\). Lemmas 3 and 4 are still valid if one replaces \(\varphi \) in (2.6) and (2.7) by any uniformly convex N-function. For more details we refer to e.g. [27].
Given some initial condition \(u_0 : \Omega \times \mathcal {O} \rightarrow {\mathbb {R}}^N\) and a stochastic force (G, W) in the sense of Assumption 1 we are interested in the system
with boundary and initial conditions given by
The system (2.11a) is a perturbed version of the gradient flow of the energy \(\mathcal {J}: W^{1,p}_{0,x} \rightarrow [0,\infty )\) given by
2.4 Weak and strong solutions
We fix the concept of solutions as follows.
Definition 7
Let \(u_0 \in L^2_\omega L^2_x\) be \(\mathcal {F}_0\)-measurable, \(p \in (1,\infty )\) and (G, W) be given by Assumption 1. An \((\mathcal {F}_t)\)-adapted process \(u \in L^2_x\) is called weak solution to (2.11) if
-
(a)
\(u \in L^2_\omega C_t L^2_x \cap L^p_\omega L^p_t W^{1,p}_{0,x}\),
-
(b)
for all \(t \in I\), \(\xi \in C^\infty _{0,x}\) and \(\mathbb {P}-a.s.\) it holds
$$\begin{aligned} \int _{\mathcal {O}} (u(t) - u_0) \cdot \xi \,\textrm{d}x+ \int _0^t \int _{\mathcal {O}} S(\nabla u) :\nabla \xi \,\textrm{d}x\,\textrm{d}s= \int _{\mathcal {O}} \int _0^t G(u)\,\textrm{d}W_s \cdot \xi \,\textrm{d}x. \end{aligned}$$(2.13)
The process u is called strong solution if it is a weak solution and additionally satisfies
-
(a)
\({\text {div}}S(\nabla u) \in L^2_\omega L^2_t L^2_x\) ,
-
(b)
for all \(t \in I\) and \(\mathbb {P}-a.s.\) it holds
$$\begin{aligned} u(t) - u_0 - \int _0^t {\text {div}}S(\nabla u) \,\textrm{d}s= \int _0^t G(u) \,\textrm{d}W_s \end{aligned}$$(2.14)as an equation in \(L^2_x\).
The next step is to answer the question about existence of weak or even strong solutions. Weak solutions can be constructed through a variational approach that uses the monotonicity of the nonlinear diffusion operator S as presented in [61, Example 4.1.9]. They rely on the Gelfand triple \(W^{1,p}_{0,x} \hookrightarrow L^2_x \hookrightarrow W^{-1,p'}_x\) for \(p \ge 2\). A quick observation shows that the Gelfand triple remains valid for \(p \ge 2n/(n+2)\). Below the threshold we need to modify the energy space. A good choice is the triple \(W^{1,p}_{0,x} \cap L^2_x \hookrightarrow L^2_x \hookrightarrow \big ( W^{1,p}_{0,x} \cap L^2_x \big )'\). Now, similar arguments as done in [61, Section 4] lead to the existence of a unique weak solution. Proving existence of a strong solution is more delicate and usually requires, not only assumptions on the growth of G, but also on the gradient of G, e.g. for all \(x \in \mathcal {O}\), \(\xi \in {\mathbb {R}}^N\)
In [43] a general approach for the construction of strong solutions to gradient flow like equations is presented. In particular, it includes the case of the p-Laplace equations.
Theorem 8
([43], Theorem 4.12) Assume \(p\ge 2\) and \(\mathcal {O}\) is a bounded convex domain. Let (G, W) be given by Assumption 1. Additionally, assume (2.15) and \(u_0 \in L^{p+\varepsilon }_\omega W^{1,p}_x\cap L_\omega ^{2+\varepsilon } L^2_x\) be \(\mathcal {F}_0\)-measurable for some \(\varepsilon > 0\). Then there exists a unique strong solution u to (2.11). Moreover,
2.5 Regularity of strong solutions
A key ingredient in the error analysis of numerical algorithms are the improved regularity properties of strong solutions in comparison to those of weak solutions. Concerning time regularity, we prove in [76] that the strong solution enjoys 1/2-time differentiability in an exponential Besov space and even the nonlinear greadient \(V(\nabla u)\) obeys 1/2-time differentiability in a Nikolskii space. The proof uses an assumption on the boundary condition of the noise coefficient G, i.e. for all \(x \in \partial \mathcal {O}\),
Theorem 9
([76] Theorem 3.8 & Theorem 3.11) Let the assumptions of Theorem 8 be satisfied. Additionally, assume (2.17). Let u be the strong solution to (2.11). Then
where \(\varPhi _2(t) = e^{t^2}-1\).
Spatial regularity is closely connected to the existence of strong solutions. Local regularity has been obtained in [15].
Theorem 10
([15], Theorem 4) Let \(u_0 \in L^2_\omega W^{1,2}_x\) be \(\mathcal {F}_0\)-measurable. Let Assumption 1 be satisfied. Additionally, assume (2.15) and denote by u the strong solution of (2.11). Then,
On sufficiently regular domains, it is possible to relate the divergence of the nonlinear operator S to the full gradient as presented in [2, 21], i.e.
Precisely, given some bounded Lipschitz domain \(\mathcal {O}\) such that \(\partial \mathcal {O} \in W^{2,1}\), i.e. \(\mathcal {O}\) is locally the subgraph of a Lipschitz continuous function of \(n-1\) variables, which is also twice weakly differentiable. Denote by \(\mathcal {B}\) the second fundamental form on \(\partial \mathcal {O}\), by \(\left|\mathcal {B}\right|\) its norm and define
where \(B_r(x)\) denotes the ball of radius r around x, \(\text {cap}_{B_1(x)}(E)\) is the capacity of the set E relative to the ball \(B_1(x)\) and \(\mathcal {H}^{n-1}\) is the \(n-1\) dimensional Hausdorff measure. The following result follows from [21, Theorem 2.1].
Lemma 11
Assume that \(\mathcal {O}\) is either
-
(a)
bounded and convex,
-
(b)
or bounded, Lipschitz and \(\partial \mathcal {O} \in W^{2,1}\) with \(\lim _{r\rightarrow 0} \mathcal {K}_\mathcal {O}(r) \le c\).
Let \(v \in L^p_\omega L^p_t W^{1,p}_{0,x}\) with \({\text {div}}S(\nabla v) \in L^2_\omega L^2_t L^2_x\). Then \(\nabla S(\nabla v) \in L^2_\omega L^2_t L^2_x\) and
In the non-degenerate setting, \(\kappa > 0\), this allows to deduce global spatial regularity.
Corollary 12
([76], Corollary 2.14) Let the assumptions of Theorem 8 be satisfied and \(\kappa > 0\). Let u be the strong solution to (2.11). Then,
3 Numerical scheme for the averaged system
In this section we will first present the discrete structures. Afterwards we construct two algorithms that approximate the solution to (2.11). Finally, we prove convergence of the approximation towards the analytic solution with optimal rates.
3.1 Space discretization
Let \(\mathcal {O} \subset {\mathbb {R}}^n\) be a bounded polyhedral domain. By \(\mathcal {T}_h\) denote a regular partition (triangulation) of \(\mathcal {O}\) (no hanging nodes), which consists of closed n-simplices called elements. For each element (n-simplex) \(T\in \mathcal {T}_h\) we denote by \(h_T\) the diameter of T, and by \(\rho _T\) the supremum of the diameters of inscribed balls.
We assume that \(\mathcal {T}_h\) is shape regular, that is there exists a constant \(\gamma \) (the shape regularity constant) such that
We define the maximal mesh-size by
We assume further that our triangulation is quasi-uniform, i.e.
For \(\ell \in \mathbb {N}_0\) we denote by \({\mathscr {P}}_\ell (\mathcal {O})\) the polynomials on \(\mathcal {O}\) of degree less than or equal to \(\ell \).
For fixed \(r \in \mathbb {N}\) we define the vector valued finite element space \(V_h\) as
Moreover, let \(\Pi _2: L^2_x \rightarrow V_h\) be the \(L^2_x\)-orthogonal projection defined by
or equivalently
We will need some classical results on the stability properties of the \(L^2_x\)-projection for finite elements.
Lemma 13
Let \(r \ge 1\), \(V_h\) be given by (3.3) and \(\mathcal {T}_h\) be quasi-uniform. Then
Due to the nonlinear structure of the p-Laplace system we also need an adapted stability result.
Proposition 14
([16], Theorem 7) Let \(r \ge 1\), \(V_h\) be given by (3.3) and \(\mathcal {T}_h\) be quasi-uniform. Then
3.2 Time discretization
Let \(\{0=t_0<\cdots <t_M=T\}\) be a uniform partition of [0, T] with mesh size \(\tau =T/M\). For \(m \ge 1\) define \(I_m := [t_{m-1},t_m]\). By \(\left|I_m\right|\) we denote the Lebesgue measure of \(I_m\). By 
we denote the mean value integral over the set \(I_m\). We also abbreviate 
for the mean value.
Let us discuss the stability properties of piecewise constant approximations generated by the mean values in terms of Nikolskii spaces.
Lemma 15
Let \(\alpha \in (0,1)\) and \(r \in [1,\infty )\). Additionally, assume \(u \in B^{\alpha }_{r,\infty } L^2_x\). Then
Proof
Due to Jensen’s inequality
A change of variables \(t = s+h\) and Fubini’s theorem yield
The proof is finished. \(\square \)
Remark 16
Lemma 15 is also valid for \(r=\infty \). Here we need to substitute the left hand side in (3.4) by
Note that \(B^{\alpha }_{\infty ,\infty }(I) = C^{\alpha }(I)\) is the space of \(\alpha \)-Hölder continuous functions.
In our application the process u only belongs to \(L^2_\omega B^{1/2}_{\varPhi _2,\infty } L^2_x \backslash L^2_\omega C^{1/2}_t L^2_x \). Therefore we need to adjust the stability result in terms of exponentially integrable Nikolskii spaces.
Lemma 17
Let \(\alpha \in (0,1)\) and \(u \in B^{\alpha }_{\varPhi _2,\infty } L^2_x\) where \(\varPhi _2(t) = e^{t^2}-1\). Then
Proof
Define \(I_m^k:= [t_{m} - 2^{-k}\tau , t_m]\). Due to the embedding \(B^{\alpha }_{\varPhi _2,\infty } \hookrightarrow C_t\) for all \(\alpha > 0\) we find u to be continuous as an \(L^2_x\)-valued process. Therefore, it holds
Observe, \(\left|I_m^{k}\right| = 2^{-k} \tau \) and \(I_m^{k-1} = [t_m - 2^{-(k-1)}\tau ,t_m - 2^{-k}\tau ] \cup [t_m - 2^{-k}\tau ,t_m]\). A shift in the integral results in
Jensen’s inequality implies
Choose \(\lambda _k \) by
In particular, since \(I_m^k \subset I \cap I- {\{{2^{-k}\tau }\}}\),
Thus,
Since \(u \in B^{\alpha }_{\varPhi _2,\infty } L^2_x \) it holds
It follows
Note that
Therefore,
The assertion follows by using the estimate (3.7) in (3.6). \(\square \)
3.3 The averaged algorithm
Let u be the strong solution of (2.11), i.e. for all \(t \in I\) and \(\mathbb {P}\)-a.s.
Take the average over \(I_1\) in (3.8) to get
To obtain the general evolution of the time averaged values of the solution we subtract (3.8) for t and \(t-\tau \) and take the average over \(I_m\)
We denote by \((\cdot ,\cdot )\) the \(L^2_x\)-inner product. The corresponding weak formulation reads for all \(\xi \in L^2_x \cap W^{1,p}_{0,x}\)
Due to the (stochastic) Fubini’s Theorem we can equivalently write
where the weights are given by
The above considerations motivate the construction of the following numerical scheme:
We initialize the algorithm as
In order to accurately reflect the special character of the first step (3.9), we define \(v_1\) via
Moreover the evolution for \(m \in {\{{2,\ldots ,M}\}}\) is determinded via
for all \(\xi _h \in V_h\) and \(\mathbb {P}\)-a.s.
The need of a special step size in the first step (3.14b) might be undesirable. It can be overcome by performing a full initial step. We propose the following second algorithm.
Initialize
Full initial step
Time stepping, for \(m \ge 2\),
The additional difficulties in the error analysis only occur in the estimate of the initial error.
The next theorem ensures that the numerical schemes are well defined. The arguments are rather standard and we only refer to [40] Section 3 for a detailed analysis of a more general system.
Theorem 18
Let \(u_0 \in L^2_\omega L^2_x\) be \(\mathcal {F}_0\)-measurable and Assumption 1 be satisfied. Then for all \(M \in \mathbb {N}\) there exists a unique \({\textbf{v}}= (v_0,\ldots ,v_M) \in (V_h)^{M+1}\) such that \(v_m\) is \(\mathcal {F}_{t_m}\)-measurable for all \(m \in {\{{0,\ldots ,M}\}}\) and \({\textbf{v}}\) solves (3.14) respectively (3.15).
3.4 Error analysis
We are ready to state and prove the main result.
Theorem 19
(Convergence of Algorithm (3.14)) Let \(u_0 \in L^2_\omega L^2_x\) be \(\mathcal {F}_0\)-measurable and Assumption 1 be satisfied. Additionally, assume \(\mathcal {T}_h\) to be quasi-uniform. Moreover, let u be the weak solution to (2.11) and \({\textbf{v}}\) be the numerical solution of (3.14).
-
(a)
Assume
$$\begin{aligned} u&\in L^2_\omega B^{1/2}_{2,\infty } L^2_x \cap L^2_\omega L^\infty _t W^{1,2}_x, \end{aligned}$$(3.16a)$$\begin{aligned} V(\nabla u)&\in L^2_\omega B^{1/2}_{2,\infty } L^2_x \cap L^2_\omega L^2_t W^{1,2}_x. \end{aligned}$$(3.16b)Then it holds
$$\begin{aligned} \begin{aligned}&\mathbb {E} \left[ \max _{m=1,\ldots ,M} \left\Vert \langle {u}\rangle _m - v_m\right\Vert _{L^2_x}^2 + \sum _{m=1}^M\int _{I_m} \left\Vert V(\nabla u_\nu ) - V(\nabla v_m)\right\Vert _{L^2_x}^2 \,\textrm{d}\nu \right] \\&\quad \lesssim \tau + h^2. \end{aligned} \end{aligned}$$(3.17) -
(b)
If additionally
$$\begin{aligned} u \in L^2_\omega B^{1/2}_{\varPhi _2,\infty } L^2_x, \end{aligned}$$(3.18)then
$$\begin{aligned} \begin{aligned} \mathbb {E} \left[ \max _{m=1,\ldots ,M} \left\Vert u(t_m)- v_m\right\Vert _{L^2_x}^2 \right] \lesssim \tau \ln (1+\tau ^{-1}) + h^2. \end{aligned} \end{aligned}$$(3.19)
Remark 20
The regularity assumption (3.18) is optimal in the sense, that it is the limiting space of the time regularity of the Wiener process W. Hytönen and Veraar prove in [50] that a Banach space valued Brownian motion has full 1/2-differentiability only on the Nikolskii-scale. More precisely, they show for all \(q,r \in [1,\infty )\) and \(\mathbb {P}\)-a.s.
The logarithmic term in (3.19) has already been used in the context of rough stochastic differential equations [23, 55]. It quantifies the distance of \(L^{\varPhi _2}\) and \(L^\infty \).
Theorem 19 can be generalized to cover fractional time and space regularity assumptions as done for the deterministic p-Laplace system in [16].
Proof
Part (a): Fix \(m \in {\{{2,\ldots ,M}\}}\), subtract (3.11) from (3.14c) and choose \(\xi _h = \Pi _2 e_m := \Pi _2 \langle {u}\rangle _m - v_m \)
Now use the symmetry of the \(L^2_x\)-projection and the algebraic identity \(2a(a-b) = a^2 - b^2 + (a-b)^2\) to get
The second term, due to the V-coercivity (2.8),
Summation over \(m \in {\{{2,\ldots ,m^*}\}}\) with \(m^* \le M\) results in
Take the maximum over \(m^* \in {\{{2,\ldots ,M}\}}\) and expectation
Step 1: We start by estimating the initial error. Similarly as before, we subtract the weak formulation of (3.9) from (3.14b) and choose \(\xi _h = \Pi _2 e_1\)
where \(e_0 := u_0 - v_0\). By (3.14a) we have \(\Pi _2 e_0 = 0\). Therefore,
Due to Hölder’s and Young’s inequalities and Itô isometry
Absorb the second term to the left hand side. For the first term we use \(a_0 \le 1\) and the Lipschitz assumption (2.3). Now, since the operator norm of the \(L^2_x\)-projection is bounded by one,
Thanks to the generalized Young’s inequality, cf. Lemma 4,
Absorb the first term to the left hand side in (3.21). The second is split up into a space and a time error. Using the nonlinear stability of the \(L^2_x\)-projection, cf. Proposition 14,
Overall, we arrive at the estimate from (3.21)
Step 2: The stochastic part in (3.20) needs some refined analysis
The first, due to Hölder’s and Young’s inequalities and an index shift,
The second term can be absorbed to the left hand side in (3.20). The third term is nothing but the initial error \(J_1\). For the first term we invoke Itô isometry, the Lipschitz condition (2.3) and \(a_{m-1} \le 1\)
Decomposition in time and space error, applying Lemmas 13 and 15 and the estimate (3.22)
Next, we analyze the second term \(J_{2,b}\) in the upper bound for \(J_2\). Define the discrete real-valued stochastic process
It is convenient to use stochastic Fubini’s theorem to rewrite
In the following we abbreviate \(\mathcal {F}_m:= \mathcal {F}_{t_m}\). Note, (3.25) does not define a martingale with respect to \(\mathcal {F}_{t_m}\). In fact, the discrepancy of not being a martingale can be quantified. The general strategy is to split up the sum into a martingale and an error term. The error term is called compensator. To determine how the compensator looks, we first compute the conditional expectations of \(K_M\) with respect to \(\mathcal {F}_{m^*}\), i.e.,
Due to the tower property of conditional expectation, the measurability of \(e_m\) with respect to \(\mathcal {F}_m\) together with the martingale property of the stochastic integral
Again using the \(\mathcal {F}_m\)-measurability of \(e_m\), we conclude that \(K_{m^*}\) is \(\mathcal {F}_{m^*}\)-measurable. Thus,
It remains to compute the conditional expectation in \(\mathcal {M}_2\). Since \(t_{m^*}\) is an interior point of \(I_{m^*+1} \cup I_{m^*} \) we split up the stochastic integral into a part that only sees values above the threshold \(t_{m^*}\) and into a lower part that only sees values below the threshold \(t_{m^*}\), i.e.,
The first vanishes due to the martingale property of the stochastic integral, while the second is measurable with respect to \(\mathcal {F}_{m^*}\). Overall,
\(\mathcal {M}_2\) is called compensator and quantifies the error of not being a martingale, i.e.,
Furthermore, increments of the discrete stochastic process K satisfy
(3.26) together with (3.27) allow to identify increments of the conditional expectations,
Observe \(\mathbb {E}\left[ K_M \big | \mathcal {F}_1 \right] = 0\), since \(\Pi _2 e_0 = 0\). The Burkholder-Davis-Gundy’s inequality implies
Now Young’s inequality, Itô isometry and the Lipschitz condition (2.3) imply
The second term is estimated as in (3.24)
Similarly, due to (3.26), Hölder’s inequality, Young’s inequality and \(\ell ^1 \hookrightarrow \ell ^\infty \)
Together we can estimate
This concludes the bound for \(J_2\) in (3.20)
Step 3: In this step we estimate the nonlinear gradient. Jensen’s inequality, the generalized Young’s inequality (2.9) and the nonlinear stability result Proposition 14 imply
Step 4: We aim at applying a Gronwall type argument. Collecting all estimates we get
Choosing \(\varepsilon \) sufficiently small and applying Gronwall’s Lemma ensures
This implies
Step 5: We artificially introduce the desired error quantities and use the stability of mean-value time projections and the \(L^2\)-space projection. Let us denote 
. A slight modification of Lemma 15 implies
Additionally, Jensen’s inequality ensures
Therefore, invoking (3.29),
The assertion (3.17) follows by an application of Lemma 13
Part (b): Now, let us assume \(u \in L^2_\omega B^{1/2}_{\varPhi _2,\infty } L^2_x\). We apply Lemma 17 to bound
This verifies (3.19) and the proof is finished. \(\square \)
One can also use time averages on the nonlinear gradient term to measure the error of the approximation.
Corollary 21
Let the assumptions of Theorem 19(a) be satisfied. Then
and
Proof
The estimate (3.30) immediately follows by an application of Jensen’s inequality and the bound (3.17),
In order to prove the second estimate (3.31) we use Lemmas 3 and 4
Absorbing the second term to the left hand side and applying the bound (3.17) verifies the assertion. \(\square \)
Remark 22
Although both (3.30) and (3.31) enjoy the same convergence rates, it is not clear whether one dominates the other. In the linear case, \(p=2\), the terms coincide.
The averaged error quantities (3.30) and (3.31) are equivalent up to oscillation to the error quantity (3.17).
Lemma 23
Let \(u \in L^p_t W^{1,p}_x\) and \(A \in L^p_x\). Then
and
Proof
The Eq. (3.32) follows by using
The estimate (3.33) is obtained trivially,
\(\square \)
Remark 24
In [28, Lemma 6.2] the authors prove the equivalence (although it is done purely in space but can be extended to time) of
If \(V(\nabla u) \in B^{1/2}_{2,\infty } L^2_x \), then Lemma 15 implies
Theorem 25
(Convergence of Algorithm 3.15) Let the assumptions of Theorem 19 be satisfied. Denote by \({\textbf{w}}\in (V_h)^{M+1}\) the solution to (3.15) and by u the weak solution to (2.11). Then
and
Proof
The proof proceeds similarly to the proof of Theorem 19. We will only prove the bound for the initial error. Instead of comparing \(w_1\) to \(\langle {u}\rangle _1\), we rather choose 
. The equation for the latter one reads
Subtracting (3.15b) from (3.36) and choosing \(\xi _h = \Pi _2 \langle {u}\rangle _{I_1 \cup I_2} - w_1\) results in
Due to Hölder’s and Young’s inequalities, Itô isometry and the growth assumption (2.2)
The boundedness of u as an \(L^2_x\) valued-process implies
The forth term is estimated similarly,
Using the \(L^2\)-stability of the \(L^2\)-projection we obtain
It remains to check the second term. Here we use the same arguments as in step 3 of the proof of Theorem 19 to conclude
Overall, we have established
Lastly, using Lemmas 13 and 15
The bound for the initial error is complete. \(\square \)
Remark 26
We want to remark, although the first step (3.15b) is announced to be a full step, it is not a full first step. The deterministic drift is scaled by a full step \(\tau \), but the stochastic term is only scaled by \(\langle {W}\rangle _1\). In Corollary 33 we find that \(\langle {W}\rangle _1\) has only half the variance of a full stochastic step.
To overcome this problem we can introduce a stochastic dummy variable \(\mathcal {W}_0 \sim \mathcal {N}(0,\tau /3)\)Footnote 1 that artificially makes up for the missing randomness. Then (3.15b) needs to be substituted by
The error analysis of Theorem 25 can be extended to this algorithm.
4 Discrete stochastic processes
In this section we investigate the law of averaged Wiener processes and propose an implementable sampling algorithm.
4.1 Wiener process
We introduce the concept of Hilbert space valued Gaussian processes.
Definition 27
A stochastic process Y is called U-valued Gaussian process with mean operator \(m: I \times U \rightarrow {\mathbb {R}}\) and variance operator \(\Sigma : I \times U \times U \rightarrow {\mathbb {R}}\), if for all \(t \in I\) and \(u \in U\) it holds
In short, we write \(Y_t \sim \mathcal {N}\left( m_t, \Sigma _t \right) \).
It follows that the stochastic forcing W defined by (2.1) is an U-valued Gaussian process with mean operator \(m_t(u) = 0\) and variance operator \(\Sigma _t(u,v) := t (u,v)_U\) for all \(t\in I\) and \(u,v \in U\). Moreover, the series (2.1) converges in the weak topology of U, due to the hypercontractivity of normally distributed random variables it holds for any \(q >0\), \(u \in U\) and \(t,s \in I\)
In fact, as soon as the index set is infinite we lose the norm convergence, since
4.2 Averaged Wiener process
In this section we compute the distributions of the random variables \((\langle {W}\rangle _m)_{m=1}^M\) and the joint distribution of the averaged increments.
A key tool in the derivation of the distribution of the averaged Wiener process is the decomposition of the process W adjusted to the equidistant time partition \({\{{I_m}\}}_{m=1}^M\). We decompose \(W|_{I_m}\) into a Brownian bridge \(\mathcal {B}_m\) and its nodal values \(W(t_{m-1})\) and \(W(t_{m})\), i.e., for \(t \in I_m\)
where
Brownian bridges have nice independency properties. They do not look into the past nor future. The following result can be found in [63, Section 1.2].
Proposition 28
Let \((\mathcal {B}_m)_{m=1}^M\) be given by (4.2). Then for all \(m \in {\{{1,\ldots ,M}\}}\)
i.e. all finite dimensional distributions of the generators of each sigma algebra are independent of each other.
Corollary 29
Let \((\mathcal {B}_m)_{m=1}^M\) be given by (4.2). Then \(\mathcal {B}_1, \ldots , \mathcal {B}_M, \Delta _1 W, \ldots , \Delta _M W\) are independent.
Next, we take the time average over the interval \(I_m\) in (4.1) and obtain
Now, it is our choice whether we want to compute the distribution of \(\langle {W}\rangle _m\) or \(\langle {\mathcal {B}_m}\rangle _m\). The formula (4.4) provides an easy way to compute the remaining one. We choose to compute the distribution of \(\langle {W}\rangle _m\). An application of Itô’s formula for \(f(s,W_s) = \frac{s-t_m}{\tau } W_s\) implies \(\mathbb {P}\)-a.s.
A stochastic integral that is driven by a Wiener process and a deterministic integrand stays Gaussian.
Lemma 30
([4] Prop. 7.1) Let \(f \in L^2(I)\). Then
is an U-valued Gaussian process with zero mean and variance
Corollary 31
\(\langle {W}\rangle _m\) is an U-valued Gaussian random variable with zero mean and variance \(\Sigma (u,v) = \frac{2t_{m-1} + t_m}{3}(u,v)_U \).
Proof
Let us define
Note, that \(W_a\) and \(W_b\) are independent. Moreover, \(W_a\) is Gaussian with variance \(\Sigma _a(u,v) = t_{m-1} \left( u,v \right) _U\) and due to Lemma 30\(W_b\) is also Gaussian. Therefore \(\langle {W}\rangle _m\) is Gaussian and it suffices to compute the mean and the variance operators.
Let \(u,v \in U\). Then \(\mathbb {E} \left[ \left( \langle {W}\rangle _m, u \right) _{U} \right] = 0\) and
The assertion is proved. \(\square \)
Corollary 32
\(\langle {\mathcal {B}_m}\rangle _m\) is an U-valued Gaussian random variable with zero mean and variance \(\Sigma (u,v) = \frac{\tau }{12}(u,v)_U\).
Proof
The distribution of a random variable is uniquely determined by its characteristic function. Corollary 31 implies
Classically, we find
The characteristic function of \(\varphi _{\langle {W}\rangle _m}\) factors due to the independence of the decomposition (4.4), i.e.,
Rearranging implies
Overall, \(\langle {\mathcal {B}_m}\rangle _m\) has the characteristic function of a Gaussian random variable with zero mean and variance \(\Sigma (u,v) = \frac{\tau }{12}(u,v)_U\). \(\square \)
At this point it is a simple task to find the distribution of the increments of the averaged Wiener process. Let us subtract (4.4) for m and \(m-1\)
where we define \(\langle {W}\rangle _0 := \Delta _0 W:= \langle {\mathcal {B}_0}\rangle _0:= 0\).
Corollary 33
\(\Delta _m \mathbb {W}\) is an U-valued Gaussian random variable with zero mean and variance \(\Sigma (u,v) = \left( \frac{2\tau }{3} \chi _{{\{{m \ge 2}\}}} + \frac{\tau }{3} \chi _{{\{{m=1}\}}} \right) \left( u, v \right) _U\).
Proof
The right hand side of (4.6) is a sum of independent, centered Gaussian random variables. Thus the left hand side is centered Gaussian. Now, it suffices to compute the variance operator.
Note, in the case \(m=1\) we have \(\Delta _1 \mathbb {W} = \langle {W}\rangle _1\) and the result follows by Corollary 31.
Let \(m \ge 2\) and \(u,v \in U\). Due to the independence,
\(\square \)
So far we have identified how each averaged increment \(\Delta _m \mathbb {W}\) is distributed. However, we also need to know what the joint distribution is, i.e., the distribution of a random vector.
Lemma 34
The random vector \((\Delta _m \mathbb {W})_{m=1}^M\) is an \(U^M\)-valued centered Gaussian random variable with variance operator \(\Sigma :U^M \times U^M \rightarrow {\mathbb {R}}\) given by
for \({\textbf{u}}, {\textbf{v}}\in U^M\) and
Proof
The Eq. (4.6) implies, that the random vector \((\Delta _m \mathbb {W})_{m=1}^M\) can be constructed via a linear transformation of independent Gaussian random vectors \((\Delta _m W)_{m=1}^M\) and \((\langle {\mathcal {B}_m}\rangle _m)_{m=1}^M\), i.e.,
where \(\mathbb {K}_1, \mathbb {K}_2 \in {\mathbb {R}}^{M \times M}\) are given by
Therefore, \((\Delta _m \mathbb {W})_{m=1}^M\) is itself a centered Gaussian vector. It remains to compute the covariance matrix. Let \(u,v \in U\) and \(m, l \in {\{{1, \ldots , M}\}}\). If \(m = l\) Corollary 33 implies
If \(\left|m-l\right|> 1\), then Eq. (4.6) and the independence imply
It remains to consider the case \(\left|m-l\right|=1\). Without loss of generality \(l = m+1\). Now, using (4.6), the independence and Corollary 32
The proof is finished. \(\square \)
4.3 Sampling algorithm
On the computer we are forced to approximate the continuous measure induced by the Wiener process W by an empirical measure. Additionally, if we want to compare different numerical schemes we need to specify how to sample the random input needed for the involved algorithms jointly. More specifically, if we want to compare our algorithms (3.14) and (3.15) and the classical Euler–Maruyama discretization (5.1), we need to sample according to the law of the random vector
Based on the decomposition (4.6) we propose the following sampling algorithm. The time discretization is achieved by the parameter \(M \in \mathbb {N}\) and the series truncation in (2.1) is done by the parameter \(J \in \mathbb {N}\).
Given \(M,J \in \mathbb {N}\).
-
(a)
(Sampling) Compute i.i.d. random variables \(\zeta _{m}^j, \eta _{m}^j \sim \mathcal {N}(0,1)\) for \(m \in {\{{1,\ldots ,M}\}}\) and \(j \in {\{{1,\ldots ,J}\}}\).
-
(b)
(Lift to Hilbert space U) For \(m \in {\{{1,\ldots ,M}\}}\) define the random variables
$$\begin{aligned} Z_m := \sqrt{\tau } \sum _{j=1}^J u^j \zeta _{m}^j, \quad \tilde{\mathbb {Z}}_m := \sqrt{ \frac{\tau }{12} }\sum _{j=1}^J u^j \eta _{m}^j, \end{aligned}$$(4.9a)where \({\{{u^j}\}}_{j\in \mathbb {N}}\) is an orthonormal system of U.
-
(c)
(Adjusting correlation) For \(m \in {\{{1,\ldots ,M}\}}\) define the random variables
$$\begin{aligned} \mathbb {Z}_m := \frac{Z_m + Z_{m-1}}{2} + \tilde{\mathbb {Z}}_m - \tilde{\mathbb {Z}}_{m-1}, \end{aligned}$$(4.9b)where \(Z_0 := \tilde{\mathbb {Z}}_0 := 0\).
Let \(\Pi _J\) be the U-orthogonal projection onto \(U_J := \text {span}(u_1,\ldots , u_J)\). The following proposition guarantees that the sampling algorithm (4.9) approximates the desired random variables.
Proposition 35
Let \({\textbf{Z}}:= \left( Z_1, \ldots , Z_M,\mathbb {Z}_1,\ldots ,\mathbb {Z}_M \right) \in U^{2M}\) be generated by (4.9a) and (4.9b). Then,
where \(\Pi _J\) is the U-orthogonal projection onto \(\text {span}(u^1,\ldots , u^J)\).
Proof
First, we need to observe that \(Z_m \sim \Pi _J \Delta _m W\) and \({\tilde{Z}}_m \sim \Pi _J \langle {\mathcal {B}_m }\rangle _m\). Then, the statement follows similarly to the proof of Lemma 34. \(\square \)
Remark 36
Proposition 35 ensures that we can compare our algorithm to the classical Euler–Maruyama discretization (5.1) on an equidistant time grid. It can be adjusted to also match non-equidistant grids.
5 Simulations
In this section we perform numerical simulations to test our algorithms (3.14) and (3.15). We denote the solution of (3.14) as \({\textbf{v}}^{\text {Half}}\) and the solution of (3.15) as \({\textbf{v}}^{\text {Full}}\). Additionally, we compare our algorithms to the classical Euler–Maruyama discretization of (2.11). That reads, find \({\textbf{v}}\in (V_h)^{M+1}\) such that for all \(\xi _h \in V_h\), \(m \ge 1\) and \(\mathbb {P}\)-a.s.
The solution to (5.1) is called \({\textbf{v}}^{\text {EM}}\).
A plain convergence result for the Euler–Maruyama scheme without any rate has been obtained in [65]. A more sophisticated analysis has been done in [39]. However, they only obtain convergence for the \(L^2_\omega L^\infty _t L^2_x\)-error with rate 1/4.
Originally the Euler–Maruyama scheme has been introduced for stochastic ordinary differential equations. In this context much more is known, see e.g. the book of Kloeden and Platen [52, Section 9.5]. However, when dropping the Lipschitz assumption on the coefficients, divergence with positive probability has been obtained in [49].
We are particularly interested in the experimental study of the following questions:
-
(a)
Do the algorithms (3.14) and (3.15) approximate mean-values or point-values?
-
(b)
How does the Euler–Maruyama scheme (5.1) compares to (3.14) and (3.15) in terms of time and space convergence?
-
(c)
How do the different error quantities (3.17), (3.30) and (3.31) for the gradient relate to each other?
-
(d)
How sensitive are the algorithms with respect to the parameter p?
All simulations are done with the help of the open source tool for solving partial differential equations FEniCS [62].
5.1 An explicit solution
In the linear case, \(p=2\), with a linear right hand side
for some \(\lambda \in {\mathbb {R}}\), it is possible to find an explicit solution. If we start the evolution defined by (2.11) in an eigenfunction of the Laplace operator, the dynamics become simpler. Let \(\mathcal {O} = (0,1)^2\), \(T = 1\) and \(u_0(x) = \sin (\pi x_1) \sin ( \pi x_2)\). Note that \(u_0\) is an eigenfunction of the 2-Laplacian with homogeneous Dirichlet data and corresponding eigenvalue \(\mu = 2 \pi ^2\). The unique solution to (2.11) is given by
Similarly, we can give a closed expression for the solution \(u_h\) to the space-discrete equation
This is equivalent to the system of linear stochastic differential equations
where \(M_{i,j} = \left( \xi _h^j, \xi _h^i \right) \) and \(S_{i,j} = \left( \nabla \xi _h^j, \nabla \xi _h^i \right) \) is the mass respectively the stiffness matrix and \({\{{\xi _h^i}\}}\) form a basis of \(V_{h}\). The eigenpairs of the discretized Laplacian are related to the linear system
Let \((\mu _h,{\textbf{u}}_h) \in (0,\infty ) \times {\mathbb {R}}^{\left|V_{h}\right|}\) be a solution to (5.6), then the solution to (5.4) started in \(u_h(0) = {\textbf{u}}_h \cdot {\varvec{\xi }}_h\) is given by
The main advantage of having an analytic solution of the space discrete equation is, that it rules out any space discretization errors.
To accurately compare continuous processes and discrete vectors, one needs to either lift the vector to a process or project the process to a vector. We do the latter approach and evaluate the continuous process \(u_h\) on the equidistant partition of I. This leads to the definition
Since the algorithms (3.14) and (3.15) approximate the mean value of the analytic solution, we define
Although we know the exact solution, the time averages of the exact solution are non-treatable without knowledge of the full trajectory of the Brownian motion \(\beta ^1\). Numerically, we substitute \(\langle {u_h}\rangle _{I_m}\) by a Riemann sum approximation, i.e. we fix an equidistant partition \({\{{[t_{m,k-1},t_{m,k}]}\}}_{k=1}^r\) of \(I_m\) with resolution \(r \in \mathbb {N}\) and define
The approximation quality is measured in the error quantities
where \({\textbf{u}}\in {\{{ {\textbf{u}}^{\text {Point}}, {\textbf{u}}^{\text {Aver}}}\}}\) and \({\textbf{v}}\in {\{{{\textbf{v}}^{\text {EM}},{\textbf{v}}^{\text {Full}},{\textbf{v}}^{\text {Half}} }\}} \).
In Fig. we plot the time convergence of the error quantities (5.11a) and (5.11b). We approximate the expectation by the Monte-Carlo method with 20 samples. Additionally, we let \(\lambda = 1\) and \(V_h\) to be the space of piecewise linear, continuous elements on a uniform mesh with \(\left|V_h\right| = 121\). The average values (5.10) are approximated by \(r = 10\).
The numerical results support that \({\textbf{v}}^{\text {EM}}\) approximates the solution on the grid points, while \({\textbf{v}}^{\text {Full}}\) and \({\textbf{v}}^{\text {Half}}\) approximate the average values of the solution. The gap is due to the difference
Initially, we observe a preasymptotic effect. It stabilizes at the time scale \(\tau \approx 10^{-3}\). Afterwards the predicted convergence speed of order 1 is achieved.
Time convergence of \({\textbf{v}}^{\text {EM}}\) (red), \({\textbf{v}}^{\text {Half}}\) (blue) and \({\textbf{v}}^{\text {Full}}\) (green) towards \({\textbf{u}}^{\text {Point}}\) (dashed) respectively \({\textbf{u}}^{\text {Aver}}\) (solid) measured in the error terms \(\mathcal {E}\) (left) and \(\mathcal {V}\) (right) (color figure online)
5.2 Beyond known solutions
In general, a major obstacle is the absence of an analytic solution to the Eq. (2.11). In particular, the distance of the numerical solution and the analytic solution as presented in (3.17), (3.19) respectively (3.34) and (3.35) are non-computable.
To overcome this difficulty we measure the error of a fine reference approximation \({\textbf{v}}_f\) and a coarse approximation \({\textbf{v}}_c\). Both, \({\textbf{v}}_f\) and \({\textbf{v}}_c\), are generated via the same algorithm (either (3.14), (3.15) or (5.1)) on a fine respectively coarse scale. Let \(h_c \ge h_f >0\) be coarse respectively fine space mesh sizes. Similarly, let \(M_c, M_f \in \mathbb {N}\) be coarse respectively fine time discretization parameters with corresponding timestep sizes \(\tau _c\), \(\tau _f\). For simplicity, we assume \(M_c /M_f = r \in \mathbb {N}\). Now, the coarse intervals are generated by the fine ones, i.e.
The coarse averaging operator can be decomposed into the fine averaging operator
To accurately substitute the analytic averaging operator, we define the discrete time averaging operator
Additionally, we define the error quantities
and
5.3 Joint sampling on fine and coarse scales
It is crucial to use the same stochastic input when computing \({\textbf{v}}_f\) and \({\textbf{v}}_c\). This can be done in two different ways. Either one first samples coarse stochastic data and a posteriori samples the fine data based on the conditional probabilities of the coarse one, or we can sample the fine stochastic input and try to reconstruct the coarse stochastic data. The latter approach is more suitable for the averaged increments.
Lemma 37
Let \(M_f, M_c \in \mathbb {N}\). Assume \(M_f = r M_c\) for some \(r \in \mathbb {N}\). Then
for \(j \in {\{{2,\ldots , M_c}\}}\).
Proof
Since \(M_f = r M_c\), we have \(\tau _c = r \tau _f\). Thus,
Due to the discrete Fubini’s theorem
The second term is easily computed
Overall,
Let \(j \in {\{{2,\ldots ,M_c}\}}\). Since \(t_{j}^c = t_{rj}^f\), it holds
The discrete Fubini’s theorem implies
Therefore,
The claim is proved. \(\square \)
Remark 38
The reconstruction formula (5.14) is the key ingredient, why it is possible to compare fine and coarse numerical solutions in an efficient way. If one tries to establish a corresponding formula for randomized algorithms as proposed in [18], this task becomes more challenging.
5.4 Simulation: unknown solution
Let \(\mathcal {O} = (0,1)^2\) and \(T = 1\). We choose \(p \in {\{{1.5,3}\}}\), \(u_0(x) = \sin (\pi x_1) \sin (\pi x_2)\) and
\(V_h\) denotes the space of piece wise linear elements with zero boundary values on a triangulation \(\mathcal {T}_h\) of \(\mathcal {O}\). We initialize the coarse triangulation \(\mathcal {T}_{h_c}\) as a uniform triangulation of \(\mathcal {O}\) such that \(\left|V_{h_c}\right| = 121\) and \(h_c \approx 1.4*10^{-1}\). \(\mathcal {T}_{h_f}\) is generated by three uniform refinements of \(\mathcal {T}_{h_c}\). Then \(\left|V_{h_f}\right| = 6561\) and \(h_f \approx 1.7*10^{-2}\). For the time discretization we use \(M_f = 1280\) and \(M_c = 40\). Therefore \(\tau _f \approx 7.8*10^{-4}\) and \(\tau _c \approx 2.5* 10^{-2}\). We measure the error of the fine numerical solution \({\textbf{v}}_f \in {\{{{\textbf{v}}_f^{\text {EM}},{\textbf{v}}_f^{\text {Half}},{\textbf{v}}_f^{\text {Full}}}\}}\) versus the coarse numerical solution \({\textbf{v}}_c \in {\{{{\textbf{v}}_c^{\text {EM}},{\textbf{v}}_c^{\text {Half}},{\textbf{v}}_c^{\text {Full}}}\}} \) of the same algorithm in the error quantities (5.12) and (5.13). The expectation is approximated by the Monte-Carlo method with 20 samples.
In Fig. respectively Fig. we plot the evolution of the error for \(p=3\) respectively \(p = 1.5\). In both cases we observe linear convergence. This indicates that on the used discretization scale the space error dominates the time error. We do not see a substantial difference in the gradient error quantities. Additionally, \({\textbf{v}}^{\text {EM}}\) measured in the point distance \(d_{L^\infty L^2}^{\,\text {point}}\) and \({\textbf{v}}^{\text {Half}}\) measured in the averaged distance \(d_{L^\infty L^2}^{\,\text {aver}}\) perform equally well.
If \(p=3\) then \({\textbf{v}}^{\text {Full}}\) behaves similarly in both error terms and performs slightly worse than \({\textbf{v}}^{\text {EM}}\) and \({\textbf{v}}^{\text {Half}}\). Contrary in the case \(p=1.5\), while still performing slightly worse compared to \({\textbf{v}}^{\text {EM}}\) and \({\textbf{v}}^{\text {Half}}\), \({\textbf{v}}^{\text {Full}}\) is better approximated in \(d_{L^\infty L^2}^{\,\text {aver}}\) than in \(d_{L^\infty L^2}^{\,\text {point}}\). This indicates that at least in the singular setting \({\textbf{v}}^{\text {Full}}\) also approximates average values.
Convergence for \(p=3\) of \({\textbf{v}}_c^{\text {EM}}\) (red), \({\textbf{v}}_c^{\text {Half}}\) (blue) and \({\textbf{v}}_c^{\text {Full}}\) (green) towards \({\textbf{v}}_f\) measured in \(d_{L^\infty L^2}^{\,\text {aver}}\) (solid, left), \(d_{L^\infty L^2}^{\,\text {point}}\) (dashed, left), \(d_{L^2 V}^{\,\text {classic}}\) (solid, right), \(d_{L^2 V}^{\,\text {inner}}\) (dashed, right) and \(d_{L^2 V}^{\,\text {outer}}\) (dash dotted, right). In the top row we use \(\tau \approx h\) and in the bottom row \(\tau \approx h^2\) (color figure online)
Convergence for \(p=1.5\) of \({\textbf{v}}_c^{\text {EM}}\) (red), \({\textbf{v}}_c^{\text {Half}}\) (blue) and \({\textbf{v}}_c^{\text {Full}}\) (green) towards \({\textbf{v}}_f\) measured in \(d_{L^\infty L^2}^{\,\text {aver}}\) (solid, left), \(d_{L^\infty L^2}^{\,\text {point}}\) (dashed, left), \(d_{L^2 V}^{\,\text {classic}}\) (solid, right), \(d_{L^2 V}^{\,\text {inner}}\) (dashed, right) and \(d_{L^2 V}^{\,\text {outer}}\) (dash dotted, right). In the top row we use \(\tau \approx h\) and in the bottom row \(\tau \approx h^2\) (color figure online)
5.5 Conclusion
Our algorithms achieve optimal linear convergence in space and optimal 1/2-convergence in time with minimal regularity assumptions. Experimentally, the Euler–Maruyama scheme (5.1) and our algorithms (3.14) and (3.15) share the same computational complexity and accuracy. However, the convergence of the Euler–Maruyama scheme is only proven with rate 1/4, cf. [39]. The construction of the random inputs is fast and can be done with the simple sampling algorithm (4.9).
Notes
denotes equality in distribution.
References
Andersson, A., Larsson, S.: Weak convergence for a spatial approximation of the nonlinear stochastic heat equation. Math. Comput. 85(299), 1335–1358 (2016). https://doi.org/10.1090/mcom/3016
Balci, A.K., Cianchi, A., Diening, L., Maz’ya, V.: A pointwise differential inequality and second-order regularity for nonlinear elliptic systems. In: Mathematische Annalen (Aug. 2021). https://doi.org/10.1007/s00208-021-02249-9
Balci, A.K., Diening, L., Weimar, M.: Higher order Calderón–Zygmund estimates for the \(p\)-Laplace equation. J. Differ. Equ. 268(2), 590–635 (2020). https://doi.org/10.1016/j.jde.2019.08.009
Baldi, P.: Stochastic Calculus. Universitext. An Introduction Through Theory and Exercises, pp. xiv+627. Springer, Cham (2017). https://doi.org/10.1007/978-3-319-62226-2
Baňas, L’., Röckner, M., Wilke, A.: Convergent numerical approximation of the stochastic total variation flow. Stoch. Partial Differ. Equ. Anal. Comput. 9(2), 437–471 (2021). https://doi.org/10.1007/s40072-020-00169-4
Barrett, J.W., Liu, W.B.: Finite element approximation of the \(p\)-Laplacian. Math. Comput. 61(204), 523–537 (1993). https://doi.org/10.2307/2153239
Barrett, J.W., Liu, W.B.: Finite element approximation of the parabolic \(p\)-Laplacian. SIAM J. Numer. Anal. 31(2), 413–428 (1994). https://doi.org/10.1137/0731022
Bartels, S., Diening, L., Nochetto, R.H.: Unconditional stability of semiimplicit discretizations of singular flows. SIAM J. Numer. Anal. 56(3), 1896–1914 (2018). https://doi.org/10.1137/17M1159166
Bauzet, C., Vallet, G., Wittbold, P., Zimmermann, A.: On a \(p(t, x)\)-Laplace evolution equation with a stochastic force. Stoch. Partial Differ. Equ. Anal. Comput. 1(3), 552–570 (2013). https://doi.org/10.1007/s40072-013-0017-z
Belenki, L., Diening, L., Kreuzer, C.: Optimality of an adaptive finite element method for the \(p\)-Laplacian equation. IMA J. Numer. Anal. 32(2), 484–510 (2012). https://doi.org/10.1093/imanum/drr016
Berselli, L.C., R\(\mathring{{\rm u}}\)žička, M.: Space–time discretization for nonlinear parabolic systems with \(p\)-structure. IMA J. Numer. Anal. 42(1), 260–299 (2022) https://doi.org/10.1093/imanum/draa079
Bertini, L., Cancrini, N.: The stochastic heat equation: Feynman–Kac formula and intermittence. J. Stat. Phys. 78(5–6), 1377–1401 (1995). https://doi.org/10.1007/BF02180136
Bojarski, B., Iwaniec, T.: p-Harmonic equation and quasiregular mappings. Banach Center Publications. 1(19), 25–38 (1987)
Breit, D.: Existence theory for stochastic power law fluids. J. Math. Fluid Mech. 17(2), 295–326 (2015). https://doi.org/10.1007/s00021-015-0203-z
Breit, D.: Regularity theory for nonlinear systems of SPDEs. Manuscr. Math. 146(3–4), 329–349 (2015). https://doi.org/10.1007/s00229-014-0704-8
Breit, D., Diening, L., Storn, J., Wichmann, J.: The parabolic \(p\)-Laplacian with fractional differentiability. IMA J. Numer. Anal. 41(3), 2110–2138 (2021). https://doi.org/10.1093/imanum/draa081
Breit, D., Gmeineder, F.: Electro-rheological fluids under random influences: martingale and strong solutions. Stoch. Partial Differ. Equ. Anal. Comput. 7(4), 699–745 (2019). https://doi.org/10.1007/s40072-019-00138-6
Breit, D., Hofmanová, M., Loisel, S.: Space–time approximation of stochastic \(p\)-Laplace-type systems. SIAM J. Numer. Anal. 59(4), 2218–2236 (2021). https://doi.org/10.1137/20M1334310
Breit, D., Mensah, P.R.: Space–time approximation of parabolic systems with variable growth. IMA J. Numer. Anal. (2019). https://doi.org/10.1093/imanum/drz039
Cianchi, A., Maz’ya, V.G.: Second-order two-sided estimates in nonlinear elliptic problems. Arch. Ration. Mech. Anal. 229(2), 569–599 (2018). https://doi.org/10.1007/s00205-018-1223-7
Cianchi, A., Maz’ya, V.G.: Optimal second-order regularity for the \(p\)-Laplace system. J. Math. Pures Appl. 9(132), 41–78 (2019). https://doi.org/10.1016/j.matpur.2019.02.015
Cianchi, A., Maz’ya, V.G.: Second-order regularity for parabolic p-Laplace problems. J. Geom. Anal. (2020). https://doi.org/10.1007/s12220-019-00213-3
Dareiotis, K., Gerencsér, M.: On the regularisation of the noise for the Euler–Maruyama scheme with irregular drift. Electron. J. Probab. 25, Paper No. 82, 18 (2020). https://doi.org/10.1214/20-ejp479
Debussche, A., Printems, J.: Weak order for the discretization of the stochastic heat equation. Math. Comput. 78(266), 845–863 (2009). https://doi.org/10.1090/S0025-5718-08-02184-4
Di Pietro, D.A., Droniou, J., Harnist, A.: Improved error estimates for hybrid high-order discretizations of Leray–Lions problems. Calcolo 58(2), Paper No. 19, 24 (2021). https://doi.org/10.1007/s10092-021-00410-z
DiBenedetto, E.: \(C^{1+\alpha }\) local regularity of weak solutions of degenerate elliptic equations. Nonlinear Anal. 7(8), 827–850 (1983). https://doi.org/10.1016/0362-546X(83)90061-5
Diening, L., Fornasier, M., Tomasi, R., Wank, M.: A relaxed Kacanov iteration for the \(p\)-Poisson problem. Numer. Math. 145(1), 1–34 (2020). https://doi.org/10.1007/s00211-020-01107-1
Diening, L., Kaplický, P., Schwarzacher, S.: BMO estimates for the \(p\)-Laplacian. Nonlinear Anal. 75(2), 637–650 (2012). https://doi.org/10.1016/j.na.2011.08.065
Diening, L., R\(\mathring{{\rm u}}\)žička, M.: Interpolation operators in Orlicz–Sobolev spaces. Numer. Math. 107(1), 107–129 (2007). https://doi.org/10.1007/s00211-007-0079-9
Diening, L., Ebmeyer, C., R\(\mathring{{\rm u}}\)žička, M.: Optimal convergence for the implicit space–time discretization of parabolic systems with \(p\)-structure. SIAM J. Numer. Anal. 45(2), 457–472 (2007). https://doi.org/10.1137/05064120X
Diening, L., Ettwein, F.: Fractional estimates for non-differentiable elliptic systems with general growth. Forum Math. 20(3), 523–556 (2008). https://doi.org/10.1515/FORUM.2008.027
Diening, L., Harjulehto, P., Hästö, P., R\(\mathring{{\rm u}}\)žička, M.: Lebesgue and Sobolev Spaces with Variable Exponents. Vol. 2017. Lecture Notes in Mathematics, pp. x+509. Springer, Heidelberg (2011). https://doi.org/10.1007/978-3-642-18363-8
Diening, L., Kreuzer, C.: Linear convergence of an adaptive finite element method for the \(p\)-Laplacian equation. SIAM J. Numer. Anal. 46(2), 614–638 (2008). https://doi.org/10.1137/070681508
Diening, L., Stroffolini, B., Verde, A.: Everywhere regularity of functionals with \(\varphi \)-growth. Manuscr. Math. 129(4), 449–481 (2009). https://doi.org/10.1007/s00229-009-0277-0
Droniou, J., Goldys, B., Le, K.-N.: Design and convergence analysis of numerical methods for stochastic evolution equations with Leray–Lions operator. IMA J. Numer. Anal. 42(2), 1143–1179 (2022). https://doi.org/10.1093/imanum/draa105
Dunst, T., Prohl, A.: The forward-backward stochastic heat equation: numerical analysis and simulation. SIAM J. Sci. Comput. 38(5), A2725–A2755 (2016). https://doi.org/10.1137/15M1022951
Ebmeyer, C., Liu, W.: Quasi-norm interpolation error estimates for the piecewise linear finite element approximation of \(p\)-Laplacian problems. Numer. Math. 100(2), 233–258 (2005). https://doi.org/10.1007/s00211-005-0594-5
Eisenmann, M., Kovács, M., Kruse, R., Larsson, S.: On a randomized backward Euler method for nonlinear evolution equations with time-irregular coefficients. Found. Comput. Math. 19(6), 1387–1430 (2019). https://doi.org/10.1007/s10208-018-09412-w
Eisenmann, M., Kovács, M., Kruse, R., Larsson, S.: Error estimates of the backward Euler–Maruyama method for multi-valued stochastic differential equations. BIT Numer. Math. 62(3), 803–848 (2022). https://doi.org/10.1007/s10543-021-00893-w
Emmrich, E., Šiška, D.: Nonlinear stochastic evolution equations of second order with damping. Stoch. Partial Differ. Equ. Anal. Comput. 5(1), 81–112 (2017). https://doi.org/10.1007/s40072-016-0082-1
Gajewski, H., Gröger, K., Zacharias, K.: Nichtlineare Operatorgleichungen und Operatordifferentialgleichungen. Mathematische Lehrbücher und Monographien, II. Abteilung, Mathematische Monographien, Band 38. Akademie, Berlin, ix+281 pp. (loose errata) (1974)
Geissert, M., Kovács, M., Larsson, S.: Rate of weak convergence of the finite element method for the stochastic heat equation with additive noise. BIT 49(2), 343–356 (2009). https://doi.org/10.1007/s10543-009-0227-y
Gess, B.: Strong solutions for stochastic partial differential equations of gradient type. J. Funct. Anal. 263(8), 2355–2383 (2012). https://doi.org/10.1016/j.jfa.2012.07.001
Gess, B., Tölle, J.M.: Ergodicity and local limits for stochastic local and nonlocal \(p\)-Laplace equations. SIAM J. Math. Anal. 48(6), 4094–4125 (2016). https://doi.org/10.1137/15M1049774
Gyöngy, I., Millet, A.: On discretization schemes for stochastic evolution equations. Potential Anal. 23(2), 99–134 (2005). https://doi.org/10.1007/s11118-004-5393-6
Gyöngy, I., Millet, A.: Rate of convergence of space time approximations for stochastic evolution equations. Potential Anal. 30(1), 29–64 (2009). https://doi.org/10.1007/s11118-008-9105-5
Hairer, M., Labbé, C.: Multiplicative stochastic heat equations on the whole space. J. Eur. Math. Soc.: JEMS 20(4), 1005–1054 (2018). https://doi.org/10.4171/JEMS/781
Hu, Y., Huang, J., Nualart, D., Tindel, S.: Stochastic heat equations with general multiplicative Gaussian noises: Hölder continuity and intermittency. Electron. J. Probab. 20(55), 50 (2015). https://doi.org/10.1214/EJP.v20-3316
Hutzenthaler, M., Jentzen, A., Kloeden, P.E.: Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with non-globally Lipschitz continuous coefficients. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 467(2130), 1563–1576 (2011). https://doi.org/10.1098/rspa.2010.0348
Hytönen, T.P., Veraar, M.C.: On Besov regularity of Brownian motions in infinite dimensions. Probab. Math. Stat. 28(1), 143–162 (2008)
Iwaniec, T., Manfredi, J.J.: Regularity of \(p\)-harmonic functions on the plane. Rev. Mat. Iberoam. 5(1–2), 1–19 (1989). https://doi.org/10.4171/RMI/82
Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic Differential Equations. Vol. 23. Applications of Mathematics (New York). Springer, Berlin, pp. xxxvi+632 (1992). https://doi.org/10.1007/978-3-662-12616-5
Kuusi, T., Mingione, G.: Linear potentials in nonlinear potential theory. Arch. Ration. Mech. Anal. 207(1), 215–246 (2013). https://doi.org/10.1007/s00205-012-0562-z
Kuusi, T., Mingione, G.: Vectorial nonlinear potential theory. J. Eur. Math. Soc.: JEMS 20(4), 929–1004 (2018). https://doi.org/10.4171/JEMS/780
Lê, K., Ling, C.: Taming singular stochastic differential equations: a numerical method. (2021). https://doi.org/10.48550/ARXIV.2110.01343
Lieberman, G.M.: Boundary regularity for solutions of degenerate elliptic equations. Nonlinear Anal. 12(11), 1203–1219 (1988). https://doi.org/10.1016/0362-546X(88)90053-3
Lions, J.-L.: Quelques méthodes de résolution des problèmes aux limites non linéaires, pp. xx+554. Dunod, Paris (1969)
Liu, W.B., Barrett, J.W.: A further remark on the regularity of the solutions of the \(p\)-Laplacian and its applications to their finite element approximation. Nonlinear Anal. 21(5), 379–387 (1993). https://doi.org/10.1016/0362-546X(93)90081-3
Liu, W.B., Barrett, J.W.: A remark on the regularity of the solutions of the \(p\)-Laplacian and its application to their finite element approximation. J. Math. Anal. Appl. 178(2), 470–487 (1993). https://doi.org/10.1006/jmaa.1993.1319
Liu, W.: On the stochastic \(p\)-Laplace equation. J. Math. Anal. Appl. 360(2), 737–751 (2009). https://doi.org/10.1016/j.jmaa.2009.07.020
Liu, W., Röckner, M.: SPDE in Hilbert space with locally monotone coefficients. J. Funct. Anal. 259(11), 2902–2922 (2010). https://doi.org/10.1016/j.jfa.2010.05.012
Logg, A., Mardal, K.-A., Wells, G.N. (eds.): Automated Solution of Differential Equations by the Finite Element Method. Vol. 84. Lecture Notes in Computational Science and Engineering. The FEniCS book, pp. xiv+723. Springer, Heidelberg (2012). https://doi.org/10.1007/978-3-642-23099-8
Mansuy, R., Yor, M.: Aspects of Brownian Motion. Universitext. Springer, Berlin, pp. xiv+195 (2008). https://doi.org/10.1007/978-3-540-49966-4
Marcellini, P., Papi, G.: Nonlinear elliptic systems with general growth. J. Differ. Equ. 221(2), 412–443 (2006). https://doi.org/10.1016/j.jde.2004.11.011
Ondreját, M., Prohl, A., Walkington, N.J.: Numerical approximation of nonlinear SPDE’s. In: Stochastics and Partial Differential Equations: Analysis and Computations (2022). https://doi.org/10.1007/s40072-022-00271-9
Prévôt, C., Röckner, M.: A Concise Course on Stochastic Partial Differential Equations. Vol. 1905. Lecture Notes in Mathematics, pp. vi+144. Springer, Berlin (2007)
Prohl, A., Wang, Y.: Strong rates of convergence for a space–time discretization of the backward stochastic heat equation, and of a linear-quadratic control problem for the stochastic heat equation. ESAIM Control Optim. Calc. Var. 27, Paper No. 54, 30 (2021). https://doi.org/10.1051/cocv/2021052
Sapountzoglou, N., Zimmermann, A.: Well-posedness of renormalized solutions for a stochastic \(p\)-Laplace equation with \(L^{1}\)-initial data. Discrete Contin. Dyn. Syst. 41(5), 2341–2376 (2021). https://doi.org/10.3934/dcds.2020367
Terasawa, Y., Yoshida, N.: Stochastic power law fluids: existence and uniqueness of weak solutions. Ann. Appl. Probab. 21(5), 1827–1859 (2011). https://doi.org/10.1214/10-AAP741
Tolksdorf, P.: Regularity for a more general class of quasilinear elliptic equations. J. Differ. Equ. 51(1), 126–150 (1984). https://doi.org/10.1016/0022-0396(84)90105-0
Uhlenbeck, K.: Regularity for a class of non-linear elliptic systems. Acta Math. 138(3–4), 219–240 (1977). https://doi.org/10.1007/BF02392316
Ural’ceva, N.N.: Degenerate quasilinear elliptic systems. Zap. Naucn. Sem. Leningrad. Otdel. Mat. Inst. Steklov. (LOMI) 7, 184–222 (1968)
Vallet, G., Wittbold, P., Zimmermann, A.: On a stochastic \(p(w, t, x)\)-Laplace equation. In: Thirteenth International Conference Zaragoza-Pau on Mathematics and its Applications. Vol. 40. Monogr. Mat. García Galdeano. Prensas Univ. Zaragoza, Zaragoza, pp. 125–134 (2016)
Vallet, G., Zimmermann, A.: The stochastic \(p(w, t, x)\)-Laplace equation with cylindrical Wiener process. J. Math. Anal. Appl. 444(2), 1359–1371 (2016). https://doi.org/10.1016/j.jmaa.2016.07.018
Wei, D.: Existence, uniqueness, and numerical analysis of solutions of a quasilinear parabolic problem. SIAM J. Numer. Anal. 29(2), 484–497 (1992). https://doi.org/10.1137/0729029
Wichmann, J.: On temporal regularity for strong solutions to stochastic \(p\)-Laplace systems (2021). arXiv: 2111.09601 [math.AP]
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.
The research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – SFB 1283/2 2021 – 317210226.
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
Diening, L., Hofmanová, M. & Wichmann, J. An averaged space–time discretization of the stochastic p-Laplace system. Numer. Math. 153, 557–609 (2023). https://doi.org/10.1007/s00211-022-01343-7
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00211-022-01343-7