Convergence rates for the numerical approximation of the 2D stochastic Navier–Stokes equations

We study stochastic Navier–Stokes equations in two dimensions with respect to periodic boundary conditions. The equations are perturbed by a nonlinear multiplicative stochastic forcing with linear growth (in the velocity) driven by a cylindrical Wiener process. We establish convergence rates for a finite-element based space-time approximation with respect to convergence in probability (where the error is measured in the Lt∞Lx2∩Lt2Wx1,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^\infty _tL^2_x\cap L^2_tW^{1,2}_x$$\end{document}-norm). Our main result provides linear convergence in space and convergence of order (almost) 1/2 in time. This improves earlier results from Carelli and Prohl (SIAM J Numer Anal 50(5):2467–2496, 2012) where the convergence rate in time is only (almost) 1/4. Our approach is based on a careful analysis of the pressure function using a stochastic pressure decomposition.


Introduction
In this paper we are concerned with the stochastic Navier-Stokes equations ⎧ ⎨ ⎩ du = μΔu dt − (∇u)u dt − ∇π dt + Φ(u)dW in Q, divu = 0 i n Q, and u 0 is a given (random) initial datum. Here the unknowns are the velocity field u : Ω × Q → R N and the pressure π : Ω × Q → R, where Q = (0, T ) × O and O ⊂ R N with N = 2, 3. The stochastic perturbation in the balance of momentum (1.1) 1 can take into account for physical, empirical or numerical uncertainties and thermodynamical fluctuations. In addition to that, a main reason why the stochastic Navier-Stokes equations (1.1) became so popular in fluid mechanical research is their application to turbulence theory (see, for instance, [2,21]). Its mathematical investigation started in the 70's with the pioneering paper of Bensoussan and Temam [1]. They provide a semi-deterministic approach based on the flow-transformation. A first fully stochastic theory has been developed by Flandoli and Gatarek [14] by showing the existence of a martingale solution. These solutions are weak in the stochastic sense meaning that the underlying probability space as well the Wiener process W are not a priori known but become an integral part of the solution. In two dimensions, when uniqueness is known, a stochastically strong solution exists (it is defined on a given probability space with a given Wiener process), see [12]. Nowadays there is a huge amount of literature concerning the analysis of (1.1) and most of the results from the deterministic theory found their stochastic counterpart. For an overview we refer to the recent survey article [25].
The situation about the numerical approximation of (1.1) is totally different and only very few results are available. In [8] a fully practical space-time approximation for the three-dimensional stochastic Navier-Stokes equations in a bounded domain is studied. It is shown that the sequence of approximate solutions converges in law (up to a subsequence) to a martingale solution if both discretization parameters tend to zero. This is the best result one can hope for without using some unproven hypotheses about the space-regularity of solutions (or to be content with local-in-time results). In two dimensions the situation is much better, at least if periodic boundary conditions are considered, that is The space-regularity of the unique strong solution is well-known (see for instance [20]). Based on this the convergence rates for a finite-element based space-time approximation is analysed in [9]. The result is linear convergence in space and convergence of order (almost) 1/4 in time. The precise estimate comparing the solution u and its space-time approximation u h,m reads as for any α < 1/4. Here we have Ω Δt,h ⊂ Ω with P Ω \ Ω Δt,h → 0 as Δt, h → 0. Consequently, one can infer from (1.2) convergence in probability (as defined by Printems [24]) with asymptotic rates (almost) 1/4 and 1 respectively. The aim of the present paper is to improve the convergence rates in time from (almost) 1/4 to (almost) 1/2, see Theorem 3 for the precise statement. This is certainly the optimal convergence rate in time in view of the stochastic forcing and should also be considered as the natural result because of the space-regularity of the solution.
The reason for the low convergence rate in time in [9] is the low time-regularity of the pressure function π in (1.1). Its appearance can only be avoided when working with finite element functions which are exactly divergence-free. Unfortunately, their construction is quite complicated such that the preferred descretizations (such as the Taylor-Hood, the Crouzeix-Raviart, and the MINI element, see [7,16,17]) are only asymptotically divergence-free. The low regularity of the pressure gradient arises from the stochastic forcing (in fact, ∇π behaves as dW ) and improvements do not seem possible unless the noise is divergence-free (this is quite restrictive as it only allows certain additive or linear multiplicative noise). In the general case (non-divergencefree finite elements and nonlinear multiplicative noise) a more subtle analysis of the pressure function is required. Our main idea is a decomposition of the pressure into a deterministic and a stochastic component (such a decomposition first appeared in [5]). The deterministic pressure part behaves as the convective term u ⊗ u. The latter one can be estimated along the lines of [9] (following classical deterministic arguments combined with a discrete stopping time). In addition a second stochastic integral appears which behaves similarly to the original stochastic integral. Although the timeregularity of this part has not improved, we benefit from the averaging properties of the Itô-integral. Combining these ideas finally leads to the optimal convergence rate in Theorem 3.
The paper is organized as follows. In Sect. 2 we present the mathematical framework, that is the probability setup, the concept of solutions and their qualitative properties. In particular, we give improved (compared to [9]) results on the timeregularity of ∇u, see Corollary 2 b). This is needed in Sect. 3 in order to estimate the error between the continuous solution and the time-discrete solution. The heart of the paper is Sect. 4 were we estimate the error between the time-discrete solution and the space-time discretization. Crucial tools are the space-regularity of the time-discrete solution from [8] and the decomposition of the corresponding pressure function.

Probability setup
Let (Ω, F, (F t ) t≥0 , P) be a stochastic basis with a complete, right-continuous filtration. The process W is a cylindrical Wiener process, that is, W (t) = k≥1 β k (t)e k with (β k ) k≥1 being mutually independent real-valued standard Wiener processes relative to (F t ) and (e k ) k≥1 is a complete orthonormal system in a separable Hilbert space U.
To give the precise definition of the diffusion coefficient Φ, consider z ∈ L 2 (T 2 ) and let Φ(z) : U → L 2 (T 2 ) be defined by Φ(z)e k = g k (·, z(·)). In particular, we suppose that g k ∈ C 1 (T 2 × R 2 ) and that the following conditions hold If we are interested in higher regularity some further assumptions are in place and we require additionally We remark that the first inequality of (2.1) implies and (2.1) finally yields In fact, all our results apply if we replace (2.1) and (2.2) by the corresponding norm estimates above. Furthermore, the conditions imposed on Φ, particularly the first assumption from (2.1), allow us to define stochastic integrals Given an (F t )-progressively measurable process u ∈ L 2 (Ω; L 2 (0, T ; L 2 (T 2 ))), the stochastic integral is a well defined process taking values in L 2 (T 2 ) (see [13] for the detailed construction). Moreover, we can multiply by test functions to obtain Similarly, we can define stochastic integrals with values in W 1,2 (T 2 ) and W 2,2 (T 2 ) respectively if u belongs to the corresponding class.

The concept of solutions
In dimension two, pathwise uniqueness for weak solutions is known under (2.1), we refer the reader for instance to Capiński-Cutland [12] and Capiński [11]. Consequently, we may work with the definition of a weak pathwise solution.
, P) be a given stochastic basis with a complete rightcontinuous filtration and an (F t )-cylindrical Wiener process W . Let u 0 be an F 0measurable random variable. Then u is called a weak pathwise solution to (1.1) with the initial condition u 0 provided (a) the velocity field u is (F t )-adapted and holds P-a.s. for all ϕ ∈ C ∞ div (T 2 ) and all t ∈ [0, T ].
Theorem 1 Let N = 2 and assume that Φ satisfies (2.1). Let (Ω, F, (F t ) t≥0 , P) be a given stochastic basis with a complete right-continuous filtration and an (F t )cylindrical Wiener process W . Let u 0 be an F 0 -measurable random variable such that u 0 ∈ L r (Ω; L 2 div (T 2 )) for some r > 2. Then there exists a unique weak pathwise solution to (1.1) in the sense of Definition 1 with the initial condition u 0 . Now, for φ ∈ C ∞ (T 2 ) we can insert φ − ∇Δ −1 divφ and obtain This corresponds to the stochastic pressure decomposition introduced in [5] (see also [6,Chap. 3] for a slightly different presentation). However, the situation with periodic boundary conditions we are considering here is much easier as the harmonic component of the pressure disappears.

Regularity of solutions
Lemma 2 Let the assumptions of Theorem 1 be satisfied. (2.9) Proof Part (a) is the standard a priori estimate which follows from applying Itô's formula to the functional f (u) = 1 2 u 2 L 2 x (and using Burkholder-Davis-Gundy inequality, assumption (2.1) and Gronwall's lemma). Note that this is legit in two dimensions since we have u ⊗ u ∈ L 2 (Q) P-a.s. by Ladyshenskaya's inequality.
The proof of (b) and (c) is quite similar to [20,Corollary 2.4.13]. However, we are working with a different setup. So, we decided to give a formal proof although it is certainly known to experts. The proof can be made rigorous by working with a Galerkin-type approximation and show that the following estimates are uniform with respect to the dimension of the ansatz space. Such a procedure is quite standard, so we leave the details to the reader.
In order to show (b) we apply Itô's formula to the function f γ (u) : (2.10) We take the supremum in time, the r 2 th power and apply expectations. Summing over γ , we find In two dimensions we have (I I ) 3 = 0 by elementary calculations. So we are left with estimating (I I ) 2 and obtain On account of assumption (2.1), Burkholder-Davis-Gundy inequality and Young's inequality we obtain for arbitrary δ > 0 using (2.7) in the last step. By similar arguments we gain Finally, we have by (2.1) Hence we obtain by (2.7) that Plugging all together and choosing δ small enough we have shown (2.8).
The proof of (c) is similar: we simply differentiate once more. We apply Itô's formula to the function f β (u) : where The main difference is that (I I ) 3 does not vanish. By [20, Lemma 2.1.20] (with m = 2) and Young's inequality we have x 5r 2 due to (2.8). By arguments similar to the proof of (b), using (2.2), we gain arguing again similarly to (b) and using (2.2). Again the claim follows by choosing δ small enough.
Under the assumptions of Lemma 2 (b) equation (1.1) is satisfied strongly in the analytical sense. That is we have . In the following we analyze the regularity of the pressure components π det and Φ π . In the following the subscript w * denotes Bochner-measurability with respect to the weak * -topology.

Corollary 1 (a) Under the assumptions of Lemma 2 we have
(c) Under the assumptions of Lemma 2 (c) we have The key tool is the continuity of . So, the regularity transfers from u ⊗ u to π det . We obtain by Ladyshenskaya's inequality using Lemma 2 (a). The estimates for (b) and (c) are very similar. We have to compare ∇π with |u||∇u|, where we have under the assumptions of Lemma 2 (b). Similarly, ∇ 2 π behaves as |u||∇ 2 u| + |∇u| 2 which belongs to the same class since Lemma 2 (c) applies. Now, we investigate the regularity of Φ π . We use continuity Δ −1 to obtain using (2.1) and Lemma 2 (a). Similarly, we have using (2.1) and Lemma 2 (b), as well as 2) and the Lemma 2 (c).
Finally, we investigate the time-regularity of the velocity field.

Finally, we have
by combing Lemma 1 with (2.3). The same conclusion holds for Φ π using Lemma 1 (a). Plugging all together and noting the embedding W 1,2 (0, T ; X ) → C α ([0, T ]; X ) for any separably Banach space X the claim follows. The proof of (b) follows along the same lines using the higher regularity from Lemma 2 (c), Corollary 1 (c) and (2.4).

Discretization in space
We work with a standard finite element set-up for incompressible fluid mechanics, see e.g. [7,15]. We denote by T h a quasi-uniform subdivision of T 2 into triangles of maximal diameter h > 0. For S ⊂ R 2 and ∈ N 0 we denote by P (S) the polynomials on S of degree less than or equal to . Moreover, we set P −1 (S) := {0}. Let us characterize the finite element spaces V h (T 2 ) and P h (T 2 ) as We will assume that i and j are both natural to get (2.15) below (this is different from [9], where also j = 0 is allowed). In order to guarantee stability of our approximations we relate V h (T 2 ) and P h (T 2 ) by the inf-sup condition, that is we assume that where C > 0 does not depend on h. This gives a relation between i and j (for instance the choice (i, j) = (1, 0) is excluded whereas (i, j) = (2, 0) is allowed). Finally, we define the space of discretely solenoidal finite element functions by The following results concerning the approximability of Π h are well-known (see, for instance [18]). There is c > 0 independent of h such that we have for all v ∈ W 1,2 div (T 2 ) and for all p ∈ W 1,2 (T 2 ) and for all p ∈ W 2,2 (T 2 ). Note that (2.17) requires the assumption j ≥ 1 in the definition of P h (T 2 ), whereas (2.16) also holds for j = 0.

Time-discretization
We consider an equidistant partition of [0, T ] with mesh size Δt = T /M and set t m = mΔt. Let u 0 be an F 0 -mesureable random variable with values in W 1,2 div (T 2 ). We aim at constructing iteratively a sequence of F t m -measurable random variables u m with values in W 1,2 div (T 2 ) such that for every φ ∈ W 1,2 div (T 2 ) it holds true P-a.s.
. The existence of a unique u m (given u m−1 and Δ m W ) solving (3.1) is straightforward as it is a stationary Navier-Stokes system. The following result follows from Lemma 3.1 in [9].

Remark 1
A corresponding result is shown in [9, Lemma 3.2] for the full pressure provided the noise is divergence-free (in this case Φ π m vanishes). This is quite restrictive as it means that Φ is linear in u.

Proof
The F t m -measurability of π det m follows directly from the one of u m stated in Lemma 3. By continuity of the operator ∇Δ −1 div on L 2 (T 2 ) we have x which holds in two-dimensions. Now, summing with respect to m, applying expectations and using Lemma 3 yields the claim.
Following [9] we set for ε > 0 using Lemma 3. We obtain the following result.
Theorem 2 Assume that (2.1) holds and that u 0 ∈ L 8 (Ω, W 1,2 div (T 2 )) is an F 0measureable random variable. Let u be the unique strong solution to (1.1) in the sense of Definition 1. Assume that we have for some α ∈ (0, 1 2 ); recall Corollary 2. Let (u m ) M m=1 be the solution to (3.1). Then we have the error estimate for any ε > 0.
for all φ ∈ W 1,2 div (T 2 ). Choosing φ = e m implies Iterating this equality and following the arguments in [9, proof of Thm. 3.1] yields (3.10) Note that only the first bound from (3.7) has been used for (3.10). The second bound from (3.7), not used in [9], allows us to estimates the remaining integral by which finishes the proof.

Finite element based space-time approximation
Now we consider a fully practical scheme combining the implicite Euler scheme in time (as in the last section) with a finite element approximation in space. For a given h > 0 let u h,0 be an F 0 -mesureable random variable with values in V h div (T 2 ) (for instance Π h u 0 ). We aim at constructing iteratively a sequence of random variables u h,m with values in V h div (T 2 ) such that for every φ ∈ V h div (T 2 ) it holds true P-a.s.
. We quote the following result concerning the existence of solutions u h,m to (4.1) from [8, Lemma 3.1].

Error analysis
In this subsection we establish convergence with rates of the above defined algorithm. We introduce for ε > 0 the sample set using Lemma 3 and 6. We also recall the definition of Ω ε Δt in (3.6). We are now ready to state our main result.

Remark 3
We do not expect that it is possible to avoid an indicator function in general (see [24] for the numerical approximation of stochastic PDEs with non-Lipschitz nonlinearities). But it is not clear if our choice of the sample subset is optimal.
The rest of the paper is devoted to the proof of Theorem 3. In fact Theorem 3 will follow from combining Theorem 2 with the following theorem which estimates the error between the time-discrete solution u m to (3.1) and the solution u h,m to (4.1).

Theorem 4
Let u 0 ∈ L 2 (Ω; W 1,2 div (T 2 )) ∩ L 8 (Ω; L 2 div (T 2 )) be F 0 -measurable and assume that Φ satisfies (2.1). Let (u m ) M m=1 be the solution to (3.1). Assume that LΔt ≤ (−ε log h) −1 for some L > 0. Then we have Δt which leads to a restrictive assumption between space and time-discretization. Additionally, we can remove the factor h −3ε arising from the convective term, see (4.6) below. Note that the error term h 2 Δt also appears in [10], where the finite-element based space-time discretization of the stochastic Stokes equations (that is the linearized version of (1.1) without convective term) is studied.

(4.5)
Young's inequality yields for every κ > 0 using also (2.15). The convective term I 2 (m) can be decomposed as As in [9, pages 2489-2491] we introduce the sample set where κ > 0 is arbitrary. For I 1 2 (m), however, we get a slightly better estimate than in [9] since our definition of I 1 2 (m) makes use of divu m = 0. We obtain Now, we explain how to bound line by line in expectation. The error in the initial datum can be bounded by h 2 by (2.14) and the assumption u 0 ∈ L 2 (Ω; W 1,2 div (T 2 )). The second term on the right-hand side can be handled by the discrete Gronwall lemma using the assumption LΔt ≤ (−ε log h) −1 . The expectation of the second line is bounded by h 2 using Lemmas 3 and 4 (the estimate for ∇π det m is the first main ingredient in our proof). This estimate is better than the corresponding one in [9] as only the deterministic part of the pressure appears here. The expectation of the terms in the third line is bounded by This is uniformly bounded as a consequence of Lemmas 3 and 6. The fourth line is controlled by Δt due to the estimate and the uniform bounds from Lemmas 3 and 6 (with q = 3). It remains to estimate the two stochastic terms Note that we used integration by parts and Π h e h,n ∈ V h div (T 2 ) together with the definition of Π π h to rewrite M m,2 into the form above. Finally, we write as well as These representations have the advantage that M 1 m,1 and M 2 m,2 are martingales (note the index n−1 in the indicator functions). Consequently, we can apply the Burkholder-Davis-Gundy inequality to estimate them. As far as M m,1 is concerned we have Here, we also used (2.1) as well as Young's inequality for κ > 0 arbitrary. Now, the first term can be absorbed for κ small enough. The second term can be handled by the discrete Gronwall lemma (note that Ω ε h,n ⊂ Ω ε h,n−1 such that 1 Ω ε h,n ≤ 1 Ω ε h,n−1 P a.s.) Using (2.14) the last term can be estimated by which is bounded by ch 2 using Lemma 3 (recall that u 0 ∈ L 2 (Ω; W 1,2 (T 2 ))). For the term M 2 m,1 we obtain by Cauchy-Schwartz inequality, Young's inequality, Itôisometry and (2.1) The first term can be absorbed for κ small enough whereas the second one can be estimated as for M 1 m,1 . Now, we come to the second main ingredient which is the estimate for M m,2 . We obtain using (2.17), (2.1) and continuity of ∇ 2 Δ −1 The first term is bounded by c κ h 4 using Lemmas 3 (recall that u 0 ∈ L 2 (Ω; W 1,2 (T 2 ))).
The first term can be absorbed for κ small enough. Arguing as for M 1 m,2 the last term is bounded by ch 2 E max n ∇u n−1 2 L 2 x ≤ ch 2 . Plugging all together and noting that Recalling that e h,m = u m − Π h u m + Π h e h,m and using (2.14) as well as Lemma 3 (a) gives the claim.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.