Analysis of Chorin-Type Projection Methods for the Stochastic Stokes Equations with General Multiplicative Noises†

This paper is concerned with developing and analyzing Chorin-type projection finite element methods for the following time-dependent stochastic Stokes problem:

where $D=(0,L)^{2}\subset\mathbb{R}^{d}\,(d=2,3)$ represents a period of the periodic domain in $\mathbb{R}^{d}$, ${\bf u}$ and $p$ stand for respectively the velocity field and the pressure of the fluid, ${\bf B}$ is an operator-valued random field, $\{{\bf W}(t);t\geq 0\}$ denotes an ${\bf L}^{2}(D)$-valued $Q$-Wiener process, and ${\bf f}$ is a body force function (see Section 2 for their precise definitions). Here we seek periodic-in-space solutions $({\bf u},p)$ with period $L$, that is, ${\bf u}(t,{\bf x}+L{\bf e}_{i})={\bf u}(t,{\bf x})$ and $p(t,{\bf x}+L{\bf e}_{i})=p(t,{\bf x})$ almost surely and for any $(t,{\bf x})\in(0,T)\times\mathbb{R}^{d}$ and $1\leq i\leq d$, where $\{\bf e_{i}\}_{i=1}^{d}$ denotes the canonical basis of $\mathbb{R}^{d}$.

The system (1.1a) is a stochastic perturbation of the deterministic Stokes system by introducing a multiplicative noise force term ${\bf B}(\cdot)d{\bf W}(s)$ and it has been used to model turbulent fluids (cf. [1, 2, 17, 21]). The stochastic Stokes system is a simplified model of the full stochastic Navier-Stokes equations by omitting the nonlinear term $({\bf u}\cdot\nabla){\bf u}$ in the drift part of the stochastic Navier-Stokes equations. Although the deterministic Stokes equations is a linear PDE system which has been well studied in the literature (cf. [14, 21] and the references therein), the stochastic Stokes system (1.1a) is intrinsically nonlinear because the diffusion coefficient ${\bf B}$ is nonlinear in the velocity ${\bf u}$. Due to the introduction of random forces it has been well known that the solution of problem (1.1) has very low regularities in time. We refer the reader to [1, 18, 11] and the references therein for a detailed account about the well-posedness and regularities of the solution for system (1.1).

Besides their mathematical and practical importance, the stochastic Stokes (and Navier-Stokes) equations have been used as prototypical stochastic PDEs for developing efficient numerical methods and general numerical analysis techniques for analyzing numerical methods for stochastic PDEs. In that regard several works have been reported in the literature [9, 12, 13, 5, 8, 3]. Euler-Maruyama time dsicretization and divergence-free finite element space discretization was proposed and analyzed in [9] in the case of divergence-free noises (i.e., $\mathbf{B}({\bf u})$ is divergence-free). Optimal order error estimates in strong norm for the velocity approximation were obtained. In [12, 13] the authors considered the general noise and analyzed the standard and a modified mixed finite element methods as well as pressure stabilized methods for space discretization, suboptimal order error estimates were proved in [12] for the velocity approximation in strong norm and for the pressure approximation in a time-averaged norm, all these suboptimal order error estimates were improved to optimal order for a Helmholtz projection-enhanced mixed finite element in [13] (also see [5] for a similar approach). It should be noted that the reason for measuring the pressure errors in a time-averaged norm is because the low regularity of the pressure field which is only a distribution in general and the numerical tests of [12, 13] suggest that these error estimates are sharp. In [8] the authors proposed a Chorin time-splitting finite element method for problem (1.1) and proved a suboptimal convergence rate in strong norm for the velocity approximation in the case of divergence-free noises. In [3] the authors proposed an iterative splitting scheme for stochastic Navier-Stokes equations and a strong convergence in probability was established in the 2-D case for the velocity approximation. In a recent work [4], the authors proposed another time-splitting scheme and proved its strong $L^{2}$ convergence for the velocity approximation.

Compared to the recent advances on mixed finite element methods [9, 12, 13], the numerical analysis of the well-known Chorin projection/splitting scheme for the stochastic Stokes equations lags behind. To the best of our knowledge, the only analysis result obtained in [8] is the optimal convergence in the energy norm for the velocity approximation in the case of divergence-free noises (i.e., $\mathbf{B}({\bf u})$ is divergence-free). Several natural and important questions arise and must be addressed for a better understanding of the Chorin projection scheme for problem (1.1). Among them are (i) Does the pressure approximation converge even when the noise is divergence-free? If so, in what sense and what rate? (ii) Does the Chorin projection scheme converge (for both the velocity and pressure approximations) for general noises? If so, in what sense and what rate? (iii) Could the performance of the standard Chorin projection scheme be improved one way or another in the case of general noises? The primary objective this paper is to provide a positive answer to each of the above questions.

As it was shown in [8], the adaptation of the standard deterministic Chorin projection scheme to problem (1.1) is straightforward (see Algorithm 1 of Section 3). The idea of the Chorin scheme is to separate the computation of the velocity and pressure at each time step which is done by solving two decoupled Poisson problems and the divergence-free constraint for the velocity approximation is enforced by a Helmholtz projection technique which can be easily obtained using the solutions of the two Poisson problems. The Chorin scheme also can be compactly rewritten as a pressure stabilization scheme at each time step as follows (cf. [8]):

One of advantages of the above Chorin scheme is that the spatial approximation spaces for $\tilde{{\bf u}}^{n+1}$ and $p^{n+1}$ can be chosen independently, so unlike in the mixed finite element method, they are not required to satisfy an inf-sup condition. Notice that a time lag on pressure appears in equation (1.2a) which causes most of difficulties in the convergence analysis (cf. [20, 15, 19, 8]). We also note that the term $-k\Delta p^{n+1}$ in equation (1.2b) is known as a pressure stabilization term.

To improve the convergence of the standard Chorin scheme, we adopt a Helmholtz projection technique as used in [13] (also see [5]). At each time step we first perform the Helmholtz decomposition ${\bf B}(\tilde{{\bf u}}^{n})=\boldsymbol{\eta}^{n}+\nabla\xi^{n}$ and then rewrite (1.2a) as

where $r^{n}=p^{n}-k^{-1}\xi^{n}\Delta W_{n+1}$. Our modified Chorin scheme consists of (1.3), (1.2b)–(1.2c) and the Helmholtz decomposition ${\bf B}(\tilde{{\bf u}}^{n})=\boldsymbol{\eta}^{n}+\nabla\xi^{n}$. Since $\boldsymbol{\eta}^{n}$ is divergence-free, it turns out that the finite element approximation of the modified Chorin scheme has better convergence properties. Notice that $p^{n}$ can be recovered from $r^{n}$ via the simple algebraic relation $p^{n}=r^{n}+k^{-1}\xi^{n}\Delta W_{n+1}$.

The main contributions of this paper are summarized below.

• •

  We proved the following error estimates in strong norms for the Chorin-$P_{1}$ finite element method (see Algorithm 3) for problem (1.1) with general multiplicative noises:

  	$\displaystyle\biggl{(}\mathbb{E}\Bigl{[}k\sum_{m=0}^{M}\|{\bf u}(t_{m})-\tilde{{\bf u}}^{m}_{h}\|^{2}\Bigr{]}\biggr{)}^{\frac{1}{2}}$	$\displaystyle+\max_{0\leq\ell\leq M}\biggl{(}\mathbb{E}\Bigl{[}\Bigl{\|}k\sum_{m=0}^{\ell}\nabla({\bf u}(t_{m})-\tilde{{\bf u}}^{m}_{h})\Bigr{\|}^{2}\Bigr{]}\biggr{)}^{\frac{1}{2}}$
  		$\displaystyle\qquad\qquad\qquad\leq C\Bigl{(}k^{\frac{1}{4}}+hk^{-\frac{1}{2}}\Bigr{)},$
  	$\displaystyle\biggl{(}\mathbb{E}\Bigl{[}k\sum_{m=0}^{M}\Bigl{\|}P(t_{m})$	$\displaystyle-k\sum_{n=0}^{m}p^{n}_{h}\Bigr{\|}^{2}\Bigr{]}\biggr{)}^{\frac{1}{2}}\leq C\Bigl{(}k^{\frac{1}{4}}+hk^{-\frac{1}{2}}\Bigr{)},$

  where $({\bf u}(t_{m}),P(t_{m}))$ are the solution to problem (1.1) while $(\tilde{{\bf u}}_{h}^{m},p_{h}^{m})$ are the discrete solution of Algorithm 3, see Sections 2 and 4 for their precise definitions.

• •

  We proposed a modified Chorin-$P_{1}$ finite element method (see Algorithm 4) and proved the following error estimates in strong norms for problem (1.1) with general multiplicative noises:

  	$\displaystyle\max_{1\leq m\leq M}\Bigl{(}\mathbb{E}\bigl{[}\bigl{\|}{\bf u}(t_{m})-\tilde{{\bf u}}^{m}_{h}\bigr{\|}^{2}\,\bigl{]}\Bigr{)}^{\frac{1}{2}}+\biggl{(}\mathbb{E}\biggl{[}k\sum_{m=1}^{M}\bigl{\|}\nabla({\bf u}(t_{n})-{\bf u}^{n}_{h})\bigr{\|}^{2}\,\biggr{]}\biggr{)}^{\frac{1}{2}}$
  	$\displaystyle\hskip 137.31255pt\leq C\Bigl{(}k^{\frac{1}{2}}+h+k^{-\frac{1}{2}}h^{2}\Bigr{)},$
  	$\displaystyle\biggl{(}\mathbb{E}\biggl{[}\bigg{\|}R(t_{m})-k\sum^{m}_{n=1}r^{n}_{h}\biggr{\|}^{2}\,\biggr{]}\biggr{)}^{\frac{1}{2}}+\biggl{(}\mathbb{E}\biggl{[}\bigg{\|}P(t_{m})-k\sum^{m}_{n=1}p^{n}_{h}\biggr{\|}^{2}\,\biggr{]}\biggr{)}^{\frac{1}{2}}$
  	$\displaystyle\hskip 137.31255pt\leq C\Bigl{(}k^{\frac{1}{2}}+h+k^{-\frac{1}{2}}h^{2}\Bigr{)}.$

  where $({\bf u}(t_{m}),P(t_{m}))$ is the solution to problem (1.1) and $R(t)$ is defined as the time-average of the pseudo pressure $r(t)$ while $({\bf u}^{m}_{h},r^{m}_{h},p^{m}_{h})$ is the solution of Algorithm 4, see Sections 2 and 4 for their precise definitions.

We note that all spatial error constants contain a growth factor $k^{-\frac{1}{2}}$, which explains the deteriorating performance of the standard (and modified) Chorin scheme when $k\to 0$ and the mesh size $h$ is fixed as observed in the numerical tests of [8]. The numerical experiments to be given in Section 5 indicate that the dependence on factor $k^{-\frac{1}{2}}$ is sharp.

The remainder of this paper is organized as follows. In Section 2, we first introduce some space notations and state the assumptions on the initial data and on ${\bf B}$ as well as recall the definition of solutions to (1.1). We then state and prove a Hölder continuity property for the pressure $p$ in a time-averaged norm. In Section 3, we define the standard Chorin projection scheme as Algorithm 1 for problem (1.1) in Subsection 3.1 and the modified Chorin scheme as Algorithm 2 in Subsection 3.2. The highlights of this section are to prove some uniform (in $k$) stability estimates which are very useful for error analysis later. In Section 4, we formulate the finite element spatial discretization for both the standard Chorin and modified Chorin schemes in Algorithm 3 and 4, respectively and prove the quasi-optimal error estimates for both algorithms as summarized above. In Section 5, we present several numerical experiments to gauge the performance of the proposed numerical methods and to validate the sharpness of the proved error estimates.

Standard function and space notation will be adopted in this paper. Let ${\bf H}^{1}_{0}(D)$ denote the subspace of ${\bf H}^{1}(D)$ whose ${\mathbb{R}}^{d}$-valued functions have zero trace on $\partial D$, and $(\cdot,\cdot):=(\cdot,\cdot)_{D}$ denote the standard $L^{2}$-inner product, with induced norm $\|\cdot\|$. We also denote ${\bf L}^{p}_{per}(D)$ and ${\bf H}^{k}_{per}(D)$ as the Lebesgue and Sobolev spaces of the functions that are periodic and have vanishing mean, respectively. Let $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\},\mathbb{P})$ be a filtered probability space with the probability measure $\mathbb{P}$, the $\sigma$-algebra $\mathcal{F}$ and the continuous filtration $\{\mathcal{F}_{t}\}\subset\mathcal{F}$. For a random variable $v$ defined on $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\},\mathbb{P})$, ${\mathbb{E}}[v]$ denotes the expected value of $v$. For a vector space $X$ with norm $\|\cdot\|_{X}$, and $1\leq p<\infty$, we define the Bochner space $\bigl{(}L^{p}(\Omega,X);\|v\|_{L^{p}(\Omega,X)}\bigr{)}$, where $\|v\|_{L^{p}(\Omega,X)}:=\bigl{(}{\mathbb{E}}[\|v\|_{X}^{p}]\bigr{)}^{\frac{1}{p}}$. We also define

We recall from [14] that the (orthogonal) Helmholtz projection ${\bf P}_{{\mathbb{H}}}:{\bf L}^{2}_{per}(D)\rightarrow{\mathbb{H}}$ is defined by ${\bf P}_{{\mathbb{H}}}{\bf v}=\boldsymbol{\eta}$ for every ${\bf v}\in{\bf L}^{2}_{per}(D)$, where $(\boldsymbol{\eta},\xi)\in{\mathbb{H}}\times H^{1}_{per}(D)/\mathbb{R}$ is a unique tuple such that

and $\xi\in H^{1}_{per}(D)/\mathbb{R}$ solves the following Poisson problem with the homogeneous Neumann boundary condition:

We also define the Stokes operator ${\bf A}:=-{\bf P}_{\mathbb{H}}\Delta:{\mathbb{V}}\cap{\bf H}^{2}_{per}(D)\rightarrow{\mathbb{H}}$.

Throughout this paper we assume that ${\bf B}:{\bf L}^{2}_{per}(D)\rightarrow{\bf L}^{2}_{per}(D)$ is a Lipschitz continuous mapping and has linear growth, that is, there exists a constant $C>0$ such that for all ${\bf v},{\bf w}\in{\bf L}^{2}_{per}(D)$

Since we shall not explicitly track the dependence of all constants on $\nu$, for ease of the presentation, unless it is stated otherwise, we shall set $\nu=1$ in the rest of the paper and assume that ${\bf f}\in L^{2}(\Omega;L^{2}_{per}(D))$. In addition, we shall use $C$ to denote a generic positive constant which may depend on $T$, the datum functions ${\bf u}_{0}$ and ${\bf f}$, and the domain $D$ but is independent of the mesh parameter $h$ and $k$.

In this section, we first formulate two Chorin-type semi-discrete-in-time schemes for problem (1.1). The first scheme is the standard Chorin scheme and the second one is a Helmholtz decomposition enhanced nonstandard Chorin scheme. We then present a complete convergence analysis for each scheme which include establishing their stability and error estimates in strong norms for both velocity and pressure approximations.