Analysis of a Legendre spectral element method (LSEM) for the two-dimensional system of a nonlinear stochastic advection-reaction-diffusion models

Analysis of a Legendre spectral element method (LSEM) for the two-dimensional system of a nonlinear stochastic advection-reaction-diffusion models

Department of Applied Mathematics, Faculty of Mathematics and Computer Sciences,
Amirkabir University of Technology, No. 424, Hafez Ave., 15914, Tehran, Iran
Institute of Applied Mathematics, Leibniz University of Hannover

Welfengarten 1, 30167 Hanover, Germany
April 15, 2019
Abstract

In this work, we develop a Legendre spectral element method (LSEM) for solving the stochastic nonlinear system of advection-reaction-diffusion models. The used basis functions are based on a class of Legendre functions such that their mass and diffuse matrices are tridiagonal and diagonal, respectively. The temporal variable is discretized by a Crank–Nicolson finite difference formulation. In the stochastic direction, we also employ a random variable based on the Wiener process. We inspect the rate of convergence and the unconditional stability for the achieved semi-discrete formulation. Then, the Legendre spectral element technique is used to obtain a full-discrete scheme. The error estimation of the proposed numerical scheme is substantiated based upon the energy method. The numerical results confirm the theoretical analysis.

Keywords: Nonlinear system of advection-reaction-diffusion equation, error estimate, spectral element method (SEM), stochastic PDEs, .

AMS subject Classification: 65M70, 34A34.

1 Introduction

We consider the stochastic nonlinear system of advection-reaction-diffusion models [1, 2]

 (1.1)

where , and denote the concentrations of the main ground substance, aqueous solution electrolyte and microorganism, respectively [1, 2]. In the above model is a known function, is the advection coefficient, is the diffusion coefficient, and are constant, respectively. Also, is a -Wiener process with respect to a filtered probability space . The nonlinear terms are

 f(u,v)=g(u,v)=h(u,v)=uκ1+u+vκ2+v.

Predictions of solute transport in aquifers generally have to rely on mathematical models based on groundwater flow and convection-dispersion equations. The groundwater model is employed to prevent and control the groundwater contaminant with the microbiological technology [2]. Several scholars investigated Eq. (1.1) for example using an improved finite element approach [2], meshless local approaches [3, 4], lattice Boltzmann technique [5], a front-tracking method [6], novel WENO methods [7], or a finite element method [8]. The interested readers can refer to [9, 10] to get more information for Eq. (1.1).

In the past, the groundwater models have been based only on deterministic considerations. In practice, aquifers are generally heterogeneous, i.e., their hydraulic properties (e.g., permeability) change in space. These variations are irregular and characterized by length scales significantly larger than the pore scale. These spatial fluctuations cause the flow variables such as concentration to change in space in an irregular manner. Therefore, a reliable description of the groundwater model can be explained only in a stochastic form [11].

The first stochastic equation can be rewritten as

 du(t)=(Au(t)+f(u))dt+dW, (1.2)

where is a linear, self-adjoint, positive definite operator where the domain is dense in and compactly embedded in (i.e., ) and the semigroup is generated by . Additionally, we assume that satisfies the linear growth condition and is twice continuously Frechet differentiable with bounded derivatives up to order 2 [12]. The initial value is deterministic as well. Therefore, (1.2) has a continuous mild solution [13]

 u(t)=etAu0+∫t0e(t−s)Af(u(s)) ds+∫t0e(t−s)A dW(s), (1.3)

where for and . Regarding the expected value of the solution, we can assume that . The same mild solutions can be employed for and .

The deterministic case of Eq. (1.1) has been studied by some scholars for example a new finite volume method [2], new Krylov WENO methods [7], local radial basis function collocation method [14], etc. Also, the SEM is applied to solve some important problems such as the Schrödinger equations [15], Pennes bioheat transfer model [16], the shallow water equations [17], integral differential equations [18, 19, 20], hyperbolic scalar equations [21], predator-prey problem [22], some problems in the finance mathematics [23, 24] and so forth.

The main aim of the current paper is to propose a new high-order numerical procedure for solving the two-dimensional system of a nonlinear stochastic advection-reaction-diffusion models. The used technique is based on the modified Legendre spectral element procedure. The coefficient matrix of the employed technique is more well-posed than the traditional Legendre spectral element method. The structure of this article is as follows. In Section 2, we propose and analysis the time-discrete scheme. In Section 3, we develop the new numerical technique and analysis it. We check the numerical results to solve the considered model in Section 4. Finally, a brief conclusion of the current paper is written in Section 5.

2 Temporal discretization

First of all, we briefly review some important notations used in the paper. Considering , we define the following functional spaces

 L2(Ω)={f:∫Ωf2dΩ<∞},H1(Ω)={f∈L2(Ω),∇f∈L2(Ω)},H10(Ω)={f∈H1(Ω),f|∂Ω=0},Hk(Ω)={f∈L2(Ω),Dβf∈L2(Ω)forall|β|≤k},

and the derivative

 Dαf=(∂α1f∂xα11)(∂α2f∂xα22)…(∂αpf∂xαpp),|α|=p∑i=1αi.

The corresponding inner products for and are as follows

 (f,g)=∫Ωf(x)g(x)dΩ,(f,g)1=(f,g)+(∇f,∇g),

and the associated norms are

Furthermore, associated norm for the space is as

 ∥f∥Hm(Ω)=⎛⎝∑0≤|α|≤m∥Dαf∥2L2(Ω)⎞⎠12.

To discretize the time variable, we define

 tn=nτ,∀ n=0,1,…,N,

where is the step size. We introduce additionally

 vn−12=v(x,y,tn−12)=12(vn+vn−1),δtvn−12=1τ(vn−vn−1),vn=v(x,y,tn).

The Crank-Nicolson scheme for problem (1.1) is as follows

 (2.1)

where is a positive constant such that Discretizing relation (2.1) yields

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩un−un−1τ+ξ(\mathbfitx)[∇un+∇un−12]−∇⋅[ζ(\mathbfitx)(∇un+∇un−12)]+wpe1f(un−12,vn−12)=Wn−Wn−1τ,vn−vn−1τ+ξ(\mathbfitx)[∇vn+∇vn−12]−∇⋅[ζ(\mathbfitx)(∇vn+∇vn−12)]+wpe2f(un−12,vn−12)=Wn−Wn−1τ,wn−wn−1τ+ξ(\mathbfitx)[∇wn+∇wn−12]−∇⋅[ζ(\mathbfitx)(∇wn+∇wn−12)]+wpe3f(un−12,vn−12)+r(\mathbfitx)[wn+wn−12]=Wn−Wn−1τ, (2.2)

or

 (2.3)

The vector-matrix configuration of Eq. (2.3) is

 \mathbfitH1Un+τ2\mathbfitI∇Un−τ2\mathbfitI∇⋅ζ(\mathbfitx)∇Un+Wn=\mathbfitH2Un−1−τ2\mathbfitI∇Un−1+τ2\mathbfitI∇⋅ζ(\mathbfitx)∇Un−1−τ\mathbfitNF(Un−1)+Wn−1, (2.4)

where is the identity matrix and

 \mathbfitH1=diag(1,1,1+τ2r(\mathbfitx)),\mathbfitH2=diag(1,1,1−τ2r(\mathbfitx)),\mathbfitN=diag(wpe1,wpe2,wpe3), (2.5)

and also the unknown vector is .

2.1 Error analysis of the semi-discrete formulation

Theorem 2.1.

If , then relation (2.4) will be unconditionally stable.

Proof.

Let and . We want to find such that

 \mathbfitH1(Un,χ) + τ2\mathbfitI(ζ(x)∇Un,∇χ)−τ2\mathbfitI(Un,∂∂xχ)−τ2\mathbfitI(Un,∂∂yχ)+(Wn,V) (2.6) = \mathbfitH2(Un−1,χ)−τ2\mathbfitI(ζ(x)∇Un−1,∇χ)+τ2\mathbfitI(Un−1,∂∂xχ) + τ2\mathbfitI(Un−1,∂∂yχ)−τ\mathbfitN(\mathbfitF,χ)+(Wn−1,V)∀ χ∈\mathbfitH10(Ω).

Let be an approximate solution of , then

 \mathbfitH1(˜Un,χ) + τ2\mathbfitI(ζ(x)∇˜Un,∇χ)−τ2\mathbfitI(˜Un,∂∂xχ)−τ2\mathbfitI(Un,∂∂yχ)+(Wn,V) (2.7) = \mathbfitH2(˜Un−1,χ)−τ2\mathbfitI(ζ(x)∇˜Un−1,∇χ)+τ2\mathbfitI(˜Un−1,∂∂xχ) + τ2\mathbfitI(˜Un−1,∂∂yχ)−τ\mathbfitN(\mathbfit˜F,χ)+(Wn−1,V)∀ χ∈\mathbfitH10(Ω),

where . Subtracting Eq. (2.7) for Eq. (2.6) , results

 \mathbfitH1(\mathbfitΨn,χ) + τ2\mathbfitI(ζ(x)∇\mathbfitΨn,∇χ)−τ2\mathbfitI(\mathbfitΨn,∂∂xχ)−τ2\mathbfitI(\mathbfitΨn,∂∂yχ) (2.8) = \mathbfitH2(\mathbfitΨn−1,χ)−τ2\mathbfitI(ζ(x)∇\mathbfitΨn−1,∇χ)+τ2\mathbfitI(\mathbfitΨn−1,∂∂xχ) + τ2\mathbfitI(\mathbfitΨn−1,∂∂yχ)−τ\mathbfitN(\mathbfitF−˜F,χ),∀ χ∈\mathbfitH10(Ω),

where

 \mathbfitΨn=E[Un−˜Un].

Setting in Eq. (2.8) yields

 \mathbfitH1(\mathbfitΨn,\mathbfitΨn) + τ2\mathbfitI(ζ(x)∇\mathbfitΨn,∇\mathbfitΨn)−τ2\mathbfitI(\mathbfitΨn,∂∂x\mathbfitΨn)−τ2\mathbfitI(\mathbfitΨn,∂∂y\mathbfitΨn) (2.9) = \mathbfitH2(\mathbfitΨn−1,\mathbfitΨn)−τ2\mathbfitI(ζ(x)∇\mathbfitΨn−1,∇\mathbfitΨn)+τ2\mathbfitI(\mathbfitΨn−1,∂∂x\mathbfitΨn) + τ2\mathbfitI(\mathbfitΨn−1,∂∂yχ)−τ\mathbfitN(\mathbfitF−˜F,\mathbfitΨn).

Applying the Cauchy-Schwarz inequality for Eq. (2.9), results

 ∥∥\mathbfitH1∥∥∥∥\mathbfitΨn∥∥2L2(Ω) + τ2∥ζ(x)∥∥∥∇\mathbfitΨn∥∥2L2(Ω)≤τ2(\mathbfitΨn,∂∂x\mathbfitΨn)+τ2(\mathbfitΨn,∂∂y\mathbfitΨn) + + τ2(\mathbfitΨn−1,∂∂y\mathbfitΨn)+τ2(\mathbfitΨn−1,∂∂x\mathbfitΨn)−τ\mathbfitN(\mathbfitF−˜F,\mathbfitΨn).

There exists constant such that

 ∥∥\mathbfitH2∥∥,∥∥\mathbfitH3∥∥≤C, (2.10)

and

 ∥∥\mathbfitF−˜\mathbfitF∥∥≤L\mathbfitΨn−1. (2.11)

By simplification we have

 ∥∥\mathbfitH1∥∥∥∥\mathbfitΨn∥∥2L2(Ω) + τ2∥ζ(x)∥∥∥∇\mathbfitΨn∥∥2L2(Ω)≤τ2∥∥\mathbfitΨn∥∥L2(Ω)∥∥∇\mathbfitΨn−1∥∥L2(Ω)+τ2∥∥\mathbfitΨn−1∥∥L2(Ω)∥∥∇\mathbfitΨn∥∥L2(Ω) + + τL∥∥\mathbfitN∥∥∥∥\mathbfitΨn−1∥∥L2(Ω)∥∥\mathbfitΨn∥∥L2(Ω).

So, from the following assumption and the definition of matrices and , we have

 ∥∥\mathbfitH2∥∥≤∥∥\mathbfitH1∥∥.

Now, we can get

 12∥∥\mathbfitH1∥∥∥∥\mathbfitΨn∥∥2L2(Ω) + τ4∥∥ζ(\mathbfitx)∥∥∥∥∇\mathbfitΨn∥∥2L2(Ω) ≤ 12∥∥\mathbfitH1∥∥∥∥\mathbfitΨn−1∥∥2L2(Ω)+τ2∥∥ζ(\mathbfitx)∥∥∥∥∇\mathbfitΨn−1∥∥2L2(Ω) + C1Lτ2∥∥ζ(\mathbfitx)∥∥∥∥\mathbfitΨn∥∥2L2(Ω)+C2Lτ2∥∥ζ(\mathbfitx)∥∥∥∥\mathbfitΨn−1∥∥2L2(Ω).

Using the below relation

 ∥∥\mathbfitΨn∥∥2\mathbfitHw(Ω)=∥∥\mathbfitH1∥∥∥∥\mathbfitΨn∥∥2L2(Ω)+12τ∥∥ζ(\mathbfitx)∥∥∥∥∇\mathbfitΨn∥∥2L2(Ω),

Eq. (2.1) is changed to

 ∥∥\mathbfitΨn∥∥2\mathbfitHw(Ω)≤∥∥\mathbfitΨn−1∥∥2\mathbfitHw(Ω)+C1Lτ∥∥ζ(\mathbfitx)∥∥∥∥\mathbfitΨn∥∥2\mathbfitHw(Ω)+C2Lτ∥∥ζ(\mathbfitx)∥∥∥∥\mathbfitΨn−1∥∥2\mathbfitHw(Ω). (2.13)

By summing Eq. (2.13) for from 0 to , gives

Thus, we have

 ∥∥\mathbfitΨn∥∥2\mathbfitHw(Ω)≤∥∥\mathbfitΨ0∥∥2\mathbfitHw(Ω)+2CLτ∥∥ζ(\mathbfitx)∥∥n∑m=1∥∥\mathbfitΨm∥∥2\mathbfitHw(Ω) (2.14)

Considering Gronwall’s inequality for Eq. (2.14) yields

 ∥∥\mathbfitΨn∥∥2\mathbfitHw(Ω) ≤ ∥∥\mathbfitΨ0∥∥2\mathbfitHw(Ω)+2CLτ∥∥ζ(\mathbfitx)∥∥n∑m=1∥∥\mathbfitΨm∥∥2\mathbfitHw(Ω) ≤ {∥∥\mathbfitΨ0∥∥2\mathbfitHw(Ω)}exp(2CLnτ∥∥ζ(\mathbfitx)∥∥)≤C∥∥\mathbfitΨ0∥∥2\mathbfitHw(Ω).

So, we have

 ∥∥\mathbfitΨn∥∥L2(Ω)≤∥∥\mathbfitΨn∥∥\mathbfitHw(Ω)≤C∥∥\mathbfitΨ0∥∥\mathbfitHw(Ω).

Theorem 2.2.

The convergence order of relation (2.4) is .

Proof.

Let us assume . We set

 \mathbfitXn=E[un−Un]n≥1,

where . Then, we have

 \mathbfitH1\mathbfitXn+τ2\mathbfitI∇\mathbfitXn−τ2\mathbfitI∇⋅ζ(\mathbfitx)∇\mathbfitXn=\mathbfitH2\mathbfitXn−1−τ2\mathbfitI∇\mathbfitXn−1+τ2\mathbfitI∇⋅ζ(\mathbfitx)∇\mathbfitXn−1+τ\mathbfitR−τ\mathbfitN(\mathbfitFn−1−˜\mathbfitFn−1). (2.15)

According to the Crank-Nicolson idea, we have

 |R|≤C1τ2.

Similar to Theorem 2.1, we obtain

 ∥∥\mathbfitXn∥∥2\mathbfitHw(Ω) ≤ ∥∥\mathbfitX0∥∥2\mathbfitHw(Ω)+2Lτ∥∥ζ(\mathbfitx)∥∥n∑m=1∥∥\mathbfitXm∥∥2\mathbfitHw(Ω)+max1≤m≤n∥∥\mathbfitR∥∥2L2(Ω) ≤ {max1≤m≤n∥∥\mathbfitR∥∥2L2(Ω)}exp(2Lnτ∥∥ζ(\mathbfitx)∥∥)≤Cτ2 ≤ {max1≤m≤n∥R∥2L2(Ω)}exp(2Lnτ∥ζ(x)∥) ≤ exp(2Lnτ∥ζ(x)∥)C1τ2≤Cτ2.

which completes the proof. ∎

3 Error estimation for full-discrete plane

In this section, we employ a new class of Legendre polynomial functions which were developed in [25].

Lemma 3.1.

[25] Consider the following relations

 ψk(x)=γk(Lk(x)−Lk+2(x)), (3.1)

in which and are the Legendre polynomials. Let us denote

 ajk=1∫−1dψk(x)dxdψj(x)dxdx,bjk=1∫−1ψk(x)ψj(x)dx. (3.2)

Then

 ajk=⎧⎪⎨⎪⎩1,k=j,0,k≠j,bjk=bkj=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩γkγj(22j+1+22j+5),k=j,−γkγj22k+1,k=j+2,0,Otherwise. (3.3)

The SEM as a combination of the finite element method and spectral polynomials has been developed by Patera [26]. By dividing the computational region into non-overlapping elements

 Ω=Ne⋃e=1Ωe,Ωi∩Ωj=∅,i≠j.

Now, we define the following projection operator.

 P1h:H10(Ω)→V0N, (3.4)

where

 (∇(u−P1hu),∇v)=0,u∈H10(Ω),∀v∈V0N, (3.5)

and is the spectral element approximation space

 V0N={w∈H10(Ω):w|Ωs∈PN(Ω),s=1,2,…,ns}. (3.6)
Lemma 3.2.

[27] Let (), therefore

 ∥∥u−P1hu∥∥≤C[ns∑k=1h2(min(Nk+1,υ)−1)kN2(1−υ)k∥u∥2υ]12. (3.7)

In the special cases and we get

 ∥∥u−P1hu∥∥≤Ch(min(N+1,υ)−1)kN1−υ∥u∥υ. (3.8)

We aim to find a such that

 \footnotesize\mathbfitH1(Un,χ) + τ2\mathbfitI(ζ(x)∇Un,∇χ)−τ2\mathbfitI(Un,∂∂xχ)−τ2\mathbfitI(Un,∂∂yχ)+(Wn,V) (3.9) = \mathbfitH2(Un−1,χ)−τ2\mathbfitI(ζ(x)∇Un−1,∇χ)+τ2\mathbfitI(Un−1,∂∂xχ) + τ2\mathbfitI(Un−1,∂∂yχ)−τ\mathbfitN(\mathbfitFn−1,χ)+τ(\mathbfitRnt,χ)+(Wn−1,V)χ∈\mathbfitH10(Ω).

The spectral element formulation is: find a such that

 \mathbfitH1(Unh,χh) + (3.10) = \mathbfitH2(Un−1h,χh)−τ2\mathbfitI(ζ(x)∇Un−1h,∇χh)+τ2\mathbfitI(Un−1h,∂∂xχh) + τ2\mathbfitI(Un−1h,∂∂yχ)−τ\mathbfitN(\mathbfitFn−1,χh)+(Wn−1,V)∀ χh∈\mathbfitH10(Ω).
Lemma 3.3.

Let

 (Gnr,d,χr) = (P1hˆUnr−ˆUnr,χr)+τ\mathbfitA2(P1hˆUnr−ˆUnr,∂∂xχr)+τ\mathbfitA3(P1hˆUnr−ˆUnr,∂∂yχr) − (P1hˆUn−1r−ˆUn−1r,χr)+τ\mathbfitA2(P1hˆUn−1r−ˆUn−1r,∂∂xχr