An averaged space–time discretization of the stochastic p-Laplace system

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.

1.2.Existence & well-posedness.The existence of analytically weak solutions in the space L 2 Ω; C [0, T ]; L 2 (O) ∩ L p Ω; L p 0, T ; W 1,p 0 (O) can be established by standard monotonicity arguments [LR10].It requires a linear growth assumption on the noise coefficient G and L 2 (O)-integrable initial data.
First results on the existence of strong solutions to stochastically perturbed gradient flows have been obtained by Gess [Ges12].It includes the degenerate, p ≥ 2, p-Laplace equation.
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.
Sobolev regularity provides an alternative scale of smoothness.In the singular case, p ∈ (1, 2), W 2,2 -regularity has been proven by Liu and Barrett [LB93b].The degenerate setting is more delicate and one can not expect ∇ 2 u to be well-defined on {∇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 [IM89].
The nonlinear character of the equation naturally introduces the additional quantities S(∇u) and V (∇u).In the scalar case both expressions are robust with respect to p ∈ (1, ∞).Here it is possible to prove V (∇u) ∈ W 1,2 via a difference quotient technique [BI87;LB93a] and S(∇u) ∈ W 1,2 by a functional inequality [CM18].The result V (∇u) ∈ W 1,2 generalizes to the vectorial case.However, the functional inequality fails in the vectorial case at least point-wise for p ≤ 2(2 − √ 2) ≈ 1.1715.Therefore, it is unclear whether a regularity result S(∇u) ∈ W 1,2 is achievable for small p.Nevertheless, for p > 2(2 − √ 2) it is shown in [BCDM21] that S(∇u) ∈ W 1,2 .Regularity for S(∇u) on the Besov and Triebel-Lizorkin scale in the plane for p ≥ 2 has been obtained in [BDW20].
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 [Bre15b].However, the time regularity is limited due to the presence of the stochastic forcing.In the degenerate case partial time regularity can be recovered by exploiting the strong formulation of the p-Laplace system as done by Wichmann [Wic21].Overall, for appropriate data assumptions it is possible to verify (see Section 2.5) Here Φ 2 (s) := e s 2 − 1 and B denotes a Besov-Orlicz space (see Section 2.1).1.4.Approximation.In the past many authors have studied the numerical approximation of the deterministic counterpart of (1.1), e.g.[Wei92; BL93; BL94; EL05; DER07; DK08; BDK12; BDN18; BM19; BDSW21; BR22].
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 ξ and u m,h is an approximation of u(t m ).In [DK08] it has been proved that the expression V (∇u) − V (∇u h ) This settles the question about optimal convergence for piece-wise linear continuous elements.In fact, the paper [DER07] deals with the parabolic system and optimal convergence under the regularity assumption (1.3) has been achieved, i.e., (1.6) The error quantity (1.6) relies on the fact that the mapping t → V (∇u(t)) L 2 (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 welldefined error quantities.
In [BDSW21] two of us develop a numerical scheme based on time averages to circumvent the usage of point evaluations.A more probabilistic approach has been used in [EKKL19].There the authors replace deterministic evaluation points by random ones.
First results for monotone stochastic equations were derived by Gyöngy and Millet [GM05;GM09].They developed an abstract discretization theory that also covers the system (1.1).However, they only deal with the degenerate case and 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 [BHL21] use randomized time steps to construct an algorithm that achieves almost optimal convergence in time and optimal convergence in space, i.e., for all α ∈ (0, 1/2), Here E t denotes the expectation with respect to the random time steps.They use Hölder and Sobolev-Slobodeckij spaces to measure the time regularity of the solution.This excludes the limiting case α = 1/2.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, P-a.s. for all m ≥ 2 and The randomness enters through the averaged increments 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 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 ∆ m W. We show in Corollary 32 that ∆ m W is a Gaussian random variable whose variance is slightly reduced compared to the classical increments ∆ 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 (∆ m W, ∆ m W) M m=1 .
1.6.Outline.In Section 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, Section 4 clarifies the construction of the discrete random input vectors and provides a sampling algorithm.Lastly, Section 5 contains numerical experiments.

Mathematical setup
This section contains classical definitions and preliminary results.It is structured as follows: Section 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 Section 2.4.Lastly, Section 2.5 collects regularity results.
Let O ⊂ R n for n ≥ 1 be a bounded Lipschitz domain (further assumptions on 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 O T := I × O for the time space cylinder.Moreover let (Ω, F, (F t ) t∈I , P) denote a stochastic basis, i.e. a probability space with a complete and right continuous filtration (F t ) t∈I .We write f g for two non-negative quantities f and g if f is bounded by g up to a multiplicative constant.Accordingly we define and .We denote by c a generic constant which can change its value from line to line.
2.1.Function spaces.As usual, L q (O) denotes the Lebesgue space and W 1,q (O) the Sobolev space, where 1 ≤ q < ∞.We denote by W 1,q 0 (O) the Sobolev spaces with zero boundary values.It is the closure of C ∞ 0 (O) (smooth functions with compact support) in the W 1,q (O)-norm.We denote by W −1,q (O) the dual of W 1,q 0 (O).We do not distinguish in the notation between vector-and matrix-valued functions.
For a Banach space (X, • X ) let L q (I; X) be the Bochner space of Bochnermeasurable functions u : I → X satisfying t → u(t) X ∈ L q (I).Moreover, C(I; X) is the space of continuous functions with respect to the norm-topology.We also use C α (I; X) for the space of Hölder continuous functions.Given an Orliczfunction Φ : [0, ∞] → [0, ∞], i.e. a convex function satisfying lim t→0 Φ(t)/t = 0 and lim t→∞ Φ(t)/t = ∞ we define the Luxemburg-norm u L Φ (I;X) := inf λ > 0 : The Orlicz space L Φ (I; X) is the space of all Bochner-measurable functions with finite Luxemburg-norm.For more details on Orlicz-spaces we refer to [DHHR11].Given h ∈ I and u : I → X we define the difference operator τ h : {u : I → X} → {u : I ∩ I − {h} → X} via τ h (u)(s) := u(s + h) − u(s).The Besov-Orlicz space B α Φ,r (I; X) with differentiability α ∈ (0, 1), integrability Φ and fine index r ∈ (1, ∞] is defined as the space of Bochner-measurable functions with finite Besov-Orlicz norm • B α Φ,r (I;X) , where In the case r = ∞ the integral in h is replaced by an essential supremum and the space is commonly called Nikolskii-Orlicz space.When Φ(t) = t p for some p ∈ (1, ∞) we call the space B α Φ,r (I; X) = B α p,r (I; X) Besov space.Similarly, given a Banach space (Y, • Y ), we define L q (Ω; Y ) as the Bochner space of Bochner-measurable functions u : Ω → Y satisfying ω → u(ω) Y ∈ L q (Ω).The space L q F (Ω × I; X) denotes the subspace of X-valued (F t ) t∈I -progressively measurable processes.Let (U, • U ) be a separable Hilbert space.L 2 (U ; L 2 (O)) denotes the space of Hilbert-Schmidt operators from U to L 2 (O) with the norm where {u j } j∈N is some orthonormal basis of U .We abbreviate the notation L q ω L q t L q x := L q (Ω; L q (I; L q (O))) and L q− := 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 (U, • U ) be a separable Hilbert space.We assume that W is an U -valued cylindrical Wiener process with respect to (F t ) t∈I .Formally W can be represented as where {u j } j∈N is an orthonormal basis of U and {β j } j∈N are independent 1-dimensional standard Brownian motions.
Now it is a classical result, see for example the book of Prévôt and Röckner [PR07], that we can construct a corresponding stochastic integral.
Proposition 2. Let Assumption 1 be true.Then the operator I defined through x -valued martingale with respect to (F t ) t∈I , (b) (Itô isometry) for all t ∈ I it holds where ϕ(t) := t 0 (κ + s) p−2 s ds.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.[DE08].For more details we refer to [DR07; BDK12; DFTW20].
Given some initial condition u 0 : Ω × O → 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 J : (2.12) 2.4.Weak and strong solutions.We fix the concept of solutions as follows.
The process u is called strong solution if it is a weak solution and additionally satisfies x , (b) for all t ∈ I and P − a.s. it holds 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 [LR10].In fact, Assumption 1 is sufficient to guarantee 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 In [Ges12] 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.

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 [Wic21] that the strong solution enjoys 1/2-time differentiability in an exponential Besov space and even the nonlinear greadient V (∇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 ∈ ∂O, Let the assumptions of Theorem 8 be satisfied.Additionally, assume (2.17).Let u be the strong solution to (2.11).Then Spatial regularity is closely connected to the existence of strong solutions.Local regularity has been obtained in [Bre15b].
x be 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 [BCDM21; CM19], i.e.
Precisely, given some bounded Lipschitz domain O such that ∂O ∈ W 2,1 , i.e.O is locally the subgraph of a Lipschitz continuous function of n − 1 variables, which is also twice weakly differentiable.Denote by B the second fundamental form on ∂O, by |B| its norm and define where B r (x) denotes the ball of radius r around x, cap B1(x) (E) is the capacity of the set E relative to the ball B 1 (x) and H n−1 is the n − 1 dimensional Hausdorff measure.

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 O ⊂ R n be a bounded polyhedral domain.By T h denote a regular partition (triangulation) of O (no hanging nodes), which consists of closed n-simplices called elements.For each element (n-simplex) T ∈ T h we denote by h T the diameter of T , and by ρ T the supremum of the diameters of inscribed balls.
We assume that T h is shape regular, that is there exists a constant γ (the shape regularity constant) such that max We define the maximal mesh-size by We assume further that our triangulation is quasi-uniform, i.e.
For ∈ N 0 we denote by P (O) the polynomials on O of degree less than or equal to .
For fixed r ∈ N we define the vector valued finite element space V h as We will need some classical results on the stability properties of the L 2 x -projection for finite elements.
Lemma 13.Let r ≥ 1, V h be given by (3.3) and 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 ([BDSW21], Theorem 7).Let r ≥ 1, V h be given by (3.3) and T h be quasi-uniform.Then By |I m | we denote the Lebesgue measure of I m .By − Im g ds we denote the mean value integral over the set I m .We also abbreviate g m = − Im g ds for the mean value.
Let us discuss the stability properties of the piecewise constant approximation process generated by the mean values in terms of Nikolskii spaces.
A change of variables t = s + h and Fubini's Theorem yield x ds dh The proof is finished.
Remark 16.Lemma 15 is also valid for r = ∞.Here we need to substitute the left hand side in (3.4) by is the space of α-Hölder continuous functions.In our application the process u only belongs to Therefore we need to adjust the stability result in terms of exponentially integrable Nikolskii spaces. (3.5) Due to the embedding B α Φ2,∞ → C t for all α > 0 we find u to be continuous as an L 2 x -valued process.Therefore, it holds lim Jensen's inequality implies In particular, since Thus, It follows Therefore, (3.7) The assertion follows by applying the square root in (3.6) and using the estimate (3.7).
3.3.The averaged algorithm.Let u be the strong solution of (2.11), i.e. for all t ∈ I and 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 − τ and take the average over We denote by (•, •) the L 2 x -inner product.The corresponding weak formulation reads for all ξ ∈ L 2 x ∩ W 1,p 0,x (3.11) 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 ∈ {2, . . ., M } is determinded via for all ξ h ∈ V h and 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 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 [EŠ17] Section 3 for a detailed analysis of a more general system.
3.4.Error analysis.We are ready to state and prove the main result.
Theorem 19 (Convergence of Algorithm (3.14)).Let u 0 ∈ L 2 ω L 2 x be F 0 -measurable and Assumption 1 be satisfied.Additionally, assume T h to be quasi-uniform.Moreover, let u be the weak solution to (2.11) and v be the numerical solution of (3.14).
(a) Assume (3.17) 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 [HV08] 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 ∈ [1, ∞) and P-a.s.
The logarithmic term in (3.19) has already been used in the context of rough stochastic differential equations [DG20; LL21].It quantifies the distance of L Φ2 and L ∞ .
Theorem 19 can be generalized to cover fractional time and space regularity assumptions as done for the deterministic p-Laplace system in [BDSW21].

Now use the symmetry of the L 2
x -projection and the algebraic identity The second term, due to the V -coercivity (2.8), Take the maximum over m * ∈ {2, . . ., M } and expectation (3.20) 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 where e 0 := u 0 − v 0 .By (3.14a) we have Π 2 e 0 = 0. Therefore, (3.21) 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 ≤ 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) (3.23) 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 ≤ 1 Decomposition in time and space error, applying Lemma 13 and Lemma 15 and the estimate (3.22) (3.24) 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 F m := F tm .Note, (3.25) does not define a martingale with respect to F tm .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 F m * , i.e., Due to the tower property of conditional expectation, the measurability of e m with respect to F m together with the martingale property of the stochastic integral Again using the F m -measurability of e m , we conclude that K m * is F m * -measurable.Thus, It remains to compute the conditional expectation in M 2 .Since t m * is an interior point of I m * +1 ∪ 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 F m * .Overall, 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.27) (3.26) together with (3.27) allow to identify increments of the conditional expectations, Observe E K M F 1 = 0, since Π 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 Together we can estimate This concludes the bound for J 2 in (3.20) (3.28) 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 ε sufficiently small and applying Gronwall's Lemma ensures This implies (3.29) 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 Additionally, Jensen's inequality ensures x dν dt .
Therefore, invoking (3.29), The assertion (3.17) follows by an application of Lemma 13 x .We apply Lemma 17 to bound This verifies (3.19) and the proof is finished.
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 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 Lemma 3 and Lemma 4 Absorbing the second term to the left hand side and applying the bound (3.17) verifies the assertion.
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).
x . (3.33) Proof.The equation (3.32) follows by using The estimate (3.33) is obtained trivially, Remark 24.In [DKS12, Lemma 6.2] the authors prove the equivalence (although it is done purely in space but can be extended to time) of Theorem 25 (Convergence of Algorithm 3.15).Let the assumptions of Theorem 19 be satisfied.Denote by w ∈ (V h ) M +1 the solution to (3.15) and by u the weak solution to (2.11).Then 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 u 1 , we rather choose u I1∪I2 := − I1∪I2 u ν dν.The equation for the latter one reads (3.36) Subtracting (3.15b) from (3.36) and choosing ξ h = Π 2 u I1∪I2 − w 1 results in Due to Hölder's and Young's inequalities, Itô isometry and the growth assumption (2.2) x dν .
The boundedness of u as an L 2 x valued-process implies The forth term is estimated similarly, x .Using the L 2 -stability of the L 2 -projection we obtain x .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 x .Lastly, using Lemma 13 and Lemma 15 x .The bound for the initial error is complete.

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 26.A stochastic process Y is called U -valued Gaussian process with mean operator m : I × U → R and variance operator Σ : In short, we write Y t ∼ N (m t , Σ t ) 1 .
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 Σ t (u, v) := t(u, v) U for all t ∈ I and u, v ∈ U .Moreover, the series (2.1) converges in the weak topology 1 denotes equality in distribution of U , due to the hypercontractivity of normally distributed random variables it holds for any q > 0, u ∈ U and t, s ∈ 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 ( W m ) M m=1 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 m=1 .We decompose W | Im into a Brownian bridge B m and its nodal values W (t m−1 ) and W (t m ), i.e., for t ∈ I m where i.e. all finite dimensional distributions of the generators of each sigma algebra are independent of each other.
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 W m or B m m .The formula (4.4) provides an easy way to compute the remaining one.We choose to compute the distribution of W m .An application of Itô's formula for f (s, W s ) = s−tm τ W s implies P-a.s.
A stochastic integral that is driven by a Wiener process and a deterministic integrand stays Gaussian.
Note, in the case m = 1 we have ∆ 1 W = W 1 and the result follows by Corollary 30.
Let m ≥ 2 and u, v ∈ U .Due to the independence, So far we have identified how each averaged increment ∆ m W is distributed.However, we also need to know what the joint distribution is, i.e., the distribution of a random vector.
Lemma 33.The random vector (∆ m W) M m=1 is an U M -valued centered Gaussian random variable with variance operator Σ : (4.7) Proof.The equation (4.6) implies, that the random vector (∆ m W) M m=1 can be constructed via a linear transformation of independent Gaussian random vectors Therefore, (∆ m W) M m=1 is itself a centered Gaussian vector.It remains to compute the covariance matrix.Let u, v ∈ U and m, l ∈ {1, . . ., M }.If m = l Corollary 32 implies If |m − l| > 1, then equation (4.6) and the independence imply It remains to consider the case |m − l| = 1.Without loss of generality l = m + 1.Now, using (4.6), the independence and Corollary 31 The proof is finished.
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 where Z 0 := Z0 := 0. Let Π J be the U -orthogonal projection onto U J := span(u 1 , . . ., u J ).The following proposition guarantees that the sampling algorithm (4.9) approximates the desired random variables.
Proof.First, we need to observe that Z m ∼ Π J ∆ m W and Zm ∼ Π J B m m .Then, the statement follows similarly to the proof of Lemma 33.
Remark 35.Proposition 34 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.

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 v Half and the solution of (3.15) as v Full .Additionally, we compare our algorithms to the classical Euler-Maruyama discretization of (2.11).That reads, find v ∈ (V h ) M +1 such that for all ξ h ∈ V h , m ≥ 1 and P-a.s.
The solution to (5.1) is called v EM .So far it is unknown whether the approximation v EM converges to the solution u of (2.11) in a suitable sense if the regularity of the map t → V (∇u t ) is limited.Originally the Euler-Maruyama scheme has been introduced for stochastic differential equations.In this context much more is known and convergence has been proven, see e.g. the book of Kloeden and Platen [KP92, Section 9.5].However, when dropping the Lipschitz assumption on the coefficients, divergence with positive probability has been obtained in [HJK11].
We are particularly interested in the experimental study of the following questions: for some λ ∈ 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 becomes simpler.Let O = (0, 1) 2 , T = 1 and u 0 (x) = sin(πx 1 ) sin(πx 2 ).Note that u 0 is an eigenfunction of the 2-Laplacian with homogeneous Dirichlet data and corresponding eigenvalue µ = 2π 2 .The unique solution to (2.11) is given by u(ω, t, x) = exp − λ 2 2 + µ t + λβ 1 (ω, t) u 0 (x).(5.3)Similarly, we can give a closed expression for the solution u h to the space-discrete equation d (u h , ξ h ) + (∇u h , ∇ξ h ) dt = λ (u h , ξ h ) dβ 1 (t).
(5.4)This is equivalent to the system of linear stochastic differential equations dM u + Su dt = λM u dβ 1 (t), (5.5)where M i,j = ξ j h , ξ i h and S i,j = ∇ξ j h , ∇ξ i h is the mass respectively the stiffness matrix and {ξ i h } form a basis of V h .The eigenpairs of the discretized Laplacian are related to the linear system Su = µM u ⇔ M −1 Su = µu.
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 where u ∈ {u Point , u Aver } and v ∈ {v EM , v Full , v Half }.
In Figure 1 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 λ = 1 and V h to be the space of piecewise linear, continuous elements on a uniform mesh with |V h | = 121.The average values (5.10) are approximated by r = 10.
The numerical results support that v EM approximates the solution on the grid points, while v Full and v Half approximate the average values of the solution.The gap is due to the difference x .Initially, we observe a preasymptotic effect.It stabilizes at the time scale τ ≈ 10 −3 .Afterwards the predicted convergence speed of order 1 is achieved.5.2.Beyond known solutions.In general, a major obstacle is the absence of an analytic solution to the equation (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 v f and a coarse approximation v c .Both, v f and 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 ≥ h f > 0 be coarse respectively fine space mesh sizes.Similarly, let M c , M f ∈ N be x , (5.13b) x .(5.13c) 5.3.Joint sampling on fine and coarse scales.It is crucial to use the same stochastic input when computing v f and 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 nice independency properties.They do not look into the past nor future.The following result can be found in [MY08, Section 1.2].Proposition 27.Let (B m ) M m=1 be given by (4.2).Then for all m ∈ {1, . . ., M } σ B m (t) t ∈ I m ⊥ σ W (t) t ∈ [0, ∞)\(t m−1 , t m ) , (4.3) decomposition (4.6) we propose the following sampling algorithm.The time discretization is achieved by the parameter M ∈ N and the series truncation in (2.1) is done by the parameter J ∈ N. Given M, J ∈ N. (a) (Sampling) Compute i.i.d.random variables ζ j m , η j m ∼ N (0, 1) for m ∈ {1, . . ., M } and j ∈ {1, . . ., J}.(b) (Lift to Hilbert space U ) For m ∈ {1, . . ., M } define the random variables Z m := √ τ {u j } j∈N is an orthonormal system of U .(c) (Adjusting correlation) For m ∈ {1, . . ., M } define the random variables (a) Do the algorithms (3.14) and (3.15) approximate mean-values or pointvalues?(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 [LMW12].5.1.An explicit solution.In the linear case, p = 2, with a linear right hand side G(v)∆W = λv ∆β 1 , (5.2)

u
Point := (u h (t m )) the algorithms (3.14) and (3.15) approximate the mean value of the analytic solution, we define u Aver exact := ( u h Im ) M m=1 .(5.9)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 β 1 .Numerically, we substitute u h Im by a Riemann sum approximation, i.e. we fix an equidistant partition {[t m,k−1 , t m,k ]} r k=1 of I m with resolution r ∈ N and define u is measured in the error quantities E(u, v) := E max m=1,...,M u m − v m

10 − 4 Figure 1 ..
Figure1.Time convergence of v EM (red), v Half (blue) and v Full (green) towards u Point (dashed) respectively u Aver (solid) measured in the error terms E (left) and V (right).