A randomized and fully discrete Galerkin finite element method for semilinear stochastic evolution equations
Raphael Kruse Yue Wu
math.NA math.PR 60H15, 65C30, 65M12, 65M60


A randomized Galerkin finite element method for SEE]

A randomized and fully discrete

Galerkin finite element method for

semilinear stochastic evolution equations


R. Kruse]Raphael Kruse
Raphael Kruse

Technische Universität Berlin

Institut für Mathematik, Secr. MA 5-3

Straße des 17. Juni 136

DE-10623 Berlin

Germany
kruse@math.tu-berlin.de

Y. Wu]Yue Wu
Yue Wu

Technische Universität Berlin

Institut für Mathematik, Secr. MA 5-3

Straße des 17. Juni 136

DE-10623 Berlin

Germany
wu@math.tu-berlin.de

 
[2010]60H15, 65C30, 65M12, 65M60  


  In this paper the numerical solution of non-autonomous
  semilinear stochastic evolution equations driven by an additive
  Wiener noise is investigated. We introduce a novel fully discrete
  numerical approximation that combines a standard Galerkin finite element
  method with a randomized Runge–Kutta scheme. 
  Convergence of the method to the mild solution is proven with
  respect to the L^p-norm, p ∈ [2,∞).
  We obtain the same temporal order of convergence as for
  Milstein–Galerkin finite element methods but without imposing any
  differentiability condition on the nonlinearity. The results are 
  extended to also incorporate a spectral
  approximation of the driving Wiener process. An application to a stochastic
  partial differential equation is discussed and illustrated through a
  numerical experiment. 


[
    [
    December 30, 2023
=====================





§ INTRODUCTION



In this paper we investigate the numerical solution of non-autonomous
semilinear stochastic evolution equations (SEE) driven by an additive Wiener 
noise. More precisely, let (H, (·, ·), ·) and (U,
(·, ·)_U, ·_U) be two separable -Hilbert spaces. For a given
T ∈ (0,∞) we denote by (Ω_W, ^W, (^W_t)_t ∈ [0,T],
_W) a filtered probability space satisfying the usual conditions. By
(W(t))_t ∈ [0,T] we denote a cylindrical (^W_t)_t ∈ [0,T]-Wiener
process on U with associated covariance operator Q ∈ℒ(U).

Our goal is the approximation of the mild solution to SEEs of the form

    X(t) + [ A X(t) + f(t,X(t)) ] t
        = g(t) W(t),     for  t ∈ (0,T],
    
        X(0) = X_0.

Hereby, we assume that -A (A) ⊂ H → H is the
infinitesimal generator of an analytic semigroup (S(t))_t ∈ [0,∞)⊂ℒ(H) on H. The initial value X_0 Ω_W → H
is assumed to be a p-fold integrable random variable for some p ∈
[2,∞), while the mapping g 
[0,T] →ℒ_2^0 is Hölder continuous with exponent
1/2. 
Here ℒ_2^0 =ℒ_2(Q^1/2(U),H) denotes the set of
all Hilbert–Schmidt operators from Q^1/2(U) to H. 
In addition, the mapping f is
assumed to be Lipschitz continuous. 
A complete and more precise statement of all conditions on
A, f, g, and X_0 is given in Section <ref>.

Under these assumptions the mild solution X  [0,T] ×Ω_W → H to (<ref>) is uniquely determined by the
variation-of-constants formula 

    X(t) = S(t) X_0 - ∫_0^t S(t - s) f(s,X(s)) s + ∫_0^t S(t-s) g(s)
      W(s)

which holds _W-almost surely for all t ∈ [0,T]. Our assumptions in
Section <ref> ensure the existence of a unique mild
solution to (<ref>). We refer the reader to <cit.> for a
general introduction of the semigroup approach to stochastic evolution
equations. We also refer to <cit.> for further details on
cylindrical Wiener processes and Hilbert space valued stochastic integrals. 

Due to the presence of the noise the mild solution 
is, in general, only of very low spatial and temporal regularity. 
This in turn results in low convergence rates of numerical approximations.
Examples for standard numerical methods for SEEs are found, for instance, in
the monographs <cit.> and the references
therein. Because of this, an accurate simulation of stochastic evolution
equations is often computationally expensive. 
This explains why the development of strategies to reduce 
the computational complexity has attracted a lot of attention over the last
decade. In particular, we mention the multilevel Monte-Carlo method that has
been applied to stochastic partial differential equations, for instance,
in <cit.>. However, the success of this approach depends on the
availability of efficient numerical methods which converge with a high order
with respect to the mean-square norm.  

One way to construct such higher order numerical approximations 
is based on Itō–Taylor expansions as discussed in <cit.>.
In fact, provided the coefficient functions are sufficiently smooth, numerical
methods of basically any temporal order can be constructed. However, these
methods sometimes behave unstable in numerical simulations and the necessity to
evaluate higher order derivatives or to generate multiple iterated stochastic
integrals limits their practical relevance. More severely, already 
the imposed smoothness requirements are too restrictive in most applications of
SEEs in infinite dimensions. For 
instance, the general assumption  
in <cit.> asks for the semilinearity f  [0,T]
× H → H to be
infinitely often Fréchet differentiable with bounded derivatives. This
condition already excludes any truly nonlinear Nemytskii-type
operator. Compare further with Remark <ref> below.
We also refer to <cit.> for further numerical methods with a higher order
temporal convergence rate, such as Milstein-type schemes or Wagner–Platen-type
methods. Although the smoothness conditions on f
are substantially relaxed in some of these papers, all results at least 
require the Fréchet differentiability.

The purpose of this paper is the introduction of 
a novel numerical method for the approximation of the solution to
(<ref>) that combines the drift-randomization 
technique from <cit.> for the numerical solution of
stochastic ordinary differential equations (SODEs) with a 
Galerkin finite element method from <cit.>. 
As in <cit.>, it turns out that the resulting method converges
with a higher rate with respect to the temporal discretization parameter
without requiring any (Fréchet-) differentiability of the nonlinearity. 
Our approach also relaxes the smoothness requirements
of the coefficients f and g with respect to the time variable t
considerably.

To introduce the new method more precisely, let k ∈ (0,T) denote
an equidistant temporal step size with associated grid points t_n = n k,
n ∈{1,…,N_k}. Hereby, N_k ∈ is determined by t_N_k = N_k k ≤ T
< (N_k + 1) k. In addition, let (V_h)_h ∈ (0,1)⊂ H be a suitable
family of finite dimensional subspaces, where the parameter h ∈ (0,1)
controls the granularity of the space V_h such that lim_h → 0dist(u, V_h) = 0 for all u ∈ H. The operators P_h  H →
V_h denote the orthogonal projectors onto V_h while A_h  V_h → V_h
is a suitable discrete version of the infinitesimal generator A.
For further details on the spatial discretization we refer to 
Section <ref>. 

For every k ∈ (0,T) and h ∈ (0,1) the
proposed randomized Galerkin finite element method is then given by
the recursion

    X_k,h^n, τ + τ_n k [ A_h X_k,h^n,τ + P_h f(t_n-1,
        X_k,h^n-1) ]    = X_k,h^n-1 + P_h g(t_n-1)
        Δ_τ_n k W(t_n-1),
    
        X_k,h^n + k [ A_h X_k,h^n + P_h f(t_n^τ
        , X_k,h^n, τ) ]    = X_k,h^n-1 + P_h g(t_n^τ)
        Δ_k W(t_n-1)

for all n ∈{1,…,N_k} with initial value X_k,h^0 = P_h X_0.
Hereby, τ = (τ_n)_n ∈ denotes an independent family of
𝒰(0,1)-distributed random variables defined on a further
probability space (Ω_τ, ^τ, _τ). The intermediate 
value of X_k,h^n,τ represents an approximation of X at the random
time point t_n^τ = t_n-1 + τ_n k. Observe that (t_n^τ)_n ∈{1,…,N_k} is an independent family of random variables with t_n^τ∼𝒰(t_n-1,t_n). Further, for all t ∈ [0,T) and κ∈ (0,T - t) we denote the Wiener increments by

    Δ_κ W(t) := W(t + κ) - W(t).


The method (<ref>) constitutes a two-staged 
Runge–Kutta method whose second stage has a randomized node. Compare further
with <cit.>. Further randomized numerical methods for
partial differential equations are studied in <cit.>. Moreover, in case τ_n ≡ 0 for all
n ∈{1,…,N_k} we would recover the 
linearly-implicit Euler–Galerkin finite element method studied, for instance,
in <cit.>. However, the presence of
the artificially randomized internal step
X_h,k^n,τ allows us to prove a higher order
temporal convergence rate compared to standard results in the literature.

In fact, under the assumptions of Section <ref>, the mild solution
(<ref>) to the SEE (<ref>) enjoys the regularity
X ∈ C([0,T];L^p(Ω_W;Ḣ^1+r)) ∩
C^1/2([0,T];L^p(Ω_W;H)), where p ∈ [2,∞) is determined
by the corresponding integrability of the initial value X_0 and the spaces
Ḣ^r = ( (-A)^r/2), r ∈ [0,1), of fractional powers of
the generator -A measure the spatial regularity of X. Then,
according to Theorem <ref> below, there exists C ∈ (0,∞)
such that for every h ∈ (0,1) and k ∈ (0,T) we have

    max_n ∈{1,…,N_k} X(t_n) - X_k,h^n _L^p(Ω;H)≤ C (h^1 + r + k^1/2 + min( r/2, γ)),

where the value of the parameter γ∈ (0,1/2] is determined by
the regularity of f and g with respect to the time variable t.

Note that the standard error estimate for Euler–Maruyama-type methods is only
of order 𝒪(h^1+r + k^1/2) under the same regularity
conditions. The same temporal convergence rate as in
(<ref>) is only recovered for SEEs with additive
noise if the linearly-implicit Euler–Galerkin finite element method is
treated as a Milstein-type scheme, see <cit.>.
This is possible since Milstein-type
schemes coincide with Euler–Maruyama-type methods in this case.
However, as already mentioned above, the error analysis of 
Milstein–Galerkin finite element methods typically requires the
differentiability of the nonlinearity f which is not required for the method
(<ref>).  

The remainder of this paper is organized as follows. After collecting some
notation and auxiliary results from stochastic analysis 
in Section <ref>, we give a more precise statement of all
assumptions imposed on the SEE (<ref>) 
in Section <ref>. In addition, we
also state the main result (<ref>) in this section, see
Theorem <ref>. For the proof of this
error estimate we apply the same methodology as in <cit.>. To this
end, we show in Section <ref> that the method (<ref>) is
bistable. The notion of bistability admits a two-sided estimate of the
error (<ref>) in terms of the local truncation error measured with
respect to a stochastic version of Spijker's norm. This local error is then
estimated in Section <ref>. In Section <ref> we
incorporate an approximation of the Wiener noise
into the method (<ref>).
In Section <ref> we finally apply the method
(<ref>) for the numerical solution of a more explicit stochastic
partial differential equation. 








§ NOTATION AND PRELIMINARIES



In this section we explain the notation used throughout this paper and collect
some auxiliary results from stochastic analysis.

First, we denote by  the set of all positive integers, while _0 = {0 }∪. 
Moreover, let (E_i, ·_E_i), i ∈{1,2},
be two normed -vector spaces. Then,
we denote by ℒ(E_1, E_2) the set of all bounded linear operators
mapping from E_1 to E_2 endowed with the usual operator norm. If E_1 =
E_2 we write ℒ(E_1) = ℒ(E_1,E_1).
If E_i, i ∈{1,2}, are separable Hilbert spaces, then we denote by
ℒ_2(E_1,E_2) the set of all Hilbert–Schmidt operators mapping from
E_1 to E_2. Recall that the Hilbert–Schmidt norm of L ∈ℒ_2(E_1,E_2) is given by

    L _ℒ_2(E_1,E_2) = ( ∑_j = 1^∞ L
      φ_j_E_2^2 )^1/2,

where (φ_j)_j ∈⊂ E_1 is an arbitrary orthonormal basis.
As mentioned in the introduction we use 
the short hand notation ℒ_2^0 =ℒ_2(Q^1/2(U),H)
and ℒ_2:=ℒ_2(U,H). 

Let us also recall a few function spaces which play an important role
in this paper. As usual, we denote by L^p(0,T;E), p ∈ [1,∞), the 
space of all p-fold integrable mappings v  [0,T] → E with
values in a Banach space (E, ·_E) endowed with the standard norm. When E=ℝ, we write L^p(0,T) for short.

We mostly measure temporal regularity of the exact solution and the coefficient
functions in terms of Hölder continuity, that is with respect to the norm

    v _C^α([0,T];E) = sup_t ∈ [0,T] v(t) _E +
      sup_0ptt_1, t_2 ∈ [0,T]t_1 ≠ t_2 v(t_1) - v(t_2) _E/|t_1 - t_2|^α,

where α∈ (0,1] denotes the Hölder exponent. The space of all
α-Hölder continuous mappings taking values in E
is denoted by C^α([0,T];E).

For the same purpose we also make use of the family
W^ν,p(0,T;E) ⊂ L^p(0,T;E) of (fractional) Sobolev 
spaces. Recall that for p ∈ [1,∞) and ν = 1 the Sobolev space
W^1,p(0,T;E) is endowed with the norm

    v _W^1,p(0,T;E) 
      = ( ∫_0^T  v(t) ^p_E t + ∫_0^T  v'(t)
      ^p_E t)^1/p,

where v' ∈ L^p(0,T;E) denotes the weak derivative of 
v ∈ W^1,p(0,T;E). Moreover, for p ∈ [1,∞) and ν∈ (0,1)
the Sobolev–Slobodeckij norm is given by

    v _W^ν,p(0,T;E) 
      = ( ∫_0^T  v(t) ^p_E t + ∫_0^T ∫_0^T v(t_1) -
      v(t_2) ^p_E/|t_1 - t_2|^1 + ν pt_2t_1)^1/p.

Further details on fractional Sobolev spaces are found, for
example, in <cit.> and <cit.>.


Our numerical method (<ref>) yields a discrete-time stochastic
process defined on the product probability space 

    (Ω, , ) := (Ω_W×Ω_τ, ^W⊗^τ,
      _W ⊗_τ),

where the corresponding expectations are denoted by 𝔼_W and
𝔼_τ. The additional random input τ = (τ_n)_n ∈ 
in (<ref>) induces a natural filtration (_n^τ)_n ∈_0 on
(Ω_τ, 
^τ, _τ) by setting _0^τ := {∅, Ω_τ} and
_n^τ := σ{τ_j   :   1 ≤ j ≤ n } for n ∈.

Moreover, for each k ∈ (0,T) let 

    π_k := { t_n = nk   :   n = 0,1,…,N_k}⊂ [0,T]

be the set of temporal grid points with equidistant step size k.
Hereby, N_k ∈ is uniquely determined by t_N_k = N_k k ≤ T <
(N_k + 1)k. For each temporal grid π_k a 
discrete-time filtration (^π_k_n)_n ∈{0,…,N_k}
on (Ω, , ) is given by 

    ^π_k_n := ^W_t_n⊗^τ_n,    for  n ∈{0,1,…,N_k}.


As a useful estimate for higher moments of stochastic integrals, a particular
version of a Burkholder–Davis–Gundy-type inequality is presented here for
later use. The proposition follows directly from <cit.>.


   
  For every p∈ [2,∞) there exists a constant C_p ∈ [0,∞) 
  such that for all s,t ∈ [0,T], s < t, and for all
  (_t^W)_t ∈ [0,T]-predictable stochastic processes 
  Y  [0,T] ×Ω_W →ℒ_2^0 satisfying  
  
    ( ∫_s^tY(ξ)^2_L^p(Ω_W;ℒ_2^0)ξ)^1/2 < ∞,

  we have
  
    ∫_s^t Y(ξ) W(ξ)_L^p(Ω_W;H)≤ C_p ( ∫_s^tY(ξ)^2_L^p(Ω_W;ℒ_2^0)ξ)^1/2.



A further important tool is the following Burkholder–Davis–Gundy-type
inequality for discrete-time martingales with values in 
a Hilbert space. For a proof we refer to <cit.>:


   
  For every p∈ [1,∞) there exist constants 
  c_p, C_p ∈ [0,∞) such that for every discrete-time H-valued
  martingale (Y^n)_n∈ℕ_0 and for every n ∈_0 we have  
  
    c_p  [Y]_n^1/2_L^p(Ω)≤max_j∈{0,…,n}Y^j_L^p(Ω;H)≤ C_p  [Y]_n^1/2_L^p(Ω),

  where [Y]_n = Y^0^2 + ∑_k=1^nY^k-Y^k-1^2
  is the quadratic variation of (Y^n)_n ∈_0.








§ ASSUMPTIONS AND MAIN RESULT



In this section we collect all essential conditions on the stochastic evolution
equation (<ref>). Then the main result is stated.


  
  The linear operator A (A) ⊂ H → H is densely defined,
  self-adjoint, and positive definite with compact inverse.



Assumption <ref> implies the existence of a positive, increasing
sequence (λ_i)_i∈ℕ⊂ such that 
0<λ_1 ≤λ_2 ≤… with lim_i→∞λ_i = ∞, and of an orthonormal basis (e_i)_i∈ℕ of H such that A e_i = λ_i  e_i for every i ∈ℕ.

In addition, it also follows from Assumption <ref> that -A is the
infinitesimal generator of an analytic semigroup (S(t))_t ∈ [0,∞)⊂ℒ(H) of contractions. More precisely, the family
(S(t))_t ∈ [0,∞) enjoys the properties

    S(0)    = Id∈ℒ(H),
    
      S(s + t)    = S(s) ∘ S(t) = S(t) ∘ S(s), 
        for all  s,t ∈ [0,∞),

and

    sup_t ∈ [0,∞) S(t) _ℒ(H)≤ 1.

A more detailed account on the theory of linear semigroups is found in
<cit.>. 

Further, let us introduce fractional powers of A, which are used to measure
the (spatial) regularity of the mild solution (<ref>). For any r≥
0 we define the operator A^r/2(A^r/2) =
{x∈ H   :  ∑_j=1^∞λ_j^r (x,e_j)^2 < ∞}⊂ H → H by 

    A^r/2 x
      := ∑_j=1^∞λ_j^r/2 (x,e_j) e_j,
        for all  x ∈(A^r/2).

Then, by setting (Ḣ^r,(·,·)_r, ·_r) :=
((A^r/2), (A^r/2·, A^r/2·),
A^r/2·), r ∈ [0,∞), we obtain a family of
separable Hilbert spaces.  



  The assumption on A can be relaxed such that A is not necessarily 
  self-adjoint. 
  In that case the fractional powers of A and the
  spaces Ḣ^r can be defined in a different way. For instance,
  we refer to <cit.>. 
  For the validity of our main result Theorem <ref> it is then crucial
  to find a suitable replacement for the assertions of
  Lemma <ref> whose proof depends on the self-adjointness of A.
  For example, we refer to <cit.> for such error estimates in the non-self-adjoint case.
  Compare further with <cit.> for an SPDE related result.


After these preparations we are able to state the assumptions on the initial
condition X_0 as well as on the drift and diffusion coefficient functions.


  
  There exist p ∈ [2,∞) and r ∈ [0,1] such that the initial value
  X_0 Ω_W → H satisfies X_0 ∈ L^p(Ω_W, _0^W,_W;
  Ḣ^1+r).



  
  The mapping f  [0,T] × H → H is continuous.
  Moreover, there exist γ∈ (0,1/2] and C_f ∈
  (0,∞) such that 
  
    f(t,u_1) - f(t,u_2)    ≤ C_f u_1 - u_2,
     f(t_1,u) - f(t_2,u)    ≤ C_f (1 + u) |t_1 -
        t_2|^γ

  for all u, u_1, u_2 ∈ H and t, t_1,t_2 ∈ [0,T].


From Assumption <ref> we directly deduce a linear growth bound
of the form

    f(t,u)≤Ĉ_f ( 1 + u ),  for all 
      t ∈ [0,T],  u ∈ H,

where Ĉ_f:=f(0,0)+C_f(1 +  T^γ). Moreover, we emphasize that
the regularity of f with respect to t is even weaker than in
<cit.> for the linearly-implicit Euler–Galerkin
finite element method. We refer to Section <ref> for a class of
mappings satisfying Assumption <ref>.


  
  The mapping g  [0,T] →ℒ_2^0 is continuous.
  Moreover, there exist p ∈ [2,∞),
  r ∈ [0,1], γ∈ (0,1/2], and C_g ∈
  (0,∞) such that
  
    g _C^1/2([0,T];ℒ_2^0)
        +  A^r/2 g _C([0,T];ℒ_2^0) 
        +  g _W^1/2 + γ,p(0,T;ℒ_2^0)   ≤ C_g.



Assumptions <ref> to <ref> with r∈ [0,1) and p∈
[2,∞) are sufficient to ensure the existence of a unique mild solution
X  [0,T] ×Ω_W → H to the stochastic evolution equation
(<ref>) with

    sup_t∈ [0,T]_W[ X(t)_1+r^p] < ∞,

and there exists a constant C depending on r and p such that

    ( _W [ X(t_1)-X(t_2)_r^p] )^1/p≤ C|t_1-t_2|^1/2,

for each t_1, t_2 ∈ [0,T]. For proofs of these regularity results 
we refer, for instance, to <cit.> and
<cit.>.

Next, we formulate the assumption on the spatial discretization.
To this end, let (V_h)_h ∈ (0,1)⊂Ḣ^1 be a family of
finite dimensional subspaces. Then, we introduce the 
Ritz projector R_h Ḣ^1→ V_h
as the orthogonal projector onto V_h with respect to the
inner product (·,·)_1. To be more precise, the Ritz projector is
given by 

    (R_h x, y_h)_1 = (x,y_h)_1   for all  x ∈Ḣ^1, y_h∈
      V_h.

The following assumption is used to quantify the speed of convergence
with respect to the spatial parameter h ∈ (0,1).


  
  Let a sequence (V_h)_h∈ (0,1) of finite dimensional subspaces of
  Ḣ^1 be given such that there exists a constant C ∈ (0,∞) with 
  
    R_h x -x≤ Ch^s x_s    for all 
        x∈Ḣ^s,  s∈{1,2},  h ∈ (0,1).



Similar estimates are obtained for the approximation of x ∈Ḣ^1 + r, r ∈ [0,1], by interpolation. A typical example of a
spatial discretization satisfying Assumption <ref> is the spectral
Galerkin method. This method is obtained by setting h = 1/N, N ∈, and V_h := span{ e_j   :   j = 1,…,N }, where
(e_j)_j ∈ denotes the family of eigenfunctions of A. For more details
see <cit.>. A further example is the standard
finite element method from <cit.> as we will discuss in
Section <ref>.

We are now well-prepared to formulate the main result of this paper.
The proof is deferred to the end of Section <ref>.


  
  Let Assumptions <ref> to <ref> be fulfilled for some p ∈
  [2,∞), r ∈ [0,1), and γ∈ (0,1/2]. Then there
  exists C ∈ (0,∞) such that for every h ∈ (0,1) and
  k ∈ (0,T) 
  
    max_n ∈{0,…,N_k} X(t_n) - X_k,h^n _L^p(Ω;H)
        ≤ C (1 +  X _C([0,T];L^p(Ω_W;Ḣ^1+r))
        +  X _C^1/2([0,T];L^p(Ω_W;H)))
        ( h^1+r + k^1/2 + min(r/2, γ)),

  where X denotes the mild solution (<ref>) to the stochastic
  evolution equation (<ref>) and (X_k,h^n)_n ∈{0,…,N_k}
  denotes the stochastic process generated by the randomized Galerkin finite
  element method (<ref>).



  
  In order to obtain the same temporal order of convergence as in
  Theorem <ref>, other results in the literature
  usually impose additional smoothness conditions on the
  nonlinearity. For instance, in <cit.>
  the authors require f ∈ C^∞_b(H;H), that is, f is infinitely
  often Fréchet differentiable with bounded derivatives. 
  However, this condition is too restrictive for all SPDEs with a truly
  nonlinear Nemytskii operator. In <cit.>, this problem is circumvented by instead 
  requiring f ∈ C^2_b(Ḣ^β;H) or
  f ∈ C^2_b(H;Ḣ^-β) for β∈ (0,1]. 
  Such conditions allow to treat some Nemytskii-type operators. In particular,
  we refer to <cit.>. However, the presence of the
  parameter β often results into a lower temporal 
  convergence rate. For instance, the rate is only equal to 1/2 if
  β = 1 in <cit.>.








§ BISTABILITY



In this section we show that the randomized Galerkin finite element
method (<ref>) constitutes a bistable numerical method in the
sense of <cit.>. More precisely, for each choice of
h ∈ (0,1) and k ∈ (0,T), we first observe that 
the scheme (<ref>) can be written as 
an abstract one-step method of the form 

    X_k,h^n    = S_k,h X_k,h^n-1 + Φ^n_k,h(X_k,h^n-1, τ_n),
           n ∈{1,…,N_k},
    
        X_k,h^0   = ξ_h,

in terms of a suitable family of linear operators S_k,h∈ℒ(H) and associated  
increment functions Φ_k,h^n  H × [0,1]
×Ω_W → H. Then, we 
verify the conditions of a stability theorem from <cit.> that yields
two-sided stability bounds for general one-step methods of the form
(<ref>).

For each k ∈ (0,T) let π_k := { t_n = nk   :   n = 0,1,…,N_k}⊂ [0,T] be the set of temporal grid points with equidistant step size
k as defined in (<ref>).
As in the introduction, we denote by P_h  H → V_h, h ∈
(0,1), the orthogonal projector onto the finite dimensional subspace V_h ⊂Ḣ^1 with respect to the inner product in H.

In this situation, we define ξ_h := P_h X_0 ∈ L^p(Ω_W;H) as the
initial condition for the numerical scheme (<ref>). 
Under Assumption <ref> with p ∈ [2,∞) 
it then holds 

    sup_h ∈ (0,1)ξ_h _L^p(Ω_W;H)≤ X_0 _L^p(Ω_W;H)

due to P_h _ℒ(H) = 1 for all h ∈ (0,1).

Next, for each h ∈ (0,1), we implicitly define a discrete version A_h
 V_h → V_h of the generator A by the relationship

    ( A_h x_h, y_h ) = ( x_h, y_h )_1,    for all  x_h,
      y_h ∈ V_h.

From Assumption <ref> it then follows immediately that A_h is symmetric
and positive definite. Moreover, for each h ∈ (0,1) and k ∈ (0,T) we
obtain a bounded linear operator S_k,h∈ℒ(H) defined by

    S_k,h := ( Id + k A_h )^-1 P_h.

For the error analysis it is also convenient to introduce a piecewise constant
interpolation of S_k,h to the whole time interval, which we denote by
S_k,h [0,T] →ℒ(H):  
For each h ∈ (0,1) and k ∈ (0,T) let S_k,h(t)
:= S_k,h^N_k for all t ∈ [t_N_k,T] and

    S_k,h(t):= ( Id + k A_h )^-j P_h,  t∈ [t_j-1,t_j)

for j∈{1,2,…,N_k}.  The following lemma contains some useful
stability bounds for S_k,h and S_k,h uniformly with respect
to the discretization parameters h and k. For a proof and more general
versions of these estimates we refer to <cit.>. 


  
  Let Assumption <ref> be satisfied. Then, the operator S_k,h given in (<ref>) is well-defined for all h ∈ (0,1) and k ∈ (0,T). 
  Furthermore, we have 
  
    sup_k∈ (0,T)sup_h∈ (0,1)S_k,h_ℒ(H)≤ 1.

  In addition, the continuous-time interpolation S_k,h [0,T]
  →ℒ(H) of S_k,h is right-continuous with existing
  left-limits. It holds true that
  
    sup_k∈ (0,T)sup_h∈ (0,1)sup_t∈ [0,T]S_k,h(t) _ℒ(H)≤ 1.



Now we are in a position to introduce the increment functions associated 
to the numerical method (<ref>). For each k ∈ (0,T), h ∈
(0,1), and j ∈{1,…,N_k} we define Φ_k,h^j  H ×
[0,1] ×Ω_W → H and Ψ_k,h^j  H × [0,1] ×Ω_W → H by setting

    Φ_k,h^j(x,τ) := -k S_k,h f(t_j-1 + τ k,
      Ψ^j_k,h(x,τ)) 
      + S_k,h g(t_j-1 + τ k) Δ_k W(t_j-1)

and

    Ψ_k,h^j(x,τ) := S_τ k,h x - τ k S_τ k,h f(t_j-1, x ) 
      + S_τ k,h g(t_j-1) Δ_τ k W(t_j-1)

for all x ∈ H and τ∈ [0,1]. We refer to (<ref>) for the
definition of the Wiener increments Δ_k W(t).

Observe that, under the assumptions of Section <ref>,
for each τ∈ [0,1] and j ∈{1,…,N_k} the mapping (x,ω)
↦Φ_k,h^j(x,τ)(ω) is measurable with respect to
ℬ(H) ⊗_t_j^W/ ℬ(H). Moreover, for each x ∈
H and almost all ω∈Ω we have that the mapping [0,1] ∋τ↦Φ^j(x,τ)(ω) ∈ H is continuous due to the pathwise
continuity of the Wiener process and the continuity of f and g.

Altogether, this shows that the numerical method (<ref>) can
be rewritten as a one-step method of the form (<ref>). The family of
random variables (X_k,h^n)_n ∈{0,…,N_k}, which is determined by
(<ref>), is therefore a discrete-time stochastic process on the
product probability space (Ω, , ) defined in (<ref>).
Moreover, it is adapted to the filtration 
(^π_k_n)_n ∈{0,…,N_k} from (<ref>).

Let us now recall the notion of bistability from <cit.>.
For this, we first introduce a family of linear spaces consisting of all 
(^π_k_n)_n ∈{0,…,N_k}-adapted and p-fold integrable grid
functions on π_k. To be more precise, for p ∈ [2,∞) and k ∈
(0,T) we set 𝒢^p_k := 𝒢^p_k (π_k,
L^p(Ω;H)) with

    𝒢^p_k := { (Z_k^n)_n ∈{0,…,N_k}  :   
      Z_k^n ∈ L^p(Ω, _n^π_k, ; H)  for all  
      n ∈{0,1,…,N_k}}.

In addition, we endow the spaces 𝒢_k^p with two different norms.
For arbitrary Z_k = (Z_k^n)_n ∈{0,…,N_k}∈𝒢^p_k
these norms are given by

    Z_k _∞,p := max_n ∈{0,…,N_k} Z_k^n
      _L^p(Ω;H)

and, for each h ∈ (0,1),

    Z_k _S,p,h :=  Z^0_k _L^p(Ω;H) + max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j Z_k^j_L^p(Ω;H),

where S_k,h has been defined in (<ref>). The norm ·_S,p,h is called (stochastic) Spijker norm. Deterministic versions 
of this norm are used in numerical analysis for finite
difference methods, for instance, in <cit.>. In <cit.> a more detailed discussion is
given in the context of stochastic differential equations. 

The final ingredient for the introduction of the 
stability concept is then the following family of
residual operators ℛ_k,h𝒢^p_k →𝒢_k^p associated to the numerical scheme (<ref>). For each
p ∈ [2,∞), k ∈ (0,T), and h ∈ (0,1) the residual of an
arbitrary grid function Z_k ∈𝒢_k^p is given by

    ℛ_k,h[Z_k](t_0) := Z_k^0 - ξ_h,   
    ℛ_k,h[Z_k](t_n) := Z_k^n - S_k,h Z_k^n-1 -
        Φ_k,h^n(Z_k^n-1,τ_n),      n ∈{1,…,N_k}.

It is not immediately evident if the residual operators are actually
well-defined for every given h ∈ (0,1) and k ∈ (0,T).
From Theorem <ref> it follows that indeed
ℛ_k,h[Z_k] ∈𝒢_k^p for all Z_k ∈𝒢_k^p under Assumptions <ref> to <ref>. 

The following definition of bistability is taken from 
<cit.>.


  The numerical scheme (<ref>) is called bistable with respect
  to the norms ·_∞,p and ·_S,p,h if there
  exists p ∈ [2,∞) such that the residual operators
  ℛ_k,h𝒢_k^p →𝒢_k^p are
  well-defined and bijective for all k ∈ (0,T) and h ∈ (0,1). In
  addition, there exists C_Stab∈ (0,∞) such that
  for all k ∈ (0,T), h ∈ (0,1), and Y_k, Z_k ∈𝒢^p_k we
  have
  
    1/C_Stabℛ_k,h[Y_k] - ℛ_k,h[Z_k] _S,p,h   ≤ Y_k - Z_k _∞,p
       ≤ C_Stabℛ_k,h[Y_k] - ℛ_k,h[Z_k] _S,p,h.
  


Under the assumptions of Section <ref> the following lemma shows
that the family of increment functions Φ_k,h is bounded at 0 ∈ H and
in a certain sense Lipschitz continuous, uniformly with respect to the
discretization parameters h ∈ (0,1) and k ∈ (0,T). These estimates are
required for the stability theorem further below.


  
  Under Assumptions <ref> to <ref> 
  there exist C_Φ,0, C_Φ,1∈ (0,∞) with
  
    sup_h ∈ (0,1)sup_k ∈ (0,T)∑_j = m^n S^n-j_k,hΦ_k,h^j(0,τ_j)
        _L^p(Ω;H)≤ C_Φ,0 (t_n - t_m-1)^1/2

  for all n,m ∈{1,…,N_k} with m ≤ n. Moreover, 
  it holds true that
  
    sup_h ∈ (0,1)∑_j = 1^n S_k,h^n-j( Φ_k,h^j(Y_k^j-1, τ_j) 
          - Φ_k,h^j( Z_k^j-1, τ_j) ) _L^p(Ω;H)
         ≤ C_Φ,1 k ∑_j = 1^n  
           Y_k^j-1 - Z_k^j-1_L^p(Ω;H)

  for all k ∈ (0,T), Y_k, Z_k ∈𝒢_k^p, and
  n ∈{1,…,N_k}.



  We first verify (<ref>). From (<ref>) we obtain
  
    ∑_j = m^n S^n-j_k,hΦ_k,h^j(0,τ_j)
        _L^p(Ω;H)≤∑_j = m^n k
        S^n-j+1_k,hf(t^τ_j, Ψ_k,h^j(0,τ_j))
        _L^p(Ω;H)
         +∑_j = m^n S^n-j+1_k,hg(t^τ_j)Δ_k W(t_j-1)
        _L^p(Ω;H)=:I_1+I_2,

  where t^τ_j:=t_j-1+τ_j k.

  For the estimate of I_1, we first apply the triangle inequality and
  (<ref>). Then, applying the linear growth (<ref>) of
  f and the boundedness of S_k,h in Lemma <ref>
  yields
  
    I_1
           ≤∑_j=m^n k 
        S^n-j+1_k,h f( t^τ_j, Ψ_k,h^j(0,τ_j))
        _L^p(Ω;H)
       ≤ k ∑_j=m^n 
        f(t^τ_j, Ψ_k,h^j(0,τ_j)) _L^p(Ω;H)
       ≤Ĉ_f k ∑_j=m^n( 1 + 
        Ψ_k,h^j(0,τ_j) _L^p(Ω;H))
       ≤Ĉ_f (1 + sup_j ∈{m, …,n}Ψ_k,h^j(0,τ_j)_L^p(Ω;H))
        (t_n-t_m-1).

  Thus, for the estimate of I_1 it remains to show that
  Ψ_k,h^j(0,τ_j)_L^p(Ω;H)
  can be bounded uniformly. Indeed, by definition of Ψ_k,h
  in (<ref>), the 
  linear growth of f in (<ref>), Assumption <ref>,
  Proposition <ref>, and estimate (<ref>)
  we have for each j  
  
    Ψ_k,h^j(0,τ_j)_L^p(Ω;H)
         ≤τ_j k S_τ_j k, hf(t_j-1,0)_L^p(Ω;H)+
        S_τ_j k, hg(t_j-1)Δ_τ_j k W(t_j-1)_L^p(Ω;H)
         ≤( Ĉ_f T + C_p C_g T^1/2),

  which is independent of j, h, and k. 

  For the estimate of I_2 we first define a new process
  ĝ [0,T] ×Ω_τ→ℒ_2^0 by
  
    ĝ(t):=g(t_j-1 + τ_j k),  t∈ [t_j-1,t_j).

  Then, we rewrite the sum as a stochastic integral by inserting
  (<ref>) and replacing g by ĝ. 
  An application of Proposition <ref>, estimate (<ref>), and
  Assumption <ref> yields 
  
    I_2 
             = ∫_t_m-1^t_nS_k,h(t_n-r) ĝ(r)
          W(r)_L^p(Ω;H)
       =(_τ[
          _W[∫_t_m-1^t_nS_k,h(t_n-r) ĝ(r) W(r)^p] ])^1/p
       ≤ C_p (_τ[ (∫_t_m-1^t_nS_k,h(t_n-r) ĝ(r) 
          ^2_ℒ_2^0r)^p/2])^1/p
       ≤ C_p ( ∫_t_m-1^t_nsup_s ∈ [0,T]g(s)^2_ℒ_2^0r)^1/2≤ C_p C_g (t_n-t_m-1)^1/2.

  Therefore, we obtain (<ref>) with
  
    C_Φ,0 :=
        Ĉ_f T^1/2(1 + Ĉ_f T + C_p C_g T^1/2)
        + C_p C_g.

  
  It remains to verify (<ref>). For this let Y_k, Z_k ∈𝒢^p_k be arbitrary. Then, by inserting the
  definition of Φ_k,h from (<ref>) and an application of
  Assumption <ref> we get
  
    ∑_j = 1^n S_k,h^n-j( Φ_k,h^j(Y_k^j-1, τ_j) 
        - Φ_k,h^j( Z_k^j-1, τ_j) ) _L^p(Ω;H)
          = k ∑_j = 1^n S_k,h^n-j+1( f (t^τ_j,Ψ_k,h^j(Y_k^j-1, τ_j) )
        - f(t^τ_j, Ψ_k,h^j(Z_k^j-1, τ_j) )) 
        _L^p(Ω;H)
         ≤  k C_f ∑_j = 1^n Ψ_k,h^j(Y_k^j-1, τ_j) - Ψ_k,h^j(Z_k^j-1, τ_j) 
        _L^p(Ω;H).

  In addition, from the same arguments we also deduce the bound
  
    Ψ_k,h^j(Y_k^j-1, τ_j)
        - Ψ_k,h^j(Z_k^j-1, τ_j) _L^p(Ω;H)
         ≤S_τ k,h (Y_k^j-1-Z_k^j-1)
        + τ_j k S_τ k,h ( f(t_j-1,Y_k^j-1) 
        - f(t_j-1,Z_k^j-1)) _L^p(Ω;H)
         ≤ (1+kC_f)Y_k^j-1-Z_k^j-1_L^p(Ω;H).

  Altogether, this proves (<ref>) with
  
    C_Φ,1 :=  C_f (1 + T C_f).

  This completes the proof of the lemma.


Next, we observe that (<ref>), Lemma <ref>, and
Lemma <ref> verify together all conditions of the stability theorem 
<cit.>. Therefore, we immediately obtain
the main result of this section:


  
  Let Assumptions <ref> to <ref> be satisfied with p ∈ [2,∞).
  Then, the randomized Galerkin finite element method (<ref>) is
  bistable with respect to the norms ·_∞,p and ·_S,p,h. 









§ CONSISTENCY AND CONVERGENCE



In the previous section it was proven that the randomized Galerkin finite
element method (<ref>) is bistable. In this section we complete the
error analysis by first deriving estimates for the local truncation error
of the mild solution to the
stochastic evolution equation (<ref>). Together with the stability
inequality (<ref>) these estimates then also yield estimates 
for the global discretization error with respect to the norm in 
L^p(Ω;H).

Let X  [0,T] ×Ω→ H denote the mild solution 
(<ref>) to the stochastic evolution equation (<ref>). For an
arbitrary step size k ∈ (0,T), we transform the stochastic process X into
a grid function by restricting it to the grid points in π_k. More formally,
we obtain X|_π_kπ_k → L^p(Ω;H) by defining

    X|_π_k(t_n) = X(t_n)

for all n ∈{0,1,…,N_k}.
From (<ref>) it follows that indeed 
X|_π_k∈𝒢^p_k for each k ∈ (0,T). Hence, we can apply
the residual operator ℛ_k,h from (<ref>) to
X|_π_k. The local truncation error is then given by

    ℛ_k,h[ X|_π_k ] _S,p,h.

In order to derive an estimate of the local truncation error we first 
recall the definition of the stochastic Spijker norm from (<ref>).
Then we insert the variation-of-constants formula (<ref>) and the definition of
the residual operator (<ref>). After some elementary 
rearrangements we arrive at the inequality

    ℛ_k,h[ X|_π_k ] _S,p,h
         ≤ X_0 - ξ_h _L^p(Ω_W;H)
        + max_n ∈{1,…,N_k} ( S(t_n) - S_k,h^n) X_0
        _L^p(Ω_W;H)
            + max_n ∈{1,…,N_k}∑_j = 1^n ∫_t_j-1^t_j( S(t_n -s) - S_k,h^n-j+1) f(s,X(s)) s_L^p(Ω_W;H)
            + max_n ∈{1,…,N_k}∑_j = 1^n ∫_t_j-1^t_j( S(t_n -s) - S_k,h^n-j+1) g(s) W(s)_L^p(Ω_W;H)
            + max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j( - ∫_t_j-1^t_j S_k,h f(s,X(s)) s
                          
        + ∫_t_j-1^t_j S_k,h g(s) W(s) 
        - Φ_k,h^j( X(t_j-1), τ_j) ) _L^p(Ω;H),

where the linear operators S_k,h∈ℒ(H) and the associated 
increment functions Φ_k,h, k ∈ (0,T), h ∈ (0,1),
are defined in (<ref>) and (<ref>), respectively. For a more
detailed proof of (<ref>) we refer to
<cit.>.

The following sequence of lemmas contains some bounds for the terms on the
right hand side of (<ref>). First, we are concerned with the
consistency of the initial value of the numerical scheme.

 
  
  Let Assumption <ref> and Assumption <ref> be
  satisfied with p ∈ [2,∞) and r ∈ [0,1]. Then there exist C ∈
  (0,∞) such that
  
    X_0 - ξ_h _L^p(Ω_W;H)≤ C h^r+1 X_0 _L^p(Ω_W;Ḣ^1+r)   for all 
        h ∈ (0,1).




  After inserting ξ_h = P_h X_0 we obtain
  
    X_0-ξ_h_L^p(Ω_W;H)   =  (Id - P_h) X_0 _L^p(Ω_W;H)≤(Id - R_h) X_0 _L^p(Ω_W;H)
       ≤ C h^r+1 X_0 _L^p(Ω_W;Ḣ^1+r),

  where the first inequality follows
  from the best approximation property of the orthogonal projector P_h 
  H → V_h, while the last line is due to 
  Assumption <ref> and Assumption <ref>. 


Next, we collect some well-known error estimates for the approximation of the
semigroup (S(t))_t ∈ [0,T]⊂ℒ(H). Recall the definition
of S_k,h from (<ref>).
For a proof of the first two error bounds in Lemma <ref>
we refer to <cit.>. A proof for (<ref>) and
(<ref>) is found in <cit.>.


  
  Let Assumptions <ref> and <ref> be satisfied.
  Then, for every ρ∈ [0,2] there exists C ∈ (0,∞)
  such that for all k ∈ (0,T), h ∈ (0,1), t∈ (0,T]
  we have
  
    ( S(t) - S_k,h(t) ) x ≤ C ( h^ρ+k^ρ/2) x_ρ   for all  x ∈Ḣ^ρ,

  and
  
    ( S(t) - S_k,h(t) ) x ≤ C ( h^ρ + k^ρ/2)
        t^-ρ/2x   for all  x ∈ H.

  In addition, there exists C ∈ (0,∞) such that for all 
  t ∈ [0,T], h ∈ (0,1), k ∈ (0,T), and x ∈ H we have
  
    ∫_0^t ( S(s) - S_k,h(s) ) x
        s≤ C ( h^2 + k ) x .

  Moreover, for every r ∈ [0,1] there exists C ∈ (0,∞) with
  
    ( ∫_0^t ( S(s) - S_k,h(s) ) 
        x ^2 s)^1/2≤ C (h^1 + r+ k^1+r/2) x_r

  for all k ∈ (0,T), h ∈ (0,1), x ∈Ḣ^r, and t ∈ (0,T].


By Lemma <ref> we can directly estimate several of the terms on the right hand side of (<ref>). We begin with the
error with respect to the initial condition.


  Let Assumption <ref> be satisfied with p ∈ [2,∞) and r ∈
  [0,1]. Then, it holds true that
  
    max_n ∈{1,…,N_k} ( S(t_n) - S_k,h^n) X_0
        _L^p(Ω_W;H)≤ C (h^1+r+k^1+r/2) X_0_L^p(Ω_W;Ḣ^1+r)

  for all h∈ (0,1) and k∈ (0,T).



  The assertion follows directly from 
  Assumption <ref> and the corresponding discrete-time version of
  (<ref>).



  
  Let Assumption <ref> to Assumption <ref> be fulfilled with p ∈
  [2,∞) and γ∈ (0,1/2]. Then there exists C ∈
  (0,∞) such that for all h ∈ (0,1) and k ∈ (0,T)
  we have
  
    max_n ∈{1, …, N_k}∑_j = 1^n ∫_t_j-1^t_j( S(t_n -s) - S_k,h^n-j+1) f(s,X(s)) s_L^p(Ω_W;H)
           ≤ C (1 +  
         X _C^1/2([0,T];L^p(Ω_W;H))) ( h^2 + k ).
 



  First, we replace S_k,h^n-j+1 by its piecewise constant interpolation
  S_k,h defined in (<ref>). 
  After adding and subtracting a few additional terms we arrive at
  
    ∑_j = 1^n ∫_t_j-1^t_j( S(t_n -s) - S_k,h^n-j+1) f(s,X(s)) s_L^p(Ω_W;H)
         ≤∫_0^t_n( S(t_n -s) - S_k,h(t_n-s) )
        ( f(s,X(s)) - f(s,X(t_n)) ) s_L^p(Ω_W;H)
            + ∫_0^t_n( S(t_n -s) - S_k,h(t_n-s) ) 
        ( f(s,X(t_n)) - f(t_n,X(t_n)) ) s_L^p(Ω_W;H)
            + ∫_0^t_n( S(t_n -s) - S_k,h(t_n-s) )
        f(t_n,X(t_n)) s_L^p(Ω_W;H)
          =: J_1^n + J_2^n +J_3^n

  for all n ∈{1,…,N_k}. We estimate the three terms separately.
  For J_1^n, we apply estimate (<ref>) with ρ = 2,
  Assumption <ref>, and the Hölder continuity (<ref>) of
  the exact solution. This yields  
  
    J_1^n 
           ≤  C C_f ( h^2 + k )
        ∫_0^t_n ( t_n - s )^-1 X(s)-X(t_n) _L^p(Ω_W;H)s
       ≤  C C_f 
         X _C^1/2([0,T];L^p(Ω_W;H))( h^2 + k )
        ∫_0^t_n ( t_n - s )^-1 + 1/2s
       ≤ C C_f T^1/2 X _C^1/2([0,T];L^p(Ω_W;H))( h^2+k ).

  Similarly, we obtain that
  
    J_2^n
           ≤  C C_f ( h^2+k ) ∫_0^t_n
        ( t_n - s)^-1 + γ(1 +  X(t_n) _L^p(Ω_W;H))
        s
       ≤ C C_f 1/γ  T^γ(1 + 
        X_C([0,T];L^p(Ω_W;H))) ( h^2 + k ).

  Concerning the term J_3^n,
  we apply the estimate (<ref>)
  and the linear growth bound (<ref>) of f. This yields
  
    J_3^n
           ≤ C ( h^2 + k )  f (t_n, X(t_n) ) 
        _L^p(Ω_W;H)
       ≤ C Ĉ_f ( 1 + 
        X _C([0,T];L^p(Ω_W;H)))
        ( h^2 + k ).

  After combining the estimates for J_1^n, J_2^n, J_3^n the proof is
  completed.



  
  Let Assumptions <ref>, <ref>, and <ref> be 
  fulfilled with p ∈ [2,∞) and r ∈ [0,1). Then, there exists C
  ∈ (0,∞) such that
  
    max_n ∈{1,…,N_k}∑_j = 1^n ∫_t_j-1^t_j( S(t_n -s) - S_k,h^n-j+1) g(s) W(s)_L^p(Ω_W;H)
       ≤ C ( A^r/2 g _C([0,T];ℒ_2^0)
        + T^1-r/2 g _C^1/2([0,T];ℒ_2^0))
        (h^1+r+k^1+r/2)
 
  for all h∈ (0,1) and k∈ (0,T).



  As in the proof of Lemma <ref>, we first
  replace the discrete-time operator S_k,h
  by its piecewise constant interpolation S_k,h
  defined in (<ref>).
  This enables us to apply Proposition <ref> for each n∈{1,…,N_k}. After adding and subtracting an additional term, we obtain
  
    ∑_j = 1^n ∫_t_j-1^t_j( S(t_n -s) - S_k,h^n-j+1) g(s) W(s)_L^p(Ω_W;H)
         ≤ C_p (∫_0^t_n(S(t_n -s) - S_k,h(t_n -s) ) g(t_n) 
        ^2_ℒ_2^0s)^1/2
            + C_p ( ∫_0^t_n( S(t_n -s) - S_k,h(t_n -s) )
        (g(s)-g(t_n)) ^2_ℒ_2^0s)^1/2
          =:C_p (J_4^n+J_5^n).

  For J^n_4 we first apply (<ref>).
  Then Assumption <ref> yields
  
    J_4^n ≤ C (h^1+r+k^1+r/2)A^r/2
        g(t_n)_ℒ_2^0≤ C (h^1+r+k^1+r/2) A^r/2 g
        _C([0,T];ℒ_2^0).

  For the estimate of J^n_5 we make use of (<ref>)
  with ρ = 1 + r and of the Hölder continuity of g. This gives 
  
    J_5^n 
           ≤ C (h^1+r+k^1+r/2)
        ( ∫_0^t_n (t_n-s)^-(1+r)g(s)-g(t_n)^2_ℒ_2^0s)^1/2
       ≤ C  g _C^1/2([0,T];ℒ_2^0) (h^1+r+k^1+r/2)
        ( ∫_0^t_n(t_n-s)^-rs)^1/2
       ≤ C T^1-r/2 g _C^1/2([0,T];ℒ_2^0)
        (h^1+r+k^1+r/2).

  Combining the two estimates then yields the assertion.



  As also discussed in <cit.>, 
  the result of Lemma <ref> does not
  hold true in the border case r=1. The reason for this is that the singularity
  caused by the error estimate (<ref>) is no longer integrable
  for r = 1. However, observe that this problem does not occur if 
  g is constant since the term J_5^n is then equal to zero or if g is
  Hölder continuous with an exponent larger than 1/2.


Finally, it remains to estimate the last term on the right hand side of
(<ref>). To this end, we first insert the definition (<ref>)
and obtain for every n ∈{1,…,N_k} 

    ∑_j = 1^n S_k,h^n-j( - ∫_t_j-1^t_j S_k,h f(s,X(s)) s
        + ∫_t_j-1^t_j S_k,h g(s) W(s)
                         - Φ_k,h^j( X(t_j-1), τ_j) )
        _L^p(Ω;H)
         ≤∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( f(s,X(s)) - f(t^τ_j, X(t^τ_j)) )
        s_L^p(Ω;H)
            +  k ∑_j = 1^n S_k,h^n-j+1( f(t^τ_j, X( t^τ_j) ) -
        f ( t^τ_j, Ψ_h,k^j(X(t_j-1),τ_j) ) )
        _L^p(Ω;H)
            + ∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( g(s) -g(t^τ_j) ) W(s)_L^p(Ω;H),

where we recall that t^τ_j=t_j-1+τ_j k.
In the following we derive estimates for each term on the right hand side
separately. The estimate of the first term is related to a randomized
quadrature rule for Hilbert space valued stochastic processes. We refer to
<cit.> for the origin of such quadrature rules.
The presented proof is an adaptation of similar results from <cit.>. Observe that classical methods require additional smoothness of
the mapping f  H → H in order to derive the same convergence rates.
Compare further with <cit.>.


  
  Let Assumption <ref> to Assumption <ref> be fulfilled with 
  p ∈ [2,∞) and γ∈ (0,1/2].
  Then there exists C ∈ (0,∞) such that
  for every h ∈ (0,1), k ∈ (0,T)
  
    max_n∈{1,…,N_k}∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( f(s,X(s)) -f(t^τ_j,
        X(t^τ_j)) ) s_L^p(Ω;H)
         ≤ C(1+  X _C^1/2([0,T];L^p(Ω_W;H))) 
        k^γ+1/2,

  where t^τ_j=t_j-1+τ_j k.



  Due to (<ref>) we have X∈ L^p([0,T]×Ω_W;H).
  From the linear growth (<ref>) of f it then follows 
  that there exists a null set 𝒩_0 ∈^W such that for all
  ω∈𝒩_0^c = Ω_W ∖𝒩_0 we have
  ∫_0^T f(s,X(s,ω))^p s < ∞.  
  Let us therefore fix an arbitrary realization ω∈𝒩_0^c.
  Then for every j ∈{1,…,N_k} we obtain
  
    ∫_t_j-1^t_j f(s,X(s,ω)) s    = k ∫_0^1 f (t_j-1 + s k, X(t_j-1 + s k,ω) )
        s
       = k _τ[  f ( t^τ_j, X(t^τ_j,ω)) ],

  due to t^τ_j ∼𝒰(t_j-1,t_j). 
  
  Next, we define a discrete-time
  error process (E^n)_n ∈{0,1,…,N_k} by setting E^0 ≡ 0 ∈
  H. Further, for every n ∈{1,…,N_k} we set
  
    E^n := ∑_j=1^n S_k,h^n-j+1(∫_t_j-1^t_j f( s,X(s,ω) ) s
        - k f ( t^τ_j,X( t^τ_j,ω) ) ).

  In addition, for every n ∈{1,…,N_k} and m ∈{0,…,n}
  we define M^m_n := S_k,h^n-m E^m, which is evidently
  an H-valued random variable on the product probability space (Ω, ,
  ). In particular, M_n := (M^m_n)_m ∈{0,…,n}⊂
  L^p(Ω;H). From E^n = M^n_n we immediately obtain the estimate
  
    E^n _L^p(Ω;H)   ≤max_m ∈{0,…,n}  M^m_n _L^p(Ω)
       = ( ∫_Ω_Wmax_m ∈{0,…,n}  M^m_n(ω, ·) _L^p(Ω_τ)^p 
          _W(ω) )^1/p

  for all n ∈{1,…,N_k}. Moreover, for each fixed ω∈𝒩_0^c we observe that the mapping 
  M^m_n(ω,·) Ω_τ→ H
  is _m^τ-measurable. Further, for each pair of m_1,m_2 ∈_0 with
  0 ≤ m_1 < m_2 ≤ n it holds true that 
  
    _τ[ M^m_2_n(ω, ·)-M^m_1_n(ω, ·) |
        ^τ_m_1]
          =  ∑_j=m_1+1^m_2 S_k,h^n-j+1_τ[ ∫_t_j-1^t_j f(s, X(s,ω)) s
        - k f( t^τ_j,X(t^τ_j, ω)) | ^τ_m_1]
          = ∑_j=m_1+1^m_2 S_k,h^n-j+1_τ[ ∫_t_j-1^t_j f(s, X(s,ω)) s
        - kf(t^τ_j ,X(t^τ_j,ω))
        ] = 0,

  since τ_j is independent of ^τ_m_1 for every j > m_1.
  Consequently, for every fixed ω∈𝒩_0^c, 
  the process M_n(ω, ·) = (M^m_n(ω, ·))_m ∈{0,…,n} is an (^τ_m)_m ∈{0,…,n}-adapted L^p(Ω_τ;H)-martingale.
  Thus,  the discrete-time version of the Burkholder–Davis–Gundy inequality,
  Proposition <ref>, is applicable and yields 
  
    max_m ∈{0,…,n} M^m_n(ω, ·)
        _L^p(Ω_τ)≤ C_p  [M_n(ω,·)]^1/2_n_L^p(Ω_τ)   for every ω∈𝒩_0^c.

  Next, we insert this and the quadratic variation of M_n(ω, ·) 
  into (<ref>). An application of (<ref>) then yields
  
    E^n _L^p(Ω;H)   ≤ C_p ( ∫_Ω_W_τ[ ( ∑_j = 1^nS_k,h^n-j+1^2_ℒ(H)×
               ∫_t_j-1^t_j f( s,X(s,ω) ) s
        - k f(t^τ_j,X(t^τ_j,ω))
        ^2 )^p/2] 
        _W(ω) )^1/p
       ≤ C_p ∑_j = 1^n∫_t_j-1^t_j( f (s,X(s) )
        - f ( t^τ_j, X(t^τ_j) ) )  s^2  _L^p/2(Ω)^1/2
       ≤ C_p ( ∑_j = 1^n( ∫_t_j-1^t_j f(s,X(s))
        - f (t^τ_j,X(t^τ_j)) _L^p(Ω;H)s)^2 )^1/2,

  where the last step follows from an application of the triangle
  inequality for the L^p/2(Ω)-norm. 
  Next, we make use of Assumption <ref> and obtain for every s ∈
  [t_j-1, t_j], j ∈{1,…,n}, the bound
  
    f (s,X(s)) - f (t^τ_j,X(t^τ_j) ) ≤ C_f ( 1 +  X(s) ) | s - t^τ_j|^γ
        + C_f  X(s) - X(t^τ_j) ,

  which together with (<ref>) implies 
  
    f (s,X(s)) - f (t^τ_j,X(t^τ_j) ) _L^p(Ω;H)
         ≤ C_f ( 1 + sup_t ∈ [0,T] X(t) _L^p(Ω_W;H)) k^γ  + C C_f   X _C^1/2([0,T];L^p(Ω_W;H)) 
        k^1/2.

  Altogether, this shows 
  
    E^n _L^p(Ω;H)≤ C_p C_f (1 + C ) T^1/2(1 +   X _C^1/2([0,T];L^p(Ω_W;H))) 
        k^γ + 1/2.

  This completes the proof.


Let us now turn to the second term on the right hand side of 
(<ref>).


  
  Let Assumption <ref> to Assumption <ref> be fulfilled with 
  p ∈ [2,∞) and r ∈ [0,1). Then there exists C ∈ (0,∞)
  such that for every k ∈ (0,T), h ∈ (0,1) we have
  
    max_n ∈{1,…,N_k} k ∑_j = 1^n S_k,h^n-j+1( f(t^τ_j, X(t^τ_j)) -
          f(t^τ_j,Ψ_h,k^j(X(t_j-1),τ_j)) 
          ) _L^p(Ω;H)
               ≤ C (1 + sup_t ∈ [0,T] X(t)
          _L^p(Ω_W;Ḣ^1+r)) 
          ( h^1+r + k^1+r/2),

  where t^τ_j=t_j-1+τ_j k. 



  First, we fix arbitrary parameter values for h ∈ (0,1) and k ∈ (0,T).
  Then, applications of the triangle inequality, the stability estimate
  (<ref>), and Assumption <ref> yield the estimate
  
    max_n ∈{1,…,N_k} k ∑_j = 1^n S_k,h^n-j+1( f(t^τ_j, X(t^τ_j)) 
        - f(t^τ_j,Ψ_h,k^j(X(t_j-1),τ_j)) 
        ) _L^p(Ω;H)
         ≤ C_f k ∑_j = 1^N_k X(t^τ_j) - Ψ_h,k^j(X(t_j-1),τ_j)
        _L^p(Ω;H).

  Next, let j ∈{1,…,N_k} be arbitrary. After inserting the 
  variation of constants formula (<ref>) for the mild solution
  and the definition (<ref>) of Ψ^j_h,k we obtain
  
    X(t^τ_j) - Ψ_h,k^j(X(t_j-1),τ_j)
          _L^p(Ω;H)
         ≤ ( S(τ_j k) - S_τ_j k, h ) X(t_j-1)
          _L^p(Ω;H)
            + ∫_t_j-1^t^τ_j( S(t^τ_j - s) f(s,X(s))
          - S_τ_j k,h f(t_j-1, X(t_j-1)) ) s_L^p(Ω;H)
            + ∫_t_j-1^t^τ_j( S(t^τ_j-s) g(s)
          -S_τ_j k,h g(t_j-1) ) W(s)_L^p(Ω;H)
          =: J_6^j + J_7^j + J_8^j.
 
  We estimate the three terms separately. 

  The estimate of J_6^j follows at once from Lemma <ref> by
  taking note of
  
    J_6^j    = ( _τ[   ( S(τ_j k) - S_τ_j k, h ) X(t_j-1)
        _L^p(Ω_W;H)^p ] )^1/p
        = ( 1/k∫_0^k ( S(θ) - S_θ, h ) X(t_j-1) _L^p(Ω_W;H)^p
        θ)^1/p
       ≤ C ( h^1+r + k^1+r/2) 
         X _C([0,T];L^p(Ω_W;Ḣ^1+r)).

  For the estimate of J_7^j it is sufficient to note that the integrand is
  bounded uniformly for all k ∈ (0,T) and h ∈ (0,1) due to 
  (<ref>), (<ref>), the linear growth (<ref>)
  of f, and (<ref>). From this we obtain
  
    J_7^j
             ≤( _τ[ (
          ∫_t_j-1^t^τ_j f(s, X(s) ) _L^p(Ω_W;H)
          +  f(t_j-1, X(t_j-1) ) _L^p(Ω_W;H)s)^p ]
          )^1/p
       ≤ 2 Ĉ_f 
          (1 +  X _C([0,T];L^p(Ω_W;H))) k.

  For the estimate of J_8^j we first add and subtract a term. This leads to 
  
    J_8^j
           ≤∫_t_j-1^t^τ_j 
        S(t^τ_j-s) ( g(s)-g(t_j-1) ) W(s)_L^p(Ω;H)
          + ∫_t_j-1^t^τ_j( S(t^τ_j-s) - S_τ_j k,h) g(t_j-1) W(s)_L^p(Ω;H),

  Then we apply
  Proposition <ref>, Assumption <ref>, and estimates
  (<ref>) and (<ref>). Altogether, this yields 
  
    J_8^j
           ≤( _τ[ ∫_t_j-1^t^τ_j 
        S(t^τ_j-s) ( g(s)-g(t_j-1) )
        W(s)^p_L^p(Ω_W;H)] )^1/p
          + ( _τ[ ∫_t_j-1^t^τ_j( S(t^τ_j-s) - S_τ_j k,h)
        g(t_j-1) W(s)^p_L^p(Ω_W;H)] 
        )^1/p
       ≤ C_p ( _τ[ (∫_t_j-1^t^τ_j g(s)-g(t_j-1)_ℒ_2^0^2 s)^p/2])^1/p
          + C_p ( _τ[ ( ∫_t_j-1^t^τ_j (S(t^τ_j-s)- S_τ_j k,h ) 
        g(t_j-1) _ℒ_2^0^2
        s)^p/2] )^1/p
       ≤ C_p  g _C^1/2([0,T];ℒ_2^0)
        k + C_p  A^r/2 g
        _C([0,T];ℒ_2^0)( h^1+r +  k^1+r/2).

  Inserting the estimates for J_6^j, J_7^j, and J_8^j into
  (<ref>) then completes the proof.


In contrast to the drift-randomized Milstein method studied in
<cit.>, we also randomize the diffusion term
in the method (<ref>). As our final lemma shows, this allows us to
only impose a smoothness condition on g with respect to the norm
W^1/2 + γ,p(0,T;ℒ_2^0) instead of the 
more restrictive Hölder norm C^1/2 + γ([0,T];
ℒ_2^0), γ∈ [0,1/2] usually found in the
literature. We refer to <cit.> for further quadrature rules 
which apply to stochastic integrals, whose regularity is measured in terms
of fractional Sobolev spaces.


  
  Let Assumption <ref> be fulfilled. 
  For every g ∈ W^ν,p(0,T;ℒ_2^0) with ν∈ (0,1], p
  ∈ [2,∞) and for every k ∈ (0,T), h ∈ (0,1) it holds true that
  
    max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( g(s) -  g( t_j-1 + τ_j k) ) W(s)_L^p(Ω;H)
           ≤ C_p T^p-2/2p g
          _W^ν,p(0,T;ℒ_2^0) k^ν.




  Fix arbitrary parameter values k ∈ (0,T), h ∈ (0,1). As in the proof
  of Lemma <ref> we introduce the process ĝ [0,T] ×Ω_τ→ℒ_2^0 defined by
  
    ĝ(t):=g(t_j-1 + τ_j k),    for  t ∈ [t_j-1,t_j),
          j ∈{1,…,N_k}.

  After inserting this and the piecewise constant interpolation
  S_k,h of S_k,h into the left hand side of
  (<ref>), we obtain for 
  every n ∈{1,…,N_k}
  
    ∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( g(s) - g(t_j-1 + τ_j k) ) W(s)_L^p(Ω;H)
           = ∫_0^t_nS_k,h(t_n-s) 
          ( g(s) - ĝ(s) ) W(s)_L^p(Ω;H)
         ≤ C_p ( ∑_j = 1^n ∫_t_j-1^t_j g(s) - g(t_j-1 + τ_j k)
          ^2_L^p(Ω_τ;ℒ_2^0)s)^1/2,

  where we applied Proposition <ref> and the stability estimate
  (<ref>) in the last step. After two applications of the
  Hölder inequality with exponents 1/ρ + 1/ρ' = 1
  and ρ = p/2 we arrive at
  
    ∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( g(s) - g(t_j-1 + τ_j k) ) 
          W(s)_L^p(Ω;H)
         ≤ C_p 
          ( ∑_j = 1^n k^1 - 2/p( ∫_t_j-1^t_j g(s) - g(t_j-1 + τ_j k) 
          ^p_L^p(Ω_τ;ℒ_2^0)s)^2/p)^1/2
         ≤ C_p T^p - 2/2p( ∑_j = 1^n
          ∫_t_j-1^t_j g(s) - g(t_j-1 + τ_j k) 
          ^p_L^p(Ω_τ;ℒ_2^0)s)^1/p
          = C_p T^p - 2/2p( ∑_j = 1^n
          1/k∫_t_j-1^t_j∫_t_j-1^t_j g(s) - g(θ) _ℒ_2^0^p θs)^1/p.

  First, we discuss the
  case ν = 1, that is g ∈ W^1,p(0,T;ℒ_2^0).
  Under this condition, there exists an absolutely continuous representative in
  the equivalence class of g, for which we obtain the estimate
  
    1/k∫_t_j-1^t_j∫_t_j-1^t_j g(s) - g(θ) _ℒ_2^0^p θs   = 1/k∫_t_j-1^t_j∫_t_j-1^t_j∫_θ^s g'(z) z_ℒ_2^0^p 
        θs
       ≤ k^p-2∫_t_j-1^t_j∫_t_j-1^t_j∫_t_j-1^t_j g'(z) _ℒ_2^0^p zθs
       = k^p ∫_t_j-1^t_j g'(z) _ℒ_2^0^p z

  for all j ∈{1,…,n}.
  Inserting this into (<ref>) then yields the desired estimate for
  ν = 1.

  For the case ν∈ (0,1) we recall the definition
  (<ref>) of the fractional Sobolev–Slobodeckij norm. Then the
  estimate of (<ref>) is continued as follows
  
    1/k∑_j = 1^n ∫_t_j-1^t_j∫_t_j-1^t_j g(s) - g(θ) _ℒ_2^0^p θs   ≤ k^ν p∑_j = 1^n ∫_t_j-1^t_j∫_t_j-1^t_j g(s) - g(θ) _ℒ_2^0^p/|s -
        θ|^1 + ν pθs
       ≤ k^ν p g ^p_W^ν,p(0,T;ℒ_2^0).

  Together with (<ref>) this completes the proof.


Combining Lemma <ref> to Lemma <ref> gives immediately
an estimate for (<ref>), which is essentially the
consistency of the numerical scheme (<ref>): 


  
  Let Assumptions <ref> to <ref> be fulfilled for some p ∈
  [2,∞), r ∈ [0,1), and γ∈ (0,1/2]. 
  Then there exists a constant C ∈ (0,∞) such that for every h ∈
  (0,1) and k ∈ (0,T) 
  
    ℛ_k,h[ X|_π_k ] _S,p,h
        ≤ C (1 +  X _C([0,T];L^p(Ω_W;Ḣ^1+r))
        +  X _C^1/2([0,T];L^p(Ω_W;H)))
        ( h^1+r + k^1/2 + min(r/2, γ)),

  where X denotes the mild solution (<ref>) to the stochastic
  evolution equation (<ref>).


Now we are ready to address the proof of the convergence result, 
Theorem <ref>. 


  From Theorem <ref>, we obtain the bistability of the numerical scheme
  (<ref>). By simply choosing Y_k=X|_π_k and
  Z_k=X_k,h in (<ref>) we get
  
    X|_π_k -X_k,h_∞,p≤ C_Stabℛ_k,h[X|_π_k] - ℛ_k,h[X_k,h] 
        _S,p,h.

  Next, we take note of ℛ_k,h[X_k,h] =0.
  The final assertion then follows from an application 
  of Theorem <ref>. 


 
  Let us briefly discuss the case of multiplicative noise.
  To be more precise, instead of (<ref>), we want to approximate the
  mild solution to a semilinear stochastic evolution equation of the form 
  
    X(t) + [ A X(t) + f(t,X(t)) ] t
          = g(X(t)) W(t),     for  t ∈ (0,T],
    
          X(0) = X_0,

  with all assumptions in Section <ref> remaining the same expect
  for Assumption <ref> on g, which is
  replaced, for instance, by <cit.>. 
  Due to the multiplicative noise, we cannot randomize the stochastic integral
  in the same way as in the additive noise case.
  However, one can still benefit from a randomization of the semilinear drift
  part. More precisely, for each k ∈ (0,T)
  and h ∈ (0,1), the drift-randomized Milstein–Galerkin finite
  element method is given by 
  
    X_k,h^n, τ   = S_τ_n k, h[ X_k,h^n-1 
          - τ_n k f(t_n-1, X_k,h^n-1) + g(X_k,h^n-1)
          Δ_τ_n k W(t_n-1) ],
    
          X_k,h^n   = S_k,h[ X_k,h^n-1
          - k f(t_n^τ, X_k,h^n, τ) + g(X_k,h^n-1)
          Δ_k W(t_n-1) 
          +∫_t_n-1^t_n g'(X_k,h^n-1) 
          [ ∫_t_n-1^σ_1g(X_k,h^n-1) W(σ_2)] W(σ_1)]

  for all n ∈{1,…,N_k} with initial value X_k,h^0 = P_h X_0.
  All lemmas on the consistency remain valid with the exception of 
  Lemma <ref>, Lemma <ref>, and Lemma <ref>.
  Instead, we can borrow <cit.>. An additional modification
  of the estimate of J_8^j in Lemma <ref>
  then yields the same rate of convergence as in Theorem <ref>.
  Hence, the convergence result in Theorem <ref> carries over to the
  multiplicative noise case for the scheme (<ref>). As in
  finite dimensions in <cit.>, the drift-randomization
  technique therefore reduces the regularity conditions on the drift
  semilinearity f 
  significantly. In particular, we do not impose a differentiability condition
  on the mapping f required in <cit.>.
  However, note that the conditions imposed on g in <cit.> are rather restrictive. 








§ INCORPORATING A NOISE APPROXIMATION



In this section we discuss the numerical approximation of a 
Q-Wiener process W [0,T]×Ω_W → U with values in
a separable Hilbert space (U, (·, ·)_U, ·_U). 
In particular, we investigate how the stability and consistency of the
numerical scheme (<ref>) will be
affected if the noise approximation is incorporated.
Hereby, we follow similar arguments as in <cit.>, which
in turn is based on <cit.>. 

Our analysis relies on a spectral approximation of the Wiener process. Because
of this we impose the following stronger assumption on the covariance operator 
Q.


  
  The covariance operator Q∈ℒ(U) is symmetric,
  nonnegative, and of finite trace.


It directly follows from Assumption <ref> that Q is a compact operator.
Moreover, the spectral theorem for compact operators then ensures
the existence of an orthonormal basis (φ_j)_j∈ℕ of the
separable Hilbert space U such that 

    Qφ_j =μ_j φ_j, j∈ℕ,

where (μ_j)_j ∈ are the eigenvalues of Q. 

For each M ∈ we introduce a truncated version Q_M ∈ℒ(U) of the covariance operator Q determined by 

    Q_M φ_j 
        := 
        μ_j φ_j,       if  j ∈{1,…,M},
    
          0,       else.
 
In the following, we further use the abbreviation Q_cM:=Q-Q_M. Note that
Q_M is of finite rank. We then define a Q_M-Wiener process
W^M[0,T]×Ω_W → U by 

    W^M(t) := ∑_j=1^M √(μ_j)β_j(t)φ_j, t∈
      [0,T],

where β_j [0,T]×Ω_W →ℝ,
j∈ℕ, is an independent family of standard real-valued Brownian
motions. The link between W^M and the original Q-Wiener process W is 
established by the relationship β_j(t) = 1/√(μ_j)
( W(t), φ_j)_U, as in <cit.>.
Moreover, we also define W^cM:=W-W^M. Note that W^cM 
is a Q_cM-Wiener process and possesses the spectral representation

    W^cM(t) := ∑_j=M+1^∞√(μ_j)β_j(t) φ_j,
       t∈ [0,T].


We now introduce a modification of the numerical scheme (<ref>),
which only uses the increments of the truncated Wiener process W^M. 
For every k ∈ (0,T), h ∈ (0,1), and M ∈
let the initial value be given by X_k,h,M^0 := P_h X_0. 
In addition, the random variables
X_k,h,M^n, n ∈{1,…,N_k}, are defined by the recursion   

    X_k,h,M^n, τ   = S_τ_n k,h[ X_k,h,M^n-1
        - τ_n k f(t_n-1, X_k,h,M^n-1)
        + g(t_n-1) Δ_τ_n k W^M(t_n-1) ],
    
        X_k,h,M^n   = S_k,h[ X_k,h,M^n-1 
        - k f(t_n-1 + τ_n k, X_k,h,M^n, τ) 
        + g(t_n-1 + τ_n k) Δ_k W^M(t_n-1) ]

for n ∈{1,…,N_k}, where the linear operators S_k,h∈ℒ(H) are defined
in (<ref>) and the random variables (τ_n)_n ∈ are
𝒰(0,1)-distributed and independent from each other as well as from
the Wiener process W. As in (<ref>) the Wiener increments are given
by  

    Δ_κ W^M(t) := W^M(t + κ) - W^M(t)

for all t ∈ [0,T) and κ∈ (0,T-t).

Analogously to (<ref>) and
(<ref>) the increment functions
Φ_k,h,M^j  H × [0,1] ×Ω_W → H and
Ψ_k,h,M^j  H × [0,1] ×Ω_W → H
are then given by

    Φ_k,h,M^j(x,τ) 
           := -k S_k,h f(t_j-1 + τ k,
        Ψ^j_k,h,M(x,τ)) 
        
        + S_k,h g(t_j-1 + τ k) Δ_k W^M(t_j-1)

and

    Ψ_k,h,M^j(x,τ) := S_τ k,h[ x - τ k f(t_j-1,
      x ) + g(t_j-1) Δ_τ k W^M(t_j-1) ]

for all x ∈ H and τ∈ [0,1].

We first study the stability of the truncated scheme
(<ref>). 


  
  Let Assumptions <ref> to <ref> and Assumption <ref> be
  satisfied. Then, for every k ∈ (0,T), h ∈ (0,1), and M ∈
  the numerical scheme (<ref>) is bistable with respect to
  the norms ·_∞,p and ·_S,p,h.  
  In particular, the stability constant C_Stab can be
  chosen independently of M ∈. 




  For the proof we observe that all estimates in Lemma <ref> hold also
  true for the noise truncated scheme (<ref>). To prove 
  the independence of the stability constant of the parameter M ∈ 
  we recall that the family of eigenfunctions (φ_j)_j ∈ of Q
  is an orthonormal basis of U. The Hilbert–Schmidt norm of g(t) ∘
  Q^1/2_M is therefore bounded by 
  
    g(t)Q_M^1/2^2_ℒ_2
          = ∑_j=1^Mμ_j  g(t) φ_j^2
          ≤∑_j=1^∞μ_j g(t) φ_j^2
          =  g(t) _ℒ_2^0^2,

  for all t ∈ [0,T] and M ∈. Inserting this into
  (<ref>) yields the assertion.



It remains to address the consistency of the noise truncated scheme
(<ref>).
For this we first introduce the associated residual operator
ℛ_k,h,M𝒢_k^p →𝒢_k^p given by 

    ℛ_k,h,M[Z_k](t_0) := Z_k^0 - ξ_h,   
    ℛ_k,h,M[Z_k](t_n) := Z_k^n - S_k,h Z_k^n-1 -
        Φ_k,h,M^n(Z_k^n-1,τ_n),    n ∈{1,…,N_k},

for each grid function Z_k ∈𝒢_k^p.
In order to control the truncation error with respect to the parameter
M∈ℕ we need the following additional assumption:


  
  Let (φ_m)_m ∈⊂ U and (μ_m)_m ∈⊂
  [0,∞) be the families of 
  eigenfunctions and eigenvalues of Q from 
  (<ref>). We assume the existence of constants
  C_Q, α∈ (0,∞) such that
  
    sup_t ∈ [0,T](∑_m=1^∞ m^2 αμ_m
        g(t)φ_m^2 )^1/2≤ C_Q.
 
 

We now state the consistency result for the noise truncated scheme
(<ref>).


  
  Let Assumptions <ref> to <ref> be fulfilled for some p ∈
  [2,∞), r ∈ [0,1), and γ∈ (0,1/2]. Let Assumptions
  <ref> and <ref> be fulfilled with α∈ (0,∞).
  Then there exists a constant C ∈ (0,∞) such that
  for every k ∈ (0,T), h ∈ (0,1), and M ∈ we have 
  
    ℛ_k,h,M[ X|_π_k ] _S,p,h≤ C 
        ( h^1+r + k^1/2 + min(r/2, γ)+M^-α),

  where X|_π_k denotes the restriction of the 
  mild solution (<ref>) to the equidistant grid points in
  π_k. 



  An inspection shows that Lemma <ref> to
  Lemma <ref> remain valid for the truncated scheme
  (<ref>) by applying, if necessary, 
  the same argument as in (<ref>). It therefore remains to
  adapt the proofs of Lemma <ref> and Lemma <ref>.

  First, we observe that the truncation of the noise only affects
  the term J_8^j appearing in (<ref>) 
  in the proof of Lemma <ref>. Hence, we need to find 
  a corresponding estimate of the term 
  
    J_8,M^j    = ∫_t_j-1^t_j^τ S(t_j^τ -s) g(s)
        W(s) - S_τ_j k, h g(t_j-1) Δ_τ_jk W^M(t_j-1)
        _L^p(Ω;H)
       ≤∫_t_j-1^t_j^τ 
        S(t_j^τ -s) ( g(s) - g(t_j-1) )
        W(s)_L^p(Ω;H)
            + ∫_t_j-1^t_j^τ( S(t_j^τ -s) - S_τ_j k, h(t_j^τ  -s ) )
        g(t_j-1) W(s)_L^p(Ω;H)
             +  S_τ_j k, h g(t_j-1) Δ_τ_j k
        W^cM(t_j-1) _L^p(Ω;H)

  for all j ∈{1,…,N_k}, where t_j^τ = t_j-1 + τ_j k. The
  first two terms are estimated in the same way as J_8^j in the proof of
  Lemma <ref>. 
  In order to give a bound for the last term we first take note of
  
    g(t) Q_cM^1/2^2_ℒ_2   = ∑_m=M+1^∞μ_m g(t)φ_m^2
          ≤∑_m=M+1^∞m^2α/M^2αμ_mg(t)φ_m^2
       ≤1/M^2α∑_m=1^∞ 
          m^2αμ_m g(t)φ_m^2
          ≤C^2_Q/M^2α

  for all t ∈ [0,T],
  which follows from Assumption <ref>. Together with applications
  of Proposition <ref> and the stability estimate (<ref>)
  from Lemma <ref> we then obtain
  
    S_τ_j k, h g(t_j-1) Δ_τ_j k
        W^cM(t_j-1) _L^p(Ω;H)   ≤ C_p ( 
        _τ[ ( ∫_t_j-1^t_j^τ g(t_j-1) Q_cM^1/2^2_ℒ_2s)^p/2] )^1/p
       ≤ C_p C_Q M^- α k^1/2.

  Altogether, this yields
  
    J_8,M^j ≤  C_p ( 
         g _C^1/2([0,T];ℒ_2^0) k 
        +  A^r/2 g
        _C([0,T];ℒ_2^0)( h^1+r +  k^1+r/2)
        + C_Q M^- α k^1/2).

  In the same way we derive a modification of Lemma <ref>. To be more
  precise, instead of (<ref>) we need to find a bound for
  the norm
  
    max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j+1[
        ∫_t_j-1^t_j g(s) W(s)
        - g(t_j-1 + τ_j k) Δ_k W^M(t_j-1) ]
        _L^p(Ω;H)
         ≤max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j+1∫_t_j-1^t_j( g(s) - g(t_j-1 + τ_j k) ) W(s)_L^p(Ω;H)
            + max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j+1 g(t_j-1 + τ_j k)
        Δ_k W^cM(t_j-1) _L^p(Ω;H).

  Lemma <ref> is applicable to the first term on the right hand
  side and it remains to give a bound for the second term. To this end, we
  recall the operators S_k,h and ĝ defined in
  (<ref>) and (<ref>), respectively. Inserting these
  operators then yields
  
    max_n ∈{1,…,N_k}∑_j = 1^n S_k,h^n-j+1 g(t_j-1 + τ_j k)
        Δ_k W^cM(t_j-1) _L^p(Ω;H)
          = max_n ∈{1,…,N_k}∫_0^t_nS_k,h(t_n - s)
        ĝ(s) W^cM(s)_L^p(Ω;H)
         ≤ C_p ( ∫_0^T ĝ(s) Q^1/2_cM^2_L^p(Ω_τ;ℒ_2)s)^1/2,

  where we also applied Proposition <ref> and the stability estimate
  (<ref>)  in the last step. After reinserting the definition
  (<ref>) of ĝ we again make use of (<ref>) and
  obtain
  
    C_p ( ∫_0^T ĝ(s) Q^1/2_cM^2_L^p(Ω_τ;ℒ_2)s)^1/2   = C_p ( k ∑_j = 1^N_k g(t_j-1 + τ_j k) Q_cM^1/2_L^p(Ω_τ;ℒ_2)^2 
        )^1/2
       ≤ C_p C_Q T^1/2 M^-α.

  Together with a simple modification of (<ref>), a combination 
  of these estimates with Lemma <ref> 
  to Lemma <ref> then yields the assertion.



  
  From the proof of Theorem <ref> it follows that
  the value of the parameter M ∈ does not need to be the
  same in the definitions (<ref>) and 
  (<ref>) of the two increment functions
  Φ_k,h,M and Ψ_k,h,M. In fact, we obtain the same order of
  convergence if we replace M in the definition of Ψ_k,h,M by
  M̃ := ⌊√(M)⌋ + 1.
  
  The reason for this is that
  the parameter M̃ only appears in the estimate (<ref>)
  of J_8,M^j with k^1/2 as a pre-factor. 
  Due to
  
    M̃^-α k^1/2≤1/2( M̃^-2 α + k )
        ≤1/2( M^-α + k )

  the orders of convergence with respect to M and k remains indeed
  unaffected by this modification. 
  This leads to a substantial reduction of the computational cost for
  Ψ_k,h,M for large values of M. 



  Assumption <ref> can be relaxed to covariance operators
  Q ∈ℒ(U) which are not of finite trace but still
  possess a spectral representation of the form (<ref>). For example,
  we mention the case of a space-time white noise W, where Q =
  Id. In that case the approximation W^M in 
  (<ref>) does not, in general, converge to W with respect
  to the norm in the Hilbert space U. 
  Nevertheless, if the mapping g still satisfies Assumption <ref>
  then Theorem <ref> and Theorem <ref>
  remain valid. In particular, observe that (<ref>) simply 
  yields the Hilbert–Schmidt norm of g(t) ∈ℒ_2(U,H)
  if μ_m ≡ 1 and α = 0.

  For further details on the analytical treatment of cylindrical Wiener
  processes we refer to <cit.>.










§ APPLICATION TO STOCHASTIC PARTIAL DIFFERENTIAL EQUATIONS



In this section we apply the randomized method (<ref>) for
the numerical solution of a semilinear stochastic partial differential
equation (SPDE). First, we reformulate the SPDE as a stochastic evolution
equation of the form (<ref>). Then we verify the conditions of
Theorem <ref>. Finally, we perform a numerical 
experiment. 

Let us first introduce the semilinear SPDE that we want to solve numerically in
this section. The goal is to find a measurable mapping u  [0,T] ×
[0,1] ×Ω_W →ℝ satisfying    

    u(t,x) = [∂^2/∂ x^2 u(t,x)
          + η( t, u(t,x) )
          ] t + σ( t ) W(t,x),
      u(0,x) = u_0(x):= 2(1-x)x,   u(t,0)=u(t,1)=0,

for t ∈ (0,T] and x∈ (0,1). The Wiener process W is assumed
to be of trace class and will be specified in more detail further below.
The coefficient function σ [0,T] → [0,∞) is used to
control the noise intensity, where we require that σ∈
C^1/2(0,T) ∩ W^1/2+ γ,p(0,T) for some γ∈
(0,1/2] and p ∈ [2,∞).
The drift function η [0,T] ×→ is
assumed to be continuous. In addition, there exists L ∈ (0,∞) and
γ∈ (0,1] such that

    | η(t,v_1) - η(t,v_2) |    ≤ L |v_1 - v_2|
    
        | η(t_1,v) - η(t_2,v) |    ≤ L (1 + |v|) |t_1 - t_2|^γ

for all t, t_1, t_2 ∈ [0,T], v, v_1, v_2 ∈.

In order to rewrite the SPDE (<ref>) as a stochastic evolution
equation we consider the separable Hilbert space H=L^2(0,1). 
Then, the operator -A = ∂^2/∂ x^2 is the Laplace
operator on (0,1) with homogeneous Dirichlet conditions. It is well-known,
see <cit.> or 
<cit.>, that
this operator
satisfies Assumption <ref>. We have that (A) = H^1_0(0,1) ∩ 
H^2(0,1). Moreover, the operator A has the eigenfunctions e_j =
√(2)sin(jπ·) and eigenvalues λ_j = j^2π^2 for j∈ℕ.   

The initial condition U_0 ∈ H is then given by

    U_0(x) := u_0(x) = x ( 1-x),    x ∈ (0,1).

Evidently, U_0 satisfies Assumption <ref> for any value of p ∈
[2,∞) and r ∈ [0,1]. In particular, we have U_0 ∈(A).

Further, let f  [0,T] × H → H be the Nemytskii operator
induced by η. To be more precise, f is defined by  

    f(t,v)(x) = -η(t,v(x)),   for all  v ∈ H,   x ∈ (0,1).

Then,with the same constant L ∈ (0,∞) as in (<ref>) we have

    f(t,v_1)- f(t,v_2)^2    = 
      ∫_0^1 |η(t,v_1(x)) - η(t,v_2(x)) |^2 x
       ≤  L^2 ∫_0^1 | v_1(x)-v_2(x)|^2 x
      = L^2 v_1-v_2^2

for all v_1, v_2 ∈ H and t ∈ [0,T]. Analogously, we get

    f(t_1,v)- f(t_2,v)≤ L (1 + v ) |t_1 - t_2 |^γ

for all v ∈ H and t_1, t_2 ∈ [0,T].
Thus, Assumption <ref> is satisfied with the same values for 
γ and C_f = L as in (<ref>).

Next, we specify the Wiener process W appearing in (<ref>).
For this we choose U = H = L^2(0,1). Then, the covariance operator Q ∈ℒ(H) is defined by setting Q e_j = μ_j e_j, where μ_j = j^-3 and  
(e_j)_j ∈⊂ H is the orthonormal basis consisting of
eigenfunctions of A. Clearly, we have Tr(Q) = ∑_j = 1^∞μ_j < ∞. Finally, the operator g [0,T]→ℒ_2^0 is 
defined by 

    g(t) := σ(t) Id

for all t ∈ [0,T], where Id∈ℒ(H) is the
identity operator. In order to verify
Assumption <ref> recall that g(t) _ℒ_2^0 =  g(t)
Q^1/2_ℒ_2. Then, for every r ∈ [0,1) we compute

    sup_t∈ [0,T]A^r/2 g(t)^2_ℒ_2^0
      = sup_t∈ [0,T]∑_j=1^∞σ(t) A^r/2√(μ_j) 
      e_j ^2 ≤σ^2_C([0,T])π^2r∑_j=1^∞
      j^2r-3<∞.

In addition, we get for all t_1, t_2 ∈ [0,T] that 

    g(t_1)-g(t_2)_ℒ_2^0^2
      = ∑_j=1^∞(σ(t_1)-σ(t_2)) √(μ_j) e_j ^2
      ≤σ^2_C^1/2([0,T])(Q)|t_1-t_2|.
 
Moreover, since σ∈ W^1/2+γ,p(0,T) 
one can easily validate that 

    ∫_0^T ∫_0^T g(t_1) - g(t_2) ^p_ℒ_2^0/|t_1 -
      t_2|^1 + p(1/2+γ) t_2t_1
       =  (Q) ∫_0^T ∫_0^T 
       | σ(t_1) - σ(t_2) |^p /|t_1 -
      t_2|^1 + p(1/2 + γ)t_2t_1
      < ∞.

This implies that g ∈ W^1/2+γ,p(0,T;ℒ_2^0).
Altogether, we have verified Assumption <ref> with
γ∈ (0,1/2) and r=1-ϵ for any ϵ∈
(0,1/2). By the same means one also verifies Assumption <ref> for
any α∈ (0,1).

Altogether, we can rewrite the SPDE (<ref>) as the following
stochastic evolution equation on H = L^2(0,1) 

    U(t) + [ A U(t) + f(t,U(t)) ] t
          = g(t) W(t),    t ∈ (0,T],
      U(0) = U_0.


Next, we turn to the numerical discretization of (<ref>). For
the spatial discretization we choose a standard finite element method
consisting of piecewise linear functions on a uniform mesh in (0,1).
It is well-known that the associated Ritz projector then satisfies
Assumption <ref>. For instance, we refer to
<cit.>.

Let us now choose the mappings η and σ in (<ref>) 
more explicitly. In the simulations below we used the function η_J → given by

    η_J(v) := ∑_n=0^J a^n cos(b^nπ v),    v ∈, J ∈,

as the semilinearity.
Note that η_J is a truncated version of the Weierstrass function
with parameters a ∈ (0,1) and b ∈ being an odd integer such that
ab>1+3/2π. For such a choice of the parameter values
the Weierstrass function (obtained for J = ∞) 
is everywhere continuous but nowhere differentiable, see <cit.>.
The mapping η_J is therefore a smooth approximation of an irregular
mapping. In particular, η_J is Lipschitz continuous and bounded for every
J ∈, but the Lipschitz constant grows exponentially with J.  

In addition, we performed the numerical experiments
with two different mappings in place of
the noise intensity σ. First, we used the function
σ_1  [0,T] → defined by

    σ_1(t) = 3 √(t),    t ∈ [0,T].

It is easily verified that indeed σ_1 ∈ C^1/2(0,T) ∩ 
W^1 - ϵ,2(0,T) for any ϵ∈ (0,1/2). 
Second, we also made use of the mapping σ_2  [0,T] → given by

    σ_2(t) = 4 √( | sin( 16 π t) | ),    t ∈ [0,T].

Note that this mapping resembles a so called fooling function, that is
particularly tailored to misguide the classical Euler–Galerkin finite element
method. See <cit.> and the references therein 
for a more detailed discussion of fooling functions. 

With these choices of η and σ, all conditions of the randomized
Galerkin finite element method (<ref>) with truncated noise
are satisfied. In particular, we expect a temporal order as high as
1/2+min(r/2,γ) ≈ 1 -ϵ by 
Theorem <ref>. 

For the simulation displayed in Figure <ref>
we chose the parameter values a=0.9, b=7, and J=5 for η_J
in (<ref>). As the final time
we set T = 1. For the spatial
approximation we fixed the equidistant step size h = 1/1001, that is,
we had N_h = 1000 degrees of freedom in the finite element space. For the 
simulation of the Wiener increments we followed the approach
in <cit.>. More precisely, the truncated
Karhunen-Loéve expansion

    W^M(t_n+1,x_i) - W^M(t_n,x_i)
      = ∑_j=1^M√(2 μ_j)sin( π j x_i ) (
      β_j(t_n+1) - β_j(t_n) )
  
can be evaluated efficiently by the discrete sine transformation 
on the nodes x_i := ih, i =1, …,N_h, of the equidistant spatial mesh. 
For simplicity we chose the value M = N_h = 1000 for the expansion
(<ref>). 

Since the purpose of the randomization technique is to improve the temporal
convergence rate compared to Euler–Maruyama-type methods, we focused
on measuring this rate in our numerical experiments. 
For its approximation we first generated a reference solution with a small step
size of k_ref = 2^-12. This reference solution was then compared 
to numerical solutions with larger step sizes k ∈{2^-i  :   i =
4,…,9}. Instead of evaluating directly the norm ·_∞,2
defined in (<ref>) we replaced the integral with respect to Ω
by a Monte Carlo approximation with 100 independent samples and the norm in
H = L^2(0,1) is approximated by a trapezoidal rule. This procedure was used
for the randomized Galerkin finite element method (<ref>) 
as well as for the classical linearly-implicit Euler–Galerkin finite element
method without any artificial randomization. 

   

The results of our simulations are shown in Figure <ref>. We plot
the Monte Carlo estimates of the root-mean-squared errors 
versus the underlying temporal step size, i.e., the
number i on the x-axis indicates the corresponding simulation is based on
the temporal step size k = 2^-i. The figure on the left hand side shows the
results for σ = σ_1 = 3 √(t), while the right hand shows
σ_2(t).
In both subfigures, the sets of data points on the black dotted curves with
triangle markers show the errors of the classical Galerkin finite
element method, while
the red-dotted error curves with diamond markers correspond to 
the simulations of the randomized Galerkin finite element method. 
In addition, we draw a
dashed blue reference line representing a method of order one.
The order numbers displayed in the figure correspond to the slope of a best
fitting function obtained by linear regression. This might be interpreted as
an average order of convergence.

For the case of σ_1, the error curves from both methods are almost
overlapping, with the same average order of convergence 0.84. Since the
randomized method is computationally up to twice as expensive as the classical
Galerkin finite element method, the latter method is clearly superior in this
example. However, as already discussed in Remark <ref>, the
results on the error analysis of the classical method currently available in
the literature are not able to  theoretically explain
the same order of convergence for
general stochastic evolution equations satisfying Assumptions <ref> to
<ref>. 

This is illustrated by the more academic example of σ_2. Here the
coefficient function of the noise intensity is chosen in such a way that the
classical Galerkin finite element method cannot distinguish between a
deterministic PDE without the Wiener noise and the SPDE
(<ref>). In fact, for all step sizes k = 2^i with i ∈{4, 5 } the classical method only evaluates σ_2 on its zeros which
explains the large errors for those step sizes seen on the right hand side
of Figure <ref>. The randomized method is less severely affected by
the highly oscillating coefficient function. For smaller step sizes this
advantage then decays.



Table <ref> contains the numerical values of the error data
displayed in Figure <ref>. In addition, we computed the corresponding
experimental orders of convergence defined by

    = log((2^-i)) - log((2^-i+1)) log(2^-i)-log(2^-i+1)

for i∈{5,6,7,8,9}, where the term (2^-i) denotes the
error of step size 2^-i. 



§ ACKNOWLEDGEMENT


The authors like to thank Sebastian Zachrau for assisting with the numerical 
experiments.
This research was carried out in the framework of Matheon
supported by Einstein Foundation Berlin. The authors also gratefully
acknowledge financial support by the German Research Foundation through the
research unit FOR 2402 – Rough paths, stochastic partial differential
equations and related topics – at TU Berlin. 




10

barth2012
A. Barth and A. Lang.
Milstein approximation for advection-diffusion equations driven by
  multiplicative noncontinuous martingale noises.
Appl. Math. Optim., 66(3):387–413, 2012.

barth2013
A. Barth and A. Lang.
L^p and almost sure convergence of a Milstein scheme for
  stochastic partial differential equations.
Stochastic Process. Appl., 123(5):1563–1587, 2013.

barth2013b
A. Barth, A. Lang, and C. Schwab.
Multilevel Monte Carlo method for parabolic stochastic partial
  differential equations.
BIT, 53(1):3–27, 2013.

becker2016
S. Becker, A. Jentzen, and P. E. Kloeden.
An exponential Wagner–Platen type scheme for SPDEs.
SIAM J. Numer. Anal., 54(4):2389–2426, 2016.

beyn2010
W.-J. Beyn and R. Kruse.
Two-sided error estimates for the stochastic theta method.
Discrete Contin. Dyn. Syst. Ser. B, 14(2):389–407, 2010.

burkholder1991
D. L. Burkholder.
Explorations in martingale theory and its applications.
In École d'Été de Probabilités de Saint-Flour
  XIX—1989, volume 1464 of Lecture Notes in Math., pages 1–66.
  Springer, Berlin, 1991.

daprato1992
G. Da Prato and J. Zabczyk.
Stochastic Equations in Infinite Dimensions, volume 44 of Encyclopedia of Mathematics and its Applications.
Cambridge University Press, Cambridge, 1992.

dinezza2012
E. Di Nezza, G. Palatucci, and E. Valdinoci.
Hitchhiker's guide to the fractional Sobolev spaces.
Bull. Sci. Math., 136(5):521–573, 2012.

eisenmann2017
M. Eisenmann, M. Kovács, R. Kruse, and S. Larsson.
On a randomized backward Euler method for nonlinear evolution
  equations with time-irregular coefficients.
ArXiv preprint, arXiv:1709:01018, 2017.

eisenmann2018
M. Eisenmann and R. Kruse.
Two quadrature rules for stochastic Itō-integrals with
  fractional Sobolev regularity.
ArXiv preprint, arXiv:1712.08152, 2017.

gilbarg2001
D. Gilbarg and N. S. Trudinger.
Elliptic partial differential equations of second order.
Classics in Mathematics. Springer-Verlag, Berlin, 2001.
Reprint of the 1998 edition.

haber1966
S. Haber.
A modified Monte-Carlo quadrature.
Math. Comp., 20:361–368, 1966.

haber1967
S. Haber.
A modified Monte-Carlo quadrature. II.
Math. Comp., 21:388–397, 1967.

hardy1916
G. Hardy.
Weierstrass’s non-differentiable function.
Trans. Amer. Math. Soc, 17(3):301–325, 1916.

hofmanova2017
M. Hofmanová, M. Knöller, and K. Schratz.
Stratified exponential integrator for modulated nonlinear
  Schrödinger equations.
ArXiv preprint, arXiv:1711.01091, 2017.

jentzen2011
A. Jentzen and P. E. Kloeden.
Taylor Approximations for Stochastic Partial Differential
  Equations, volume 83 of CBMS-NSF Regional Conference Series in Applied
  Mathematics.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia,
  PA, 2011.

jentzen2015
A. Jentzen and M. Röckner.
A Milstein scheme for SPDEs.
Found. Comput. Math., 15(2):313–362, 2015.

kruse2014b
R. Kruse.
Consistency and stability of a Milstein–Galerkin finite element
  scheme for semilinear SPDE.
Stoch. Partial Differ. Equ. Anal. Comput., 2(4):471–516, 2014.

kruse2014
R. Kruse.
Strong and Weak Approximation of Stochastic Evolution
  Equations, volume 2093 of Lecture Notes in Mathematics.
Springer, 2014.

kruse2017
R. Kruse and Y. Wu.
Error analysis of randomized Runge–Kutta methods for differential
  equations with time-irregular coefficients.
Comput. Methods Appl. Math., 17(3):479–498, 2017.

kruse2017b
R. Kruse and Y. Wu.
A randomized Milstein method for stochastic differential equations
  with non-differentiable drift coefficients.
ArXiv preprint, arXiv:1706.09964, 2017.

larsson2009
S. Larsson and V. Thomée.
Partial Differential Equations with Numerical Methods,
  volume 45 of Texts in Applied Mathematics.
Springer-Verlag, Berlin, 2009.
Paperback reprint of the 2003 edition.

leonhard2015
C. Leonhard and A. Rößler.
An efficient derivative-free Milstein scheme for stochastic partial
  differential equations with commutative noise.
ArXiv preprint, arXiv:1509.08427, 2015.

lord2014
G. J. Lord, C. E. Powell, and T. Shardlow.
An Introduction to Computational Stochastic PDEs.
Cambridge Texts in Applied Mathematics. Cambridge University Press,
  New York, 2014.

lord2013
G. J. Lord and A. Tambue.
Stochastic exponential integrators for the finite element
  discretization of SPDEs for multiplicative and additive noise.
IMA J. Numer. Anal., 33(2):515–543, 2013.

pazy1983
A. Pazy.
Semigroups of Linear Operators and Applications to Partial
  Differential Equations, volume 44 of Applied Mathematical Sciences.
Springer, New York, 1983.

rockner2007
C. Prévôt and M. Röckner.
A concise course on stochastic partial differential equations,
  volume 1905 of Lecture Notes in Mathematics.
Springer, 2007.

simon1990
J. Simon.
Sobolev, Besov and Nikolskiĭ fractional spaces: imbeddings
  and comparisons for vector valued spaces on an interval.
Ann. Mat. Pura Appl. (4), 157:117–148, 1990.

spijker1968
M. N. Spijker.
Stability and convergence of finite-difference methods, volume
  1968 of Doctoral dissertation, University of Leiden.
Rijksuniversiteit te Leiden, Leiden, 1968.

spijker1971
M. N. Spijker.
On the structure of error estimates for finite-difference methods.
Numer. Math., 18:73–100, 1971/72.

stummel1973
F. Stummel.
Approximation methods in analysis.
Matematisk Institut, Aarhus Universitet, Aarhus, 1973.
Lectures delivered during the spring term, 1973, Lecture Notes
  Series, No. 35.

thomee2006
V. Thomée.
Galerkin Finite Element Methods for Parabolic Problems,
  volume 25 of Springer Series in Computational Mathematics.
Springer-Verlag, Berlin, second edition, 2006.

wang2015
X. Wang.
An exponential integrator scheme for time discretization of nonlinear
  stochastic wave equation.
J.Sci. Comput., 64(1):234–263, 2015.

wang2017
X. Wang.
Strong convergence rates of the linear implicit Euler method for
  the finite element discretization of SPDEs with additive noise.
IMA J. Numer. Anal., 37(2):965–984, 2017.



