Fully Discrete Mixed Finite Element Methods for the Stochastic Cahn-Hilliard Equation with Gradient-type Multiplicative Noise

TL;DR

This paper introduces and analyzes fully discrete mixed finite element methods for the stochastic Cahn-Hilliard equation with gradient-type multiplicative noise, proving strong convergence and optimal rates, supported by numerical experiments.

Contribution

The paper develops new finite element methods for the stochastic Cahn-Hilliard equation and establishes their strong convergence with optimal rates, addressing low regularity challenges.

Findings

Proved strong convergence with optimal rates.

Established Hölder continuity in time for solutions.

Validated results through numerical experiments.

Abstract

This paper develops and analyzes some fully discrete mixed finite element methods for the stochastic Cahn-Hilliard equation with gradient-type multiplicative noise that is white in time and correlated in space. The stochastic Cahn-Hilliard equation is formally derived as a phase field formulation of the stochastically perturbed Hele-Shaw flow. The main result of this paper is to prove strong convergence with optimal rates for the proposed mixed finite element methods. To overcome the difficulty caused by the low regularity in time of the solution to the stochastic Cahn-Hilliard equation, the H\"{o}lder continuity in time with respect to various norms for the stochastic PDE solution is established, and it plays a crucial role in the error analysis. Numerical experiments are also provided to validate the theoretical results and to study the impact of noise on the Hele-Shaw flow as well as the interplay of the geometric evolution and gradient-type noise.

Tables

Table 1. Table 1: Test 1: Temporal errors and convergence rates with ϵ = 0.1 italic-ϵ 0.1 \epsilon=0.1 , δ = 5 𝛿 5 \delta=5 and h = 2 / 2 6 ℎ 2 superscript 2 6 h=2/2^{6} .

$\tau$	$L^{\infty }\left(L^2\right)$ error	order	$L^2\left(H^1\right)$ error	order
3.2000E-03	3.93E-02		1.92E-02
1.6000E-03	1.52E-02	1.37	8.25E-03	1.22
8.0000E-04	9.56E-03	0.67	5.01E-03	0.72
4.0000E-04	6.59E-03	0.54	3.27E-03	0.61
2.0000E-04	4.77E-03	0.47	2.27E-03	0.52
1.0000E-04	3.28E-03	0.54	1.59E-03	0.52

Table 2. Table 2: Test 1: Spatial errors and convergence rates with ϵ = 0.1 italic-ϵ 0.1 \epsilon=0.1 , δ = 25 𝛿 25 \delta=25 , τ = 5 × 10 − 5 𝜏 5 superscript 10 5 \tau=5\times 10^{-5} and M = 10 4 𝑀 superscript 10 4 M=10^{4} .

$h$	$L^{\infty }\left(L^2\right)$ error	order	$L^2\left(H^1\right)$ error	order
5.0000E-01	9.30E-02		3.10E-02
2.5000E-01	3.47E-02	1.42	1.94E-02	0.67
1.2500E-01	1.10E-02	1.66	9.38E-03	1.05
6.2500E-02	3.66E-03	1.58	4.63E-03	1.02

Equations

$$\displaystyle du=\Bigl{[}-\Delta\Big{(}\epsilon\Delta u-\frac{1}{\epsilon}f(u)\Big{)}\Bigr{]}dt+\delta\nabla u\cdot X\circ dW_{t}$$

$$\displaystyle\frac{\partial u}{\partial n}=\frac{\partial}{\partial n}\Big{(}\epsilon\Delta u-\frac{1}{\epsilon}f(u)\Big{)}=0$$

$$u=u_0$$

$$F(u)=\frac{1}{4}(u^2-1{)}^2.$$

$$\displaystyle du=\left[-\Delta\Big{(}\epsilon\Delta u-\frac{1}{\epsilon}f(u)\Big{)}+\frac{\delta^{2}}{2}\text{div}(B\nabla u)\right]dt+\delta\nabla u\cdot XdW_{t},$$

$$\displaystyle du=\Bigl{[}\Delta w+\frac{\delta^{2}}{2}\text{div}(B\nabla u)\Bigr{]}dt+\delta\nabla u\cdot XdW_{t}$$

$$w=-ϵ\Delta u+\frac{1}{ϵ}f(u)$$

$$\frac{\partial u}{\partial n}=\frac{\partial w}{\partial n}=0$$

$$u=u_0$$

$$\Delta w$$

$$\frac{\partial w}{\partial n}$$

$$w$$

$$V_n$$

$${\Gamma }_0$$

$$V_n=\frac{1}{2}{\left[\frac{\partial w}{\partial n}\right]}_{{\Gamma }_t}+\delta {\stackrel{\circ }{W}}_{t}X\cdot n,$$

$$u_t=\Delta u-\frac{1}{{\epsilon }^2}f(u),$$

$$\underset{t\in [0,T]}{\sup }\mathbb{E}\left[\Vert u{\Vert }_{W^{s,\infty }(\mathcal{D})}^p\right]\le C_0=C(p,\delta ,ϵ)\forall p\ge 1,$$

$$(u(t),\varphi )$$

$$+\delta {\int }_0^t(\nabla u(s)\cdot Xd{W}_s,\varphi )\forall \varphi \in H^1(\mathcal{D})\forall t\in (0,T],$$

$$(w(t),\phi )$$

$$\mathbb{E}\left[\Vert \nabla (u(s)-u(t)){\Vert }_{L^2(\mathcal{D})}^2\right]+ϵ\mathbb{E}\left[{\int }_t^s\Vert \nabla \Delta (u(\zeta )-u(t)){\Vert }_{L^2(\mathcal{D})}^{2}d\zeta \right]\le C_1(s-t),$$

$$C_1=C\underset{t\le \zeta \le s}{\sup }\mathbb{E}\left[\Vert \nabla \Delta u(\zeta ){\Vert }_{L^2(\mathcal{D})}^2\right]+C(\delta ,\frac{1}{ϵ})\underset{t\le \zeta \le s}{\sup }\mathbb{E}\left[\Vert u(\zeta ){\Vert }_{H^2(\mathcal{D})}^4\right].$$

$${\Psi }^{\prime }(u)(v)$$

$${\Psi }^{\prime \prime }(u)(m,v)$$

$$\displaystyle\|\nabla(u(s)-u(t))\|_{L^{2}(\mathcal{D})}^{2}=2\int_{t}^{s}\bigl{(}\nabla(u(\zeta)-u(t)),\nabla(-\Delta(\epsilon\Delta u(\zeta)-\frac{1}{\epsilon}f(u(\zeta)))$$

$$\displaystyle\qquad+\frac{\delta^{2}}{2}\mbox{div}(B\nabla u(\zeta)))\bigr{)}d\zeta+2\int_{t}^{s}\bigl{(}\nabla(u(\zeta)-u(t)),\nabla(\delta\nabla u(\zeta)\cdot Xd{W}_{\zeta})\bigr{)}$$

$$+{\delta }^2{\int }_t^s(\nabla (\nabla u(\zeta )\cdot X),\nabla (\nabla u(\zeta )\cdot X))d\zeta .$$

$$\displaystyle\|\nabla(u(s)-u(t))\|_{L^{2}(\mathcal{D})}^{2}=2\int_{t}^{s}\bigl{(}\nabla\Delta(u(\zeta)-u(t)),-\epsilon\nabla\Delta(u(\zeta)-u(t))\bigr{)}d\zeta$$

$$\displaystyle\qquad-2\int_{t}^{s}\bigl{(}\nabla\Delta(u(\zeta)-u(t)),\epsilon\nabla\Delta u(t)\bigr{)}d\zeta$$

$$\displaystyle\qquad+2\int_{t}^{s}\bigl{(}\nabla\Delta(u(\zeta)-u(t)),\frac{1}{\epsilon}\nabla f(u(\zeta))\bigr{)}d\zeta$$

$$\displaystyle\qquad-\delta^{2}\int_{t}^{s}\bigl{(}\Delta(u(\zeta)-u(t)),B:D^{2}u(\zeta)+\nabla u(\zeta)\cdot\mbox{div}(B)\bigr{)}d\zeta$$

$$\displaystyle\qquad+2\int_{t}^{s}\bigl{(}\nabla(u(\zeta)-u(t)),\nabla(\delta\nabla u(\zeta)\cdot Xd{W}_{\zeta})\bigr{)}$$

$$+{\delta }^2{\int }_t^s{\int }_{\mathcal{D}}|D^{2}u(\zeta )X+(\nabla X{)}^T\nabla u(\zeta ){|}^{2}dxd\zeta ,$$

$$\mathbb{E}\left[\Vert \nabla w(s)-\nabla w(t){\Vert }_{L^2(\mathcal{D})}^2\right]\le C_2(s-t),$$

$$C_2=C\underset{t\le \zeta \le s}{\sup }\mathbb{E}\left[\Vert u(\zeta ){\Vert }_{H^7(\mathcal{D})}^2\right]+C(\delta ,\frac{1}{ϵ})\underset{t\le \zeta \le s}{\sup }\mathbb{E}\left[\Vert u(\zeta ){\Vert }_{H^6(\mathcal{D})}^6\right].$$

$$g_1(u(s))$$

$$g_2(u(s))$$

$$g_1^{\prime }(u)(v)$$

$$g_1^{\prime \prime }(u)(m,v)$$

$$g_2^{\prime }(u)(v)$$

$$\displaystyle\qquad\cdot\bigl{[}6u(s)v(s)\nabla u(s)+3u^{2}(s)\nabla v(s)-\nabla v(s)\bigr{]}dx,$$

Full text

{AMS}

65N12, 65N15, 65N30

1 Introduction

We consider the following stochastic Cahn-Hilliard (SCH) problem:

[TABLE]

where $\mathcal{D}\subset\mathbb{R}^{d}$ ($d=2,3$) is a bounded domain, $n$ stands for the unit outward normal to $\partial\mathcal{D}$, and $T>0$ is a fixed number. $W_{t}$ denotes a standard real-valued Wiener process on a given filtered probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}:t\geq 0\},\mathbb{P})$, “$\circ$” refers to the Stratonovich interpretation of the stochastic integral. $X:\mathbb{R}^{d}\longrightarrow\mathbb{R}^{d}$ is a smooth divergence-free vector field defined on $\mathcal{D}$ satisfying $X\cdot n=0$ on $\partial\mathcal{D}$.

Moreover, $f=F^{\prime}$, the derivative of a smooth double equal well potential $F$ taking its global minimum zero at $\pm 1$. In this paper we focus on the following quartic potential density function:

[TABLE]

We note that the Stratonovich SPDE (1) can be equivalently rewritten as the following Itô SPDE:

[TABLE]

where $B=X\otimes X\in\mathbb{R}^{d\times d}$ with $B_{ij}=X_{i}X_{j}$ ($i,j=1,...,d$).

By introducing the so-called chemical potential $w:=-\epsilon\Delta u+\frac{1}{\epsilon}f(u)$, the above primal formulation of the SCH problem can be rewritten as the following mixed formulation:

[TABLE]

which will be used to develop fully discrete finite element numerical methods in this paper.

The deterministic Cahn-Hilliard equation (i.e., $\delta=0$) was originally introduced in [7] to describe complicated phase separation and coarsening phenomena in a melted alloy that is quenched to a temperature at which only two different concentration phases can exist stably. In the equation, $u$ represents the concentration of one of two metallic components of the alloy mixture, the small parameter $\epsilon>0$ is called the interaction length. Note that in (6)–(7), $t$ is the fast time representing $\frac{t}{\epsilon}$ in the original Cahn-Hilliard formulation. The existence of bistable states suggests that nonconvex energy is associated with the equation (cf. [2, 11, 7]). The Cahn-Hilliard equation is well-known also because it closely relates to a celebrated moving interface problem, namely the Hele-Shaw (or Mullins-Sekerka) problem/flow. It was proved in [37, 2] that, as $\epsilon\searrow 0$, the chemical potential $w:=-\epsilon\Delta u+\epsilon^{-1}f(u)$ tends to a limit, which, together with a free boundary $\Gamma:=\cup_{0\leq t\leq T}(\Gamma_{t}\times\{t\})$, satisfies the following Hele-Shaw (or Mullins-Sekerka) problem:

[TABLE]

where $\sigma=\int_{-1}^{1}\sqrt{\frac{F(s)}{2}}\,ds$, $\kappa$ and $V_{n}$ are the mean curvature and the outward normal velocity of the interface $\Gamma_{t}$, $n$ is the unit outward normal to either $\partial\mathcal{D}$ or $\Gamma_{t}$, $\left[\frac{\partial w}{\partial n}\right]_{\Gamma_{t}}:=\frac{\partial w^{+}}{\partial n}-\frac{\partial w^{-}}{\partial n}$, and $w^{+}$ and $w^{-}$ are respectively the restriction of $w$ in the exterior and interior of $\Gamma_{t}$ in $\mathcal{D}$. More details about the justification of the limit can be found in [2, 8, 40] and its numerical approximations in [13, 14, 15, 29, 41] and in [11].

In applications of the Hele-Shaw flow, uncertainty may arise and come from various sources such as thermal fluctuation, impurities of the materials and the intrinsic instabilities of the deterministic evolutions. Therefore, the evolution of the flow/interface under influence of noise is of great importance in applications, it is necessary and interesting to consider stochastic effects, and to study the impact of noise on its phase field models and solutions, especially on their long time behaviors. This then leads to considering the stochastic phase field models. However, how to incorporate noises correctly into phase field models is often a delicate issue.

In this paper, we consider the following stochastically perturbed Hale-Shaw flow:

[TABLE]

where a white-in-time noise multiplied by a smooth spatial coefficient function $X$ is added to the normal velocity of the interface $\Gamma_{t}$, and the parameter $\delta>0$ represents the noise intensity. By an heuristic argument (see [39, 11] for an analogous argument), we can formally show that equation (1) is a phase field formulation of the above stochastic Hele-Shaw flow.

It should be noted that there is another stochastic Cahn-Hilliard equation, called Cahn-Hilliard-Cook (CHC) equation, which has been extensively studied in the literature, see [10, 34, 5] for PDE analysis and [28, 25, 26, 19] and the references therein for its numerical approximations. However, the noise in the CHC equation is additive and the parameter $\epsilon=1$ in those works. Hence, there may have no connection between the CHC equation and the above stochastic Hele-Shaw flow. We also note that numerical approximations of various stochastic versions of the following Allen-Cahn equation:

[TABLE]

which is a closely related to the Cahn-Hilliard equation, have been extensively investigated in the literature [23, 21, 30, 24, 32, 16]. Most of those works focused on either additive noise or function-type multiplicative noise. Recently, finite element approximations of the stochastic Allen-Cahn (SAC) equation with gradient-type multiplicative noise had been carried out by the authors in [17]. This SAC equation was derived as and partially proved to be a phase field formulation of the stochastic mean curvature flow [22, 39, 42].

The goal of this paper is to extend the work of [17] to the stochastic Cahn-Hilliard problem (1)–(3). Specifically, we shall develop and analyze a fully discrete mixed finite element method for this problem, and establish strong convergence with rates for the proposed mixed finite element method under the assumption that the strong solution $u(\cdot,t)$ of the underlying SPDE problem belongs to $W^{s,\infty}(\mathcal{D})$ for almost every $t\in[0,T]$ and it satisfies the following high moment estimates:

[TABLE]

where $\mathbb{E}\left[{\cdot}\right]$ denotes the expectation operator. It turns out that the divergence-free property of $X$ plays a key role in our analysis, which guarantees the sample-wise mass conservation for the strong solution to problem (1)–(3). Another key ingredient for the error analysis is the Hölder continuity estimates for the strong solution. To the best of our knowledge, numerical analysis has yet been done for the stochastic Cahn-Hilliard equation with gradient-type multiplicative noise in the literature.

The rest of the paper is organized as follows. In Section 2, we define the weak formulation for problem (1)–(3) and derive several Hölder continuity estimates for the strong solution of the SPDE problem. In Section 3, a fully discretized mixed finite element method is formulated and properties of the discrete inverse Laplacian operator are presented, which will be utilized to establish the well-posedness and stability of the discrete method, and to prove the strong convergence with rates in Section 4. Finally, in Section 5, we report several numerical experiments to validate our theoretical results and to examine the interplay of the geometric parameter $\epsilon$ and the noise intensity $\delta$.

Throughout this paper we shall use $C$ to denote a generic positive constant independent of the parameters $\epsilon$, $\delta$, space and time mesh sizes $h$ and $\tau$, which can take different values at different occurrences.

2 Preliminaries

Standard functional space and function notation in [1, 6] will be adopted in this paper. In particular, $H^{k}(\mathcal{D})$ for $k\geq 0$ denotes the Sobolev space of order $k$, $(\cdot,\cdot)$ and $\|\cdot\|_{L^{2}(\mathcal{D})}$ denote the standard inner product and norm of $L^{2}(\mathcal{D})$.

In this section, we shall establish several technical lemmas about Hölder continuity estimates for the strong solution of problem (1)–(3) that play a key role in error analysis in Section 4. These estimates play the role of the time derivatives of the solution in the deterministic case.

First, we define the weak formulation for problem (1)–(3), based on the mixed formulation (6)–(9), as follows: Seeking an $\mathcal{F}_{t}$-adapted and $H^{1}(\mathcal{D})\times H^{1}(\mathcal{D})$-valued process $(u(\cdot,t),w(\cdot,t))$ such that there hold $\mathbb{P}$-almost surely

[TABLE]

Second, we derive a Hölder continuity estimate in time of the solution function $u$ with respect to the spatial $H^{1}$-seminorm.

Lemma 2.1**.**

Let $(u,w)$ be the solution to problem (6)-(9) and assume $u$ is sufficiently regular in the spatial variable. Then for any $t,s\in[0,T]$ with $t<s$, we have

[TABLE]

where

[TABLE]

Proof 2.2**.**

Apply Itô’s formula to $\Psi(u(s)):=||\nabla u(s)-\nabla u(t)||_{L^{2}(\mathcal{D})}^{2}$, and notice that

[TABLE]

then we have

[TABLE]

Then we obtain

[TABLE]

where the embedding theorem from $H^{2}(\Omega)$ to $L^{\infty}(\Omega)$ is used in estimating the nonlinear term.

*Taking the expectation on both sides of (22), and using Young’s inequality, the lemma is proved. *

It turns out we also need to control the chemical potential $w$ to handle the nonlinear terms in the error analysis. The following lemma establishes a Hölder continuity estimate in time for $w$ with respect to the spatial $H^{1}$-seminorm.

Lemma 2.3**.**

Let $(u,w)$ be the solution to problem (6)-(9) which is assumed to be sufficiently regular in the spatial variable. Then for any $t,s\in[0,T]$ with $t<s$, we have

[TABLE]

where

[TABLE]

Proof 2.4**.**

Define $g(u(s)):=g_{1}(u(s))+g_{2}(u(s))$, where

[TABLE]

Notice that

[TABLE]

and

[TABLE]

Applying Itô’s formula to $g(w(s)):=\|\nabla w(s)-\nabla w(t)\|_{L^{2}(\mathcal{D})}^{2}$, then we have

[TABLE]

where

[TABLE]

*Taking the expectation on both sides of (27), and using Young’s inequality and the embedding theorem, the lemma is proved. *

3 Formulation of mixed finite element method

In this section we define our mixed finite element method for (17)-(18) and introduce several auxiliary operators that will be used in Section 4.

Let $t_{n}=n\tau$ ($n=0,1,...,N$) be a uniform partition of $[0,T]$ with $\tau=T/N$ and $\mathcal{T}_{h}$ be a quasi-uniform triangulation of $\mathcal{D}$. Let $V_{h}$ be the finite element space given by

[TABLE]

where $\mathcal{P}_{1}(K)$ denotes the space of polynomials of degree one on $K\in\mathcal{T}_{h}$. Our fully discrete mixed finite element methods for (17)-(18) is defined as seeking $\mathcal{F}_{t_{n}}$-adapted and $V_{h}\times V_{h}$-valued process $\{(u_{h}^{n},w_{h}^{n})\}$ ($n=1,\dots,N$) such that $\mathbb{P}$-almost surely

[TABLE]

where $\bar{\Delta}$ denotes the difference operator, $\bar{\Delta}W_{n+1}:=W_{t_{n+1}}-W_{t_{n}}\sim\mathcal{N}(0,\tau)$ and $f^{n+1}:=(u^{n+1}_{h})^{3}-u^{n+1}_{h}$. The initial values $(u^{0}_{h},w^{0}_{h})$ are chosen by solving

[TABLE]

Note that $u^{0}_{h}=P_{h}u_{0}$ where $P_{h}:L^{2}(\mathcal{D})\longrightarrow V_{h}$ is the standard $L^{2}$-projection operator satisfying the following error estimates [9, 6]

[TABLE]

for all $v\in H^{2}(\mathcal{D})$. Furthermore, a direct calculation shows that the numerical solution function $u^{n}_{h}$ satisfies the sample-wise mass conservation property, i.e., $(u^{n}_{h},1)=(u_{0},1)$ almost surely for all $n=0,1,...,N$.

Let $\mathring{V}_{h}$ be the subspace of $V_{h}$ with zero mean, i.e.,

[TABLE]

We introduce the inverse discrete Laplace operator $\Delta_{h}^{-1}:\mathring{V}_{h}\rightarrow\mathring{V}_{h}$ as follows: given $\zeta\in\mathring{V}_{h}$, define $\Delta_{h}^{-1}\zeta\in\mathring{V}_{h}$ such that

[TABLE]

For any $\zeta,\Phi\in\mathring{V}_{h}$, we can define the discrete $H^{-1}$ inner product by

[TABLE]

The induced mesh-dependent $H^{-1}$ norm is given by

[TABLE]

The following properties can be easily verified (cf. [3]):

[TABLE]

and, if $\mathcal{T}_{h}$ is quasi-uniform, we further have

[TABLE]

Setting $\hat{u}^{n}_{h}=u^{n}_{h}-\bar{u}_{0}$ and $\hat{w}^{n}_{h}=w^{n}-\bar{w}^{n}_{h}$, where $\bar{v}=|\mathcal{D}|^{-1}(v,1)$, we can equivalently formulate (28)–(29) as: seeking $\mathcal{F}_{t_{n}}$-adapted and $\mathring{V}_{h}\times\mathring{V}_{h}$-valued process $\{(\hat{u}_{h}^{n},\hat{w}_{h}^{n})\}$ ($n=1,\dots,N$) such that $\mathbb{P}$-almost surely

[TABLE]

where $\hat{f}^{n+1}=(\hat{u}^{n+1}_{h}+\bar{u}_{0})^{3}-(\hat{u}^{n+1}_{h}+\bar{u}_{0})$.

The next theorem establishes the well-posedness for the proposed numerical method.

Theorem 3.1**.**

The scheme (28)–(29) (or (39)–(40)) is uniquely solvable, provided that the following mesh constraint is satisfied

[TABLE]

Proof 3.2**.**

For any $v_{h}\in\mathring{V}_{h}$, let $\eta_{h}=-\Delta^{-1}_{h}v_{h}\in\mathring{V}_{h}$ in (39), we have

[TABLE]

By (34), (40) and integration by parts, we can rewrite (42) as

[TABLE]

Now we define $B:\mathring{V}_{h}\longrightarrow\mathring{V}_{h}$ by

[TABLE]

For any $v_{h}\in\mathring{V}_{h}$, we have

[TABLE]

Notice that

[TABLE]

where we had used

[TABLE]

to obtain the last inequality. Moreover, we have

[TABLE]

by (36), and

[TABLE]

by the Cauchy-Schwarz inequality. Combining (45)–(48) yields

[TABLE]

Hence we have

[TABLE]

provided that the mesh constraint (41) holds. It follows from Brouwer’s fixed point theorem (cf. [20] and [38, Lemma 7.2]) that there exists $\hat{u}^{n+1}_{h}\in\mathring{V}_{h}$ such that

[TABLE]

which also implies the existence of the solution to (43). This $\hat{u}^{n+1}_{h}$ together with $\hat{w}^{n+1}_{h}$ determined by (40) solves (39)–(40).

Next, it suffices to establish the uniqueness of the solution to (43). Assume $\hat{u}^{n+1}_{h,1}$ and $\hat{u}^{n+1}_{h,2}$ are two solutions to (43). Denote $U^{n+1}_{h}=\hat{u}^{n+1}_{h,1}-\hat{u}^{n+1}_{h,2}$, we have

[TABLE]

Taking $v_{h}=U^{n+1}_{h}$ in (53), using the fact that

[TABLE]

and the similar estimates to (47) and (48), we obtain

[TABLE]

*Therefore, under the mesh constraint (41), we conclude that $U^{n+1}_{h}=0$. This completes the proof. *

Next theorem derives an a priori estimates for $u^{n}_{h}$.

Theorem 3.3**.**

Let $(u^{n}_{h},w^{n}_{h})\in V_{h}\times V_{h}$ be the unique solution of (28)–(29) and suppose the mesh constraint (41) is satisfied, there holds

[TABLE]

Proof 3.4**.**

It suffices to prove the estimates for the solution to (39)–(40). Taking $\eta_{h}=-\Delta^{-1}_{h}\hat{u}^{n+1}_{h}$ in (39) and $v_{h}=\hat{u}^{n+1}_{h}$ in (40), we have

[TABLE]

Notice that

[TABLE]

and the right-hand side of (57) can be estimated as follows:

[TABLE]

and

[TABLE]

by integration by parts and (36). Moreover, it follows from (36) that

[TABLE]

Taking the expectation on both sides of (57), summing over $n=0,1,...,\ell-1$ with $1\leq\ell\leq N$, using (58)–(61) and the fact that

[TABLE]

we get

[TABLE]

*Finally, (56) follows from applying the discrete Gronwall inequality to (62). The proof is complete. *

4 Strong convergence analysis

The goal of this section is to establish the strong convergence with rates for the fully discrete mixed finite element method defined in the previous section. To the end, we introduce for $n=0,1,2,...,N$,

[TABLE]

With the help of Hölder continuity estimates derived in Section 2, we are able to prove strong convergence with rates for $E^{n}$, which is stated in the following theorem.

Theorem 4.1**.**

Under the mesh constraint (41), there holds

[TABLE]

Proof 4.2**.**

Subtracting (28)–(29) from (17)–(18) after substituting [math] by $t_{n}$, and $t$ by $t_{n+1}$, we get $\mathbb{P}$-almost surely

[TABLE]

Since $\Phi^{n+1}(\omega)\in\mathring{V}_{h}$, setting $\eta_{h}=-\Delta^{-1}_{h}\Phi^{n+1}(\omega)$ in (64) and $v_{h}=\tau\Phi^{n+1}(\omega)$ in (65), it follows from the definition of $\Delta^{-1}_{h}$ (cf. (33)) that

[TABLE]

Combining (66) and (67), and taking expectation on both sides, we have

[TABLE]

The left-hand side of (68) can be rewritten as

[TABLE]

Now we estimate the right-hand side of (68). Since $P_{h}$ is the $L^{2}$-projection operator, we have

[TABLE]

For the second term on the right-hand side of (68), we have by (30) that

[TABLE]

For the third term on the right-hand side of (68), we observe that

[TABLE]

First of all, we have

[TABLE]

by the monotonicity property of the nonlinearity. Secondly, we can estimate the first term on the right-hand side of (72) by

[TABLE]

where we had applied (31) and high moment bounds (16) for the SPDE solution. Combining (72)–(74), (30) and (36), we obtain

[TABLE]

For the fourth term on the right-hand side of (68), we have

[TABLE]

by the Hölder continuity for $\nabla w$ (cf. Lemma 2.3). Similarly, we have

[TABLE]

by the Hölder continuity for $\nabla u$ (cf. Lemma 2.1).

For the fifth term on the right-hand side of (68), we have

[TABLE]

For the seventh term on the right-hand side of (68), we have by the integration by parts, the martingale property, the Itô isometry and (36) that

[TABLE]

Similarly, we have

[TABLE]

where we have used Lemma 2.1 and the following Poincaré’s inequality:

[TABLE]

Combining (68)–(71) and (75)–(80), summing over $n=0,1,...,\ell-1$ with $1\leq\ell\leq N$, we have

[TABLE]

Therefore, under the mesh constraint (41), we have by the discrete Gronwall inequality that

[TABLE]

*Finally, the estimate (63) follows from (82), the triangle inequality, and the fact that $\Phi^{0}=0$. The proof is complete. *

Remark 4.3**.**

*The error estimates in Theorem 4.1 is sub-optimal with respect to the $\|\cdot\|_{-1,h}$-seminorm, this is due to the existence of the gradient-type noise, hence, the estimate is sharp in general. Numerical results in Section 5 indeed confirm the sub-optimal convergence whenever the noise is relatively large. However, the error is optimal in the $H^{1}$-seminorm which is also confirmed by the numerical experiments. *

Remark 4.4**.**

*Because the discrete Gronwall inequality was employed near the end of the proof, the error estimate in Theorem 4.1 depends on $\frac{1}{\epsilon}$ exponentially. We note that a polynomial order dependence on $\frac{1}{\epsilon}$ of the errors was achieved in the deterministic case (cf. [14, 15, 13]) by using a PDE spectrum estimate result, however, such a spectrum estimate is yet proved to hold in the stochastic setting. *

5 Numerical experiments

In this section, we report several numerical examples to check the performance of the proposed fully discrete mixed finite element method and numerically study the impact of noise on the evolution of the solution and the stochastic Hele-Shaw flow.

We consider the SPDE (6)–(9) on the square domain $\mathcal{D}=[-1,1]^{2}$ and choose

[TABLE]

where $r=|x|$. It is clear that $\text{div}X=0$ in $\mathcal{D}$ and $X\cdot n=0$ on $\partial\mathcal{D}$.

Let $N_{h}=\text{dim}V_{h}$ and $\{\psi_{i}\}_{i=1}^{N_{h}}$ be the nodal basis of $V_{h}$. Denote by ${\bf u}^{n+1}$ (resp., ${\bf w}^{n+1}$) the coefficient vector of the discrete solution $u^{n+1}_{h}=\sum_{i=0}^{N_{h}}u^{n+1}_{i}\psi_{i}$ (resp., $w^{n+1}_{h}=\sum_{i=0}^{N_{h}}w^{n+1}_{i}\psi_{i}$) at time $t_{n+1}=(n+1)\tau$, $n=0,...,N-1$. Then (28)–(29) are equivalent to

[TABLE]

where $\bf{M}$ and $\bf{A}$ denote respectively the mass and stiffness matrices, $\bf{A_{X}}$ is the weighted stiffness matrix with $({\bf A_{X}})_{ij}=(\nabla\psi_{j}\cdot X,\nabla\psi_{i}\cdot X)$, ${\bf{N}}({{\bf u}^{n+1}})$ is the nonlinear contribution corresponding to the nonlinear term $(f^{n+1},v_{h})$, $({\bf{C_{X}}})_{ij}=(\nabla\psi_{j}\cdot X,\psi_{i})$ and $\bf{W}$ is the discrete Brownian motion with increments $\bar{\Delta}{{\bf W}_{n+1}}={{\bf W}_{n+1}}-{{\bf W}_{n}}$.

In all our tests, we use the Brownian motion generated by using step size $\tau_{\text{ref}}=5\times 10^{-5}$ and compute at least $M=1000$ Monte Carlo realizations. The first test concerns a smooth initial function, aiming to verify the rates of convergence of the proposed method with respect to the temporal mesh size $\tau$ and the spatial mesh size $h$. The second and third tests are designed to investigate the influence of the noise intensity $\delta$ and the parameter $\epsilon$ on the stochastic evolutions for two different non-smooth initial functions.

5.1 Test 1

In this test we check the rates of convergence of the method (28)-(29) with a smooth initial function

[TABLE]

We examine the errors $\sup_{0\leq n\leq N}\mathbb{E}\left[{\|E^{n}\|^{2}_{L^{2}(\mathcal{D})}}\right]$ and $\mathbb{E}\left[{\sum_{n=1}^{N}\tau\|\nabla E^{n}\|^{2}_{L^{2}(\mathcal{D})}}\right]$, where $E^{n}=u(t_{n})-u^{n}_{h}$. Since the exact solution is unknown, we approximate the errors by

[TABLE]

Here $u^{n}_{\text{ref}}$ refers to a reference solution.

First we examine the convergence rates in the time discretization by varying $\tau$ with the fixed parameters $\epsilon=0.1$, $\delta=5$ and $h=2/2^{6}$. In Table 1, we use $L^{\infty}(L^{2})$ (resp., $L^{2}(H^{1})$) to denote the norm (resp., seminorm) corresponding to the square root of the errors $\sup_{0\leq n\leq N}\mathbb{E}\left[{\|E^{n}\|^{2}_{L^{2}(\mathcal{D})}}\right]$ (resp., $\mathbb{E}\left[{\sum_{n=1}^{N}\tau\|\nabla E^{n}\|^{2}_{L^{2}(\mathcal{D})}}\right]$). We observe the half order convergence rate for the $L^{2}(H^{1})$-error as predicted in Theorem 4.1. Note that the $L^{\infty}(L^{2})$ error estimate is not available in Theorem 4.1, however it still provides us useful information about the accuracy of the numerical method.

Next we investigate convergence in the space discretization by varying $h$ with the fixed parameters $\epsilon=0.1$, $\delta=25$ and $\tau=5\times 10^{-5}$. Since $\delta$ is relatively large, we compute $M=10^{4}$ realizations. From Table 2, we see that the $L^{2}(H^{1})$-error converges with order $1$ which is consistent with the theoretical estimate of Theorem 4.1. Also, the $L^{\infty}(L^{2})$-error converges with order less than $2$, indicating a sub-optimal convergence in the lower order norm as predicted by Theorem 4.1.

5.2 Test 2

In this test, we take the initial function to be

[TABLE]

where $d_{0}(x)$ represents the signed distance function to the ellipse

[TABLE]

First, we investigate the evolution of the zero-level set with respect to the noise intensity $\delta$ with fixed $\epsilon=0.01$. Figure 1 plots snapshots at several time points of the zero-level set of $\bar{u}_{h}$ for $\delta=1,5,10$, where

[TABLE]

We observe that when the noise is relatively small ($\delta=1$), the zero-level set is close to the deterministic interface (for $\delta=0$). However, for relatively large noises ($\delta=5,10$), the zero-level sets rotate and evolve faster.

Next, we fix $\delta=1$ and study the influence of the parameter $\epsilon$ on the evolution of the numerical interfaces. In Figure 2, snapshots at four fixed time points of the zero-level set of $\bar{u}_{h}$ are depicted for three different $\epsilon=0.01,0.015,0.04$. Numerical results suggest the convergence of the numerical interface to the stochastic Hele-Shaw flow as $\epsilon\to 0$ at each of four time points. In addition, the numerical interface evolves faster in time for larger $\epsilon$.

Notice that in Figure 1–2, we only plot the evolutions on the subdomain $[-0.6,0.6]^{2}$ for a better resolution.

In Figure 3, we plot the change of the expected value of the discrete energy

[TABLE]

in time with fixed $\epsilon=0.01$. It indicates that the decay property still holds for $\delta=1,5,10$.

5.3 Test 3

In this test, we consider the case with

[TABLE]

where $d_{0}(x)=\min\{d_{1}(x),d_{2}(x)\}$, and $d_{1}(x)$ and $d_{2}(x)$ denotes respectively the signed distance function to the ellipses

[TABLE]

In Figure 4, we depict snapshots at several time points of the zero-level set of $\bar{u}_{h}$ for $\delta=1,5,10$ with fixed parameter $\epsilon=0.01$. For all cases, the two separated zero-level sets eventually merge and evolve to a circular shape. For larger noise intensity ($\delta=5,10$), the two interfaces merge faster and develop two concentric interfaces where the outer interface evolves to a circular shape and the inner interface shrinks and eventually vanishes.

Next, we plot a few snapshots of the zero-level set of $\bar{u}_{h}$ for $\epsilon=0.01,0.015,0.04$ with fixed $\delta=1$ in Figure 5. Again, the numerical interface evolves faster in time for larger $\epsilon$, and the numerical interfaces stay close for $\epsilon=0.01$ and $\epsilon=0.015$.

The decay of the expected value of the discrete energy is shown in Figure 6, where we consider three noise intensity levels $\delta=1,5,10$ with fixed $\epsilon=0.01$.