Error analysis of nonlinear time fractional mobile/immobile advection-diffusion equation with weakly singular solutions

Hui Zhang; Xiaoyun Jiang; Fawang Liu

doi:10.1515/fca-2021-0009

Artikel Öffentlich zugänglich

Error analysis of nonlinear time fractional mobile/immobile advection-diffusion equation with weakly singular solutions

Hui Zhang , Xiaoyun Jiang und Fawang Liu

Veröffentlicht/Copyright: 29. Januar 2021

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Fractional Calculus and Applied Analysis Band 24 Heft 1

Abstract

In this paper, a weighted and shifted Grünwald-Letnikov difference (WSGD) Legendre spectral method is proposed to solve the two-dimensional nonlinear time fractional mobile/immobile advection-dispersion equation. We introduce the correction method to deal with the singularity in time, and the stability and convergence analysis are proven. In the numerical implementation, a fast method is applied based on a globally uniform approximation of the trapezoidal rule for the integral on the real line to decrease the memory requirement and computational cost. The memory requirement and computational cost are O(Q) and O(QK), respectively, where K is the number of the final time step and Q is the number of quadrature points used in the trapezoidal rule. Some numerical experiments are given to confirm our theoretical analysis and the effectiveness of the presented methods.

MSC 2010: Primary 26A33; Secondary 65M06; 65M12; 65M15; 65M70; 35R11

Key Words and Phrases: two-dimensional nonlinear time fractional mobile/immobile advection-dispersion equation; correction method; stability and convergence analysis; fast method

1 Introduction

In this paper, we consider the following two-dimensional nonlinear time fractional mobile/immobile advection-diffusion equation [33, 42],

ut+β0CDtαu=−V∇u+DΔu+F(u)+f(x,y,t),x,y∈Ω,t∈I,u(x,y,0)=u0(x,y),x,y∈Ω,u(x,y,t)=0,x,y∈∂Ω,t∈I,(1.1)

where Ω = (a, b) × (c, d), I = (0, T], u denotes the solute concentration in the total (mobile+immobile) phase, and β is the fractional capacity coefficient. Here V and D are the velocity and dispersion coefficient for the mobile phase (and hence may be directly measured). The time drift term u_t is added to describe the motion time and thus helps to distinguish the status of particles conveniently (see also the discussion in [2]). F(u) is a nonlinear term satisfies the local Lipschitz condition. 0CDtα (0 < α < 1) is the Caputo fractional derivative of order α that can be defined [31] as follows,

0CDtαu(t)=1Γ(1−α)∫0tu′(s)ds(t−s)α.(1.2)

Many researchers have studied the nonlinear time fractional mobile/immobile advection-diffusion equation (1.1). A RBF meshless approach was introduced for modeling a fractal mobile/immobile transport model [22]. Zhang et al. [40] proposed a novel numerical method for the time variable fractional order mobile/immobile advection-dispersion model. For a more detailed description of this aspect, we refer the readers to [27, 42]. A variety of numerical methods have been used to approximating the fractional partial differential equations [1, 5, 6, 7, 11, 12, 13, 14, 15, 18, 20, 21, 30]. Among them, spectral method has attracted more and more attentions [3, 9, 17, 36, 37] because of its high-order accuracy. In this paper, we develop a WSGD Legendre spectral method for the considered model.

In practical applications, analytical solutions of time fractional differential equations may have strong singularity at t = 0 [28, 35]. In order to solve the singularity, some researchers [16, 32] adopted the graded mesh method. For the considered equation, we introduce some correction terms [38, 39] to deal with the non-smooth solutions. The stability and convergence analysis of the correction method are also proven.

We can rewrite the model as

ut+β0RLDtα(u−u0)=−V∇u+DΔu+F(u)+f(x,y,t),(1.3)

here 0RLDtα (0 < α < 1) is the Riemann-Liouville fractional derivative of order α [31]. Using WSGD method yields the following discrete convolution [25] as

0RLDtαu(tn)=1τα∑l=0nλn−l(α)ul,0≤n≤K,(1.4)

where τ is the time step size, λn−l(α) are the convolution quadrature weights. The implementation of (1.4) requires O(K) active memory and O(K²) arithmetic operations by direct computation. So the direct calculation of (1.4) becomes computational expensive when it is used to discretize time fractional partial differential equations. Some works [19, 24, 39] have been presented to reduce the memory requirements and computational cost of the discrete convolution for approximating the Riemann-Liouville fractional derivative. In this numerical implementation, we apply a fast method, in which the Hankel contour is used to express the quadrature weight λn(α) as an integral on the half line, see later for details.

The highlights of this paper are: (i) we introduce some correction terms to solve the two-dimensional nonlinear time fractional mobile/immobile advection-diffusion equation with non-smooth solutions; (ii) stability and convergence analysis of correction method are proven; (iii) a fast method is applied to decrease the memory requirement and computational cost.

The paper is organized as follows. Section 2 develops the WSGD Legendre spectral method to solve the two-dimensional nonlinear time fractional mobile/immobile advection-diffusion equation. And some correction terms are introduced to deal with the non-smooth solutions. The stability and convergence analysis of the correction method are presented in Section 3. In Section 4, a fast algorithm for numerical implementation is proposed. Numerical examples are provided in Section 5. Finally, some conclusions and discussions are summarized in Section 6.

2 Numerical method

For the temporal discretization, we divide the interval [0, T] into K equal subintervals with a time-step size τ = T/K. Let t_n = nτ, uⁿ = u(x, y, t_n), 0 ≤ n ≤ K. For convenience, we introduce the following notations:

δtun+1/2=un+1−un+1τ,un+1/2=un+un+12.

We approximate the fractional derivative as (1.3),

λ0(α)=1+α2g0(α),λj(α)=1+α2gj(α)−α2gj−1(α),j≥1,(2.1)

and gj(α)=(−1)kαj for j ≥ 0, and they can be evaluated recursively by

g0(α)=1,gj(α)=1−α+1jgj−1(α),j=1,2,….(2.2)

If α = 0, we set g0(α)=1,gj(α)=0 for j ≥ 1. It is obvious that

δtun+1/2+βDτα,n+1/2u~=−V∇un+1/2+DΔun+1/2+F(un+1/2)+fn+1/2+En,(2.3)

where u~n=un−u0,Dτα,nu=1τα∑l=0nλn−l(α)ul, and Eⁿ is the truncation error satisfying |Eⁿ| ≤ O(τ²).

The WSGD scheme can not preserve second-order accuracy in time approximation when u has a strong singularity. The key assumption is that the analytical solution u(x, y, t) satisfies the following form [4, 26],

u−u0=∑j=1mcj(x,y)tσj+cm+1(x,y)tσm+1+⋯,(2.4)

where 1 ≤ σ_j < σ_j+1. We develop the correction method to deal with the non-smooth solutions [38],

0RLDtαun=Dτα,nu+1τα∑j=1m1Wn,jα(uj−u0)+R1n=Dτα,n,m1u+R1n,(2.5)

where m₁ ≤ m, Wn,jα are the starting weights satisfying

Dτα,nu+1τα∑j=1m1Wn,jα(uj−u0)=Γ(1+σj)Γ(1+σj−α)tnσj−α,(2.6)

in which u = t^σ_j. We can get Wn,jα (1 ≤ j ≤ m) from (2.6) and Wn,jα is independent of τ for n > 0. For the first-order time derivative, it has

∂tun+1/2=δtun+1/2+1τ∑j=1m2Wn,j(uj−u0)+R2n=δtm2un+1/2+R2n,(2.7)

where m₂ ≤ m, W_n,j can be obtained by the following equation

σj2(tn+1σj−1+tnσj−1)=δtun+1/2+1τ∑j=1m2Wn,j(uj−u0),(2.8)

for u = t^σ_j.

Based on boundary conditions, we use the following basis functions [37] in the x and y directions for the spatial discretization,

ϕi(x)=Li(x^)−Li+2(x^),φi(y)=Li(y^)−Li+2(y^),(2.9)

where i = 0, 1, …, N − 2 and L_i(ẑ) is the Legendre polynomial [34], and x^∈(−1,1),x=(b−a)x^+a+b2∈(a,b),y^∈(−1,1),y=(d−c)y^+c+d2∈(c,d). The function space S_N can be specified as follows,

SN=span{ϕi(x)φj(y),i,j=0,1,…,N−2}.

Combing (2.3) and the above results leads to the following fully discrete form,

(δtm2uNn+1/2,v)+(Dτα,n+1/2,m1u~N,v)+V(∇uNn+1/2,v)+D(∇uNn+1/2,∇v)=F(uNn+1/2),v+(fn+1/2,v),∀v∈SN,(2.10)

and

uN0=ΠN1,0u0,(2.11)

where ΠN1,0 denotes the orthogonal projection operator [37]: H01(Ω) → S_N,

(∇(u−ΠN1,0u),∇v)=0,∀v∈SN.(2.12)

We rewrite the equation as

Gn(v)=1+βτ1−αλ0(α)2(uNn+1,v)+τV2(∇uNn+1,v)+τD2(∇uNn+1,∇v).(2.13)

Next, we derive the matrix representation. The numerical solution can be denoted by

uNn+1=∑i=1N+1∑j=1N+1cijn+1ϕi(x)φj(y).(2.14)

Letting v = ϕ_hφ_l (h, l = 0, 1, …, N − 2), we have

∑i=0N−2∑j=0N−2{1+βτ1−αλ0(α)2(ϕiφj,ϕhφl)+τV2(∂xϕiφj,ϕhφl)+τV2(ϕi∂yφj,ϕhφl)+τD2(∂xϕiφj,∂xϕhφl)+τD2(ϕi∂yφj,ϕh∂yφl)}cijn+1=Gn(ϕhφl).(2.15)

The matrix representation is shown as

(1+βτ1−αλ0(α)2Mx⊗My+τV2(Px⊗My+Mx⊗Py)+τD2(Sx⊗My+Mx⊗Sy))Cn+1=Gn,(2.16)

where

(Mx)h,i=(ϕi,ϕh),(My)l,j=(φj,φl),(Px)h,i=(∂xϕi,ϕh),(Py)l,j=(∂yφi,φj),(Sx)h,i=(∂xϕi,∂xϕh),(Sy)l,j=(∂yφj,∂yφl),(Cn+1)i,j=cijn+1,i,j=0,1,…,N−2.

Remark 2.1

When this model has non-homogeneous boundary conditions, we use the following basis functions [37] in the x and y directions for the spatial discretization,

ϕ1(x)=L0(x^)−L1(x^)2,φ1(y)=L0(y^)−L1(y^)2,ϕi(x)=Li−2(x^)−Li(x^),φi(y)=Li−2(y^)−Li(y^),ϕN+1(x)=L0(x^)+L1(x^)2,φN+1(y)=L0(y^)+L1(y^)2,(2.17)

where i = 2, 3, …, N.

3 Theoretical analysis of the correction method

3.1 Preliminaries

The following lemmas are given used in the later proof.

Lemma 3.1

([23]). As for(1.4), there exists any positive integerLand real vector (u₀, u₁, u₂, …, u_L)^T ∈ ℝ^L+1, it holds that

∑n=0L∑j=0nλj(α)un−jun≥0.

Lemma 3.2

([37]). Letτ, ganda_k, γ_k, for integersk > 0, be nonnegative numbers such that

an≤τ∑k=0n−1γkak+g,n≥0.

Then,

an≤gexpτ∑k=0n−1γk.

Lemma 3.3

([38]). Letm₁be a positive integer andWn,jα (j = 1, 2, …, m₁) be defined by(2.6). Suppose thatσ_jare a sequence of increasing positive numbers. We have

Wn,jα=O(nσ1−2−α)+O(nσ2−2−α)+⋯+O(nσm1−2−α).

Lemma 3.4

([38]). Letu(t) = t^σ (σ ≥ 0) andαbe a real number. There exists a positive and bounded constantCindependent ofnandτsuch that the errorR1ndefined in(2.5)has the following estimate,

|R1n|≤CSm1στσ−αnσmax−α−2,

whereSm1σ=Πk=1m1|σ−σk|,σmax=max{σ,σ1,…,σm1}.

Lemma 3.5

([37]). For any u ∈ S_N, there exists a positive constantCindependent onn, τ, andNsuch that the inverse inequality is obtained,

∥u∥L∞≤CN∥u∥.

Lemma 3.6

([41]). Suppose that u ∈ H0r(Ω), there exists a positive constantCindependent onn, τ, andNsuch that

∥u−ΠN1,0u∥s≤CNs−r∥u∥r,s=0,1,r≥1.

3.2 Stability and convergence analysis

We first present the following stability result.

Suppose that uNk has the error u^Nk. We obtain the error equation,

(δtm2u^Nk+1/2,v)+(Dτα,k+1/2,m1u~^N,v)+V(∇u^Nk+1/2,v)+D(∇u^Nk+1/2,∇v)=F^Nk+1/2,v+f^k+1/2,v,(3.1)

where F^Nk=F(uNk+u^Nk)−F(uNk). Let C be an appropriate positive constant and Θ be a suitable domain. Suppose that F′(y) is bounded in the appropriate domain Θ. Denoting

uNmax=max0⩽k⩽L∥uNk∥∞,FNmax=maxy∈Θ|F′(y)|.(3.2)

We assume

∥u^Nk∥∞⩽C0,0⩽k⩽n⩽m,k,m∈Z+.(3.3)

The following assumption can be given

∥u^Nk+1∥∞⩽C1,if∥u^Nk∥∞⩽C0,(3.4)

where C₁ is a positive constant independent of τ, N and k.

Letting v=u^Nk+1/2 in (2.10) yields,

∥u^Nk+1∥2−∥u^Nk∥22τ+β(Dτα,k+1/2u^N,u^Nk+1/2)+D∥∇u^Nk+1/2∥2=−∑j=1m2Wk,j((u^Nj−u^N0)/τ,u^Nk+1/2)+β(Dτα,k+1/2u^N0,u^Nk+1/2)−τ1−α2∑j=1m1(Wk,jα+Wk+1,jα)((u^Nj−u^N0)/τ,u^Nk+1/2)+V(u^Nk+1/2,∇u^Nk+1/2)+(F^Nk+1/2,u^Nk+1/2)+(f^k+1/2,u^Nk+1/2).(3.5)

Based on the result in [38] and Lemma 3.3, we obtain

qk=∑l=0kλl(α)=O(k−α),|Wk,jα|≤Ckσm1−2−α,|Wk,j|≤Ckσm2−3.(3.6)

Combining Cauchy-Schwarz inequality and Lemma 3.1, we can easily get

Since F(u) satisfies the local Lipschitz condition, it is shown that

∥F^Nk∥=∥F(uNk+u^Nk)−F(uNk)∥≤C∥u^Nk∥.(3.8)

Summing n in (3.7) from 0 to n gives

Using (3.6), it infers that,

According to the assumption, we have

A1n=∑k=1n|Wk,jα+Wk+1,jα|≤C∑k=1nkσm1−2−α≤Cmax{1,nσm1−1−α},A2n=∑k=1n|Wk,j|≤C∑k=0nkσm2−3≤Cmax{1,nσm2−2}.(3.11)

For sufficiently small τ, we can obtain the following inequality easily by (3.6) and (3.11),

∥u^Nn+1∥2≤ρn+C∑k=1nbk∥u^Nk∥2,(3.12)

where

ρn=C(∥u^N0∥2+A2n∑j=1m2∥δtu^Nj−1/2∥2+τ1−αA1n∑j=1m1∥δtu^Nj−1/2∥2+τ∑k=0n+1∥f^k∥2),

bk=τkσm2−3+τ2−αkσm1−2−α+τ1−αn−α+τ.

Based on the condition σ_m₁, σ_m₂ ≤ 3, we have

∑k=1nbk≤CT2−α.(3.13)

Theorem 3.1

Suppose thatuNn+1 (n = 0, 1, …, K − 1) are the solutions of numerical scheme (2.10). F(u) satisfies the local Lipschitz condition, C, C₀, C₁andĈare appropriate positive constants independent of n, τ, andN. Suppose∥u^Nk+1∥∞ ≤ C₁under the condition∥u^Nk∥∞ ≤ C₀ (0 ≤ k ≤ n). Ifσ_m₁, σ_m₂ ≤ 3, ρn≤(C0C^N)2exp(−CT2−α),then

∥u^Nn+1∥2≤ρnexp(CT2−α).(3.14)

Proof

We use the mathematical induction for performing this proof. The inequality (3.14) can be easily proved for n = 0. Then we verify that the inequality (3.14) holds for any 0 ≤ n ≤ K − 1. Assume that ∥u^Nk∥∞ ≤ C₀ (0 ≤ k ≤ m),

∥u^Nn+1∥2≤ρn+C∑k=1nbk∥u^Nk∥2,0≤n≤m.(3.15)

From Gronwall inequality (see Lemma 3.2),

∥u^Nn+1∥2≤ρnexp(CT2−α),0≤n≤m.(3.16)

Next, we prove that (3.14) holds for n = m + 1, by Lemma 3.5

∥u^Nm+1∥∞2≤(C^N)2∥u^Nm+1∥2≤(C^N)2ρmexp(CT2−α)≤(C0)2.(3.17)

We get ∥u^Nk∥∞ ≤ C₀ (0 ≤ k ≤ m + 1). Thus, it infers that,

∥u^Nm+2∥2≤ρm+1+C∑k=1m+1bk∥u^Nk∥2≤ρm+1exp(CT2−α).(3.18)

The inequality (3.14) is held for n = m + 1.

The proof is completed. □

Next, the convergence analysis of the numerical method can be presented.

Theorem 3.2

Assume thatr ≥ 2, let uanduNn+1 (0 ≤ n ≤ K − 1) be the solutions of(1.1)and(2.10), respectively. Ifr ≤ N + 1, u ∈ C³(I; H^r(Ω)), F(u) ∈ L²(I; H^r(Ω)) satisfies the local Lipschitz condition, u₀ ∈ H^r(Ω), andm₁, m₂ ≤ mwithσ_m₁, σ_m₂ ≤ 3, for small enoughτ, there exists a positive constantCindependent ofn, τ, andNsuch that

∥uNn+1−u(⋅,tn+1)∥≤C(N−r+τmin{2,σm1+1+0.5−α,σm2+1−0.5}).(3.19)

Proof

Setting u_* = ΠN1,0u,e = u_N − u_*andη = u_* − u. We can obtain the following error equation:

(δtm2ek+1/2,v)+(Dτα,k+1/2,m1e,v)+V(∇ek+1/2,v)+D(∇ek+1/2,∇v)=(F(uNk+1/2)−F(u∗k+1/2),v)+Hk,v,(3.20)

where

Hn=−δtm2ηk+1/2−Dτα,k+1/2,m1η−V∇ηk+1/2+F(u∗k+1/2)−F(uk+1/2)+O(τ2tnσm1+1−2−α)+O(τ2tnσm2+1−3).

Taking v = e^k+1/2. Similar to the process of Theorem 3.1,

∥en+1∥2≤C(∑j=1m1∥δtej−1/2∥2+∑j=1m2∥δtej−1/2∥2)+Cτ∑k=1n∥Hk∥2.(3.21)

From Lemma 3.1, Lemma 3.4, and Lemma 3.6

∥Hk∥≤C(N−r+τσm1+1−αkσm1+1−2−α+τσm2+1−1kσm2+1−3).

Based on the above results, we obtain

∥en+1∥2≤C(∑j=1m2∥δtej−1/2∥2+∑j=1m1∥δtu^Nj−1/2∥2+N−2r+τ2σm1+1−2α+1∑k=1nk2(σm1+1−2−α)+τ2σm2+1−1∑k=1nk2(σm2+1−3)).(3.22)

Moreover,

∑j=1m1∥δtej−1/2∥2+∑j=1m2∥δtej−1/2∥2≤C(N−2r+τmin{4,2σm1+1+1−2α,2σm2+1−1}).

It can be derived that

∥en+1∥≤C(N−r+τmin{2,σm1+1+0.5−α,σm2+1−0.5}).

Therefore, (3.19) is proven, which ends the proof. □

4 Fast method

In this section, we introduce a fast method for approximating the discrete convolution for discretizing the time fractional derivative to reduce the memory requirement and computational cost.

The convolution weights λn(α) is expressed by

λn(α)=τ1+αsinαππ∫0∞wα(1+wτ)−1−nF^λ(−w)dw,(4.1)

where

F^λ(w)=(τw)−αλ(α)(1−τw),(4.2)

in which λ^(α) are the generalized formula of order p defined by

λ(α)(z)=(1−z)α∑l=12g^l−1(1−z)l−1=∑l=0∞λl(α)zl,(4.3)

the values of ĝ_i (i = 0, 1, …, 5) can be seen in [8]. Let w = exp(z) [29], Eq. (4.1) becomes

λn(α)=τ1+α∫−∞∞ψn(z)dz,(4.4)

where

ψn(z)=(1+ezτ)−1−nψ(z),ψ(z)=sin(απ)πe(1+α)zF^λ(−ez).(4.5)

It can be observed that ψ_n(z) decays exponentially as |z| → ∞. Using the exponentially convergent trapezoidal rule [29] to approximate the integral ∫−∞∞ψn(z)dz,

λn(α)=τ1+α∑j=−∞∞pj(1+wjτ)−1−n,(4.6)

where w_j = e^z_j, z_j = jΔ z, Δ z > 0, Q > 1 is a positive number, p_j = Δ zψ(z_j). Determination of p_j and w_j can be found in [10]. We derive the following equation by truncating (4.6),

λn(α)=τ1+α∑j=0Q−1pj(1+wjτ)−1−n.(4.7)

The discrete convolution 1τα∑l=0nλn−l(α)ul is decomposed into

1τα∑l=0nλn−l(α)ul=1τα∑l=n−n0nλn−l(α)ul+1τα∑l=0n−n0−1λn−l(α)ul,(4.8)

where 1τα∑l=n−n0nλn−l(α)ul can be computed directly. For 1τα∑l=0n−n0−1λn−l(α)ul, we can have

1τα∑l=0n−n0−1λn−l(α)ul=τ∑l=0n−n0−1ul∑j=0Q−1pj(1+wjτ)−1−(n−l)=(1+wjτ)−n0−1∑j=0Q−1pjqn−n0(j),(4.9)

where qn−n0(j)=τ∑l=0n−n0−1(1+wjτ)−(n−n0−l)ul satisfies

qn(j)=11+wjτqn−1(j)+τun−1,q0(j)=0.(4.10)

Therefore, the fully discrete form using fast method can be obtained,

(δtm2uNn+1/2,v)+12τα(∑l=n−n0nλn−l(α)+∑l=n−n0n+1λn+1−l(α))(ul,v)+(1+wjτ)−n0−12(∑j=0Q−1pjqn−n0(j),v)+(1+wjτ)−n0−12(∑j=0Q−1pjqn+1−n0(j),v)+V(∇uNn+1/2,v)+D(∇uNn+1/2,∇v)=F(uNn+1/2),v+(fn+1/2,v)−12τα∑j=1m1(Wn,jα+Wn+1,jα)(uj−u0,v).(4.11)

The memory requirement and computational cost are O(Q) and O(QK).

5 Numerical experiments

In this section, we perform numerical experiments to verify the effectiveness and robustness of our methods.

The convergence orders in time and space in the L²-norm sense are defined as

order=log(∥error(τ1)∥/∥error(τ2)∥)log(τ1/τ2)intime,log(∥error(N1)∥/∥error(N2)∥)log(N2/N1)inspace,(5.1)

where error=∥u(tn)−uNn∥ is the error equation, τ₁ ≠ τ₂, and N₁ ≠ N₂.

Example 5.1

The following equation with non-smooth solution is presented,

ut+0CDt0.1u=−∇u+Δu+u2+f,(x,y,t)∈Ω×I,(5.2)

where Ω = (0, 2)², I = (0, 1], and

f=∑k=16((1+0.1(k−1))t0.1(k−1)k+1+Γ(2+0.1(k−1))t1+0.1(k−2)(k+1)Γ(2+0.1(k−2))+2π2t1+0.1(k−1)k+1)sinπxsinπy−∑k=16t1+0.1(k−1)k+1sinπxsinπy+π∑k=16t1+0.1(k−1)k+1(cosπxsinπy+sinπxcosπy),(5.3)

with the following initial-boundary conditions

u(x,y,0)=0,(x,y)∈Ω,u(⋅,t)=0,(x,y)∈∂Ω,t∈I.(5.4)

The exact solution of this problem is u=∑k=16t1+0.1(k−1)k+1sinπxsinπy.

Figure 5.1 shows the numerical solution and the exact solution of Example 5.1 when N = 50 and τ = 0.01, the numerical solution fits well with the exact solution. According to the regularity of solution, we choose N = 32, m₁ = m₂ = m and σ_j = 1 + 0.1(j − 1), (j = 1, 2,…). Table 5.1 and Table 5.2 exhibit the L² errors of the method (2.10). We can see that the adding of correction terms increases the accuracy and second-order (even higher) accuracy when three correction terms are used. By selecting the correction term appropriately, we can give the numerical solution with higher accuracy.

Figure 5.1

The numerical solution and the exact solution for Example 5.1.

Table 5.1

The errors and order of L²-norm versus τ for Example 5.1.

1/τ	m = 0		m = 1		m = 3
	Error	Order	Error	Order	Error	Order
10	4.4286e-04	–	5.5752e-04	–	4.9461e-04	–
20	1.4475e-04	1.6133	1.3051e-04	2.0949	8.9802e-05	2.4615
40	1.8065e-04	-0.3196	4.5476e-05	1.5210	1.6803e-05	2.4180
80	1.7468e-04	0.0485	2.4136e-05	0.9139	2.5634e-06	2.7126
160	1.6002e-04	0.1264	1.6798e-05	0.5229	3.1461e-07	3.0265

Table 5.2

The errors and order of L^∞-norm versus τ for Example 5.1.

1/τ	m = 0		m = 1		m = 3
	Error	Order	Error	Order	Error	Order
10	5.3630e-04	–	7.6964e-04	–	6.3155e-04	–
20	2.1836e-04	1.2963	1.8287e-04	2.0733	1.1189e-04	2.4968
40	2.6904e-04	-0.3011	6.4215e-05	1.5099	2.0676e-05	2.4360
80	2.5587e-04	0.0724	3.5178e-05	0.8682	3.5882e-06	2.5266
160	2.3311e-04	0.1344	2.4467e-05	0.5238	5.6949e-07	2.6555

Let Q = 256 and α = 0.3, we apply the fast method to solve the problem. Figure 2(a) displays the computational time of the fast method (4.11) and the direct method (2.10). Obviously, the computational time of the fast method increases linearly, while the computational time of the fast method increases superlinearly (quadratic complexity). Figure 2(b) shows the difference of the fast method solution and the direct method solution. We observe that machine precision is observed. It is also shown that the error caused by the trapezoidal rule in the fast method is independent of the time stepsize τ and the order of the polynomial space N, which is also almost independent of the fractional order α.

Figure 5.2

(a): The computational time of the fast method and direct method for Example 5.1; (b): The difference between the numerical solutions of the fast method and direct method for Example 5.1.

Example 5.2

Consider the following equation

ut+0CDtαu=−∇u+Δu+u−u3+f,(x,y,t)∈Ω×I,(5.5)

with the following initial-boundary conditions

u(x,y,0)=0,(x,y)∈Ω,u(0,y,t)=∑k=04tkey,u(1,y,t)=∑k=04tke1+y,y∈∂Ω,t∈I,u(x,0,t)=∑k=04tkex,u(x,1,t)=∑k=04tkex+1,x∈∂Ω,t∈I.(5.6)

In this example, we choose

f=∑k=14ktk−1+Γ(k+1)tk−αΓ(k+1−α)−tkex+y+(∑k=04tkex+y)3

such that the exact solution of this problem is u=∑k=04tkex+y, and I = (0, 1].

The errors and order of L²-norm versus τ for Example 5.2.

Taking α = 0.3, τ = 0.01 and N = 32. The numerical solution and the exact solution of Example 5.2 are shown in Fig. 3. From the convergence analysis of Theorem 3.2, if the solution u(x, y, t) is sufficiently smooth in time, the convergence rate is O(τ²) by considering m₁ = 2 and m₂ = 0. When m₁ = 0, the convergence rate is O(τ^1.5−α); when m₁ = 1, the convergence rate is O(τ^2.5−α) for α > 1/2, O(τ²) for α ≤ 1/2. Table 2 displays the L²-errors at t = 1 when m₂ = 0. It can be seen that second-order accuracy is observed for m₁ = 2 or m₁ = 1 (α ≤ 1/2). When m₁ = 0, we can not obtain second-order accuracy. These numerical results are agree with the above theoretical analysis.

Figure 5.3

The numerical solution and the exact solution for Example 5.2.

Similar to Example 5.1, let Q = 256 and α = 0.3. The numerical solution is obtained based on the fast method. Fig. 4(a) presents the computational time of the fast method and direct method. Fig. 4(b) gives the difference between the fast method and the direct method for different α. The validity of this fast method is further illustrated.

Figure 5.4

(a): The computational time of the fast method and direct method for Example 5.2; (b): The difference between the numerical solutions of the fast method and direct method for Example 5.2.

6 Conclusions

In this paper, we propose a numerical method for the two-dimensional nonlinear time fractional mobile/immobile advection-diffusion equation.

Table 5.3

The errors and order of L₂-norm versus τ for Example 5.2.

m₁	1/τ	α = 0.2		α = 0.5		α = 0.9
		Error	Order	Error	Order	Error	Order
	20	1.5857e-03	–	2.9847e-03	–	3.4156e-03	–
	40	5.0423e-04	1.6530	8.0024e-04	1.4228	1.4503e-03	1.2358
0	80	1.4490e-04	1.7990	2.6358e-04	1.6022	5.8713e-04	1.3046
	160	4.1056e-05	1.8194	7.9546e-05	1.7284	2.1351e-04	1.4594
	320	1.1090e-05	1.8883	2.3088e-05	1.7846	7.3409e-05	1.5403

	20	8.8367e-04	–	1.0844e-03	–	2.3154e-03	–
	40	2.6639e-04	1.7300	3.4332e-04	1.6593	7.7491e-04	1.5792
1	80	7.6374e-05	1.8024	1.0032e-04	1.7749	2.4815e-04	1.6428
	160	2.0547e-05	1.8942	2.7740e-05	1.8546	7.6899e-05	1.6902
	320	5.3126e-06	1.9514	7.2367e-06	1.9386	2.2647e-05	1.7636

	20	9.2044e-04	–	1.0625e-03	–	2.9242e-03	–
	40	2.3686e-04	1.9583	2.8836e-04	1.8815	7.8566e-04	1.8961
2	80	5.8651e-05	2.0138	7.3227e-05	1.9774	2.0268e-04	1.9547
	160	1.3883e-05	2.0788	1.8218e-05	2.0070	5.1397e-05	1.9794
	320	3.2168e-06	2.1096	4.4093e-06	2.0228	1.3539e-05	1.9900

The WSGD difference scheme is developed for the time stepping, while we apply the Legendre spectral method for the space discretization. The correction scheme is introduced to deal with the non-smooth solutions, and the stability and convergence analysis are proven. In the numerical implementation, a fast method is used based on a globally uniform approximation of the trapezoidal rule for the integral on the real line to decrease the memory requirement and computational cost. Some numerical results are provided to confirm our theoretical analysis and the effectiveness of the presented methods.

Acknowledgements

We would like to express our gratitude to the Editor for taking time to handle the manuscript. This work has been supported by the National Natural Science Foundation of China (Grants Nos. 12001326, 11771254, 11801221), Natural Science Foundation of Jiangsu Province (Grant No. BK20180586), Natural Science Foundation of Shandong Province (Grants No. ZR2020QA032, ZR2019ZD42), China Postdoctoral Science Foundation (Grant Nos. BX20190191, 2020M672038), and Australian Research Council via the Discovery Projects (DP 180103858, DP 190101889).

References

[1] L. Banjai, M. López-Fernández, Efficient high order algorithms for fractional integrals and fractional differential equations. Numer. Math. 141, No 2 (2019), 289–317; 10.1007/s00211-018-1004-0.Suche in Google Scholar

[2] D.A. Benson, M.M. Meerschaert, A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations. Adv. Water. Resources32 (2009), 532–539; 10.1016/j.advwatres.2009.01.002.Suche in Google Scholar

[3] F. Chen, J. Shen, H.J. Yu, A new spectral element method for pricing European options under the Black-Scholes and Merton jump diffusion models. J. Sci. Comput. 52 (2012), 499–518; 10.1007/s10915-011-9556-5.Suche in Google Scholar

[4] K. Diethelm, The Analysis of Fractional Differential Equations. Springer-Verlag, Berlin (2010).10.1007/978-3-642-14574-2Suche in Google Scholar

[5] K. Diethelm, J.M. Ford, N.J. Ford, M. Weilbeer, Pitfalls in fast numerical solvers for fractional differential equations. J. Comput. Appl. Math. 186, No 2 (2006), 482–503; 10.1016/j.cam.2005.03.023.Suche in Google Scholar

[6] W. Fan, F. Liu, X. Jiang, I. Turner, A novel unstructured mesh finite element method for solving the time-space fractional wave equation on a two-dimensional irregular convex domain. Fract. Calc. Appl. Anal. 20, No 2 (2017), 352–383; 10.1515/fca-2017-0019; https://www.degruyter.com/view/journals/fca/20/2/fca.20.issue-2.xml.Suche in Google Scholar

[7] L. Feng, F. Liu, I. Turner, L. Zheng, Novel numerical analysis of multi-term time fractional viscoelastic non-Newtonian fluid models for simulating unsteady and MHD and Couette flow of a generalized Oldroyd-B fluid. Fract. Calc. Appl. Anal. 21, No 4 (2018), 867–868; 10.1515/fca-2018-0058; https://www.degruyter.com/view/journals/fca/21/4/fca.21.issue-4.xml.Suche in Google Scholar

[8] L. Galeone, R. Garrappa, Fractional Adams-Moulton methods. Math. Comput. Simulation79 (2008), 1358–1367; 10.1016/j.matcom.2008.03.008.Suche in Google Scholar

[9] B.Y. Guo, T.J. Wang, Composite Laguerre-Legendre spectral method for fourth-order exterior problems. J. Sci. Comput. 44 (2010), 255–285; 10.1007/s10915-010-9367-0.Suche in Google Scholar

[10] L. Guo, F.H. Zeng, I. Turner, K. Burrage, G.E. Karniadakis, Efficient multistep methods for tempered fractional calculus: Algorithms and simulations. SIAM J. Sci. Comput. 41, No 4 (2019), A2510–A2535; 10.1137/18M1230153.Suche in Google Scholar

[11] P. Guo, C. Zeng, C. Li, Y. Chen, Numerics for the fractional Langevin equation driven by the fractional Brownian motion. Fract. Calc. Appl. Anal. 16, No 1 (2013), 123–141; 10.2478/s13540-013-0009-8; https://www.degruyter.com/view/journals/fca/16/1/fca.16.issue-1.xml.Suche in Google Scholar

[12] M. Ilic, F. Liu, I. Turner, V. Anh, Numerical approximation of a fractional-in-space diffusion equation (I). Fract. Calc. Appl. Anal. 8, No 3 (2005), 323–341; at http://www.math.bas.bg/complan/fcaa.Suche in Google Scholar

[13] M. Ilic, F. Liu, I. Turner, V. Anh, Numerical approximation of a fractional-in-space diffusion equation (II): With nonhomogeneous boundary conditions. Fract. Calc. Appl. Anal. 9, No 4 (2006), 333–349; at http://www.math.bas.bg/complan/fcaa.Suche in Google Scholar

[14] H. Jafari, H. Tajadodi, D. Baleanu, A modified variational iteration method for solving fractional Riccati differential equation by Adomian polynomials. Fract. Calc. Appl. Anal. 16, No 1 (2013), 109–122; 10.2478/s13540-013-0008-9; https://www.degruyter.com/view/journals/fca/16/1/fca.16.issue-1.xml.Suche in Google Scholar

[15] R. Lazarov, P. Vabishchevich, A numerical study of the homogeneous elliptic equation with fractional boundary conditions. Fract. Calc. Appl. Anal. 20, No 2 (2017), 337–351; 10.1515/fca-2017-0018; https://www.degruyter.com/view/journals/fca/20/2/fca.20.issue-2.xml.Suche in Google Scholar

[16] C. Li, Q. Yi, A. Chen, Finite difference methods with non-uniform meshes for nonlinear fractional differential equations. J. Comput. Phys. 316 (2016), 614–631; 10.1016/j.jcp.2016.04.039.Suche in Google Scholar

[17] C. Li, F. Zeng, F. Liu, Spectral approximations to the fractional integral and derivative. Fract. Calc. Appl. Anal. 15, No 3 (2012), 383–406; 10.2478/s13540-012-0028-x; https://www.degruyter.com/view/journals/fca/15/3/fca.15.issue-3.xml.Suche in Google Scholar

[18] C.P. Li, Q. Yi, Finite difference method for two-dimensional nonlinear time-fractional subdiffusion equation. Fract. Calc. Appl. Anal. 21, No 4 (2018), 1046–1072; 10.1515/fca-2018-0057; https://www.degruyter.com/view/journals/fca/21/4/fca.21.issue-4.xml.Suche in Google Scholar

[19] J.R. Li, A fast time stepping method for evaluating fractional integrals. SIAM J. Sci. Comput. 31 (2010), 4696–4714; 10.1137/080736533.Suche in Google Scholar

[20] H.L. Liao, T. Tang, T. Zhou, A second-order and nonuniform time-stepping maximum-principle preserving scheme for time-fractional Allen-Cahn equations. J. Comput. Phys. 414 (2020), ID 109473; 10.1016/j.jcp.2020.109473.Suche in Google Scholar

[21] F. Liu, M.M. Meerschaert, R. McGough, P. Zhuang, Q. Liu, Numerical methods for solving the multi-term time fractional wave equations. Fract. Calc. Appl. Anal. 16, No 1 (2013), 9–25; 10.2478/s13540-013-0002-2; https://www.degruyter.com/view/journals/fca/16/1/fca.16.issue-1.xml.Suche in Google Scholar PubMed PubMed Central

[22] Q. Liu, F. Liu, I. Turner, V. Anh, Y.T. Gu, A RBF meshless approach for modeling a fractal mobile/immobile transport model. Appl. Math. Comput. 226 (2014), 336–347; 10.1016/j.amc.2013.10.008.Suche in Google Scholar

[23] Y. Liu, Y.W. Du, H. Li, J.F. Wang, A two-grid finite element approximation for a nonlinear time-fractional Cable equation. Nonlinear Dyn. 85 (2016), 2535–2548; 10.1007/s11071-016-2843-9.Suche in Google Scholar

[24] M. López-Fernández, C. Lubich, A. Schädle, Adaptive, fast, and oblivious convolution in evolution equations with memory. SIAM J. Sci. Comput. 30 (2008), 1015–1037; 10.1137/060674168.Suche in Google Scholar

[25] C. Lubich, Discretized fractional calculus. SIAM J. Math. Anal. 17 (1986), 704–719; 10.1137/0517050.Suche in Google Scholar

[26] Y. Luchko, Initial-boundary-value problems for the generalized multi-term time fractional diffusion equation. J. Math. Anal. Appl. 374 (2011), 538–548; 10.1016/j.jmaa.2010.08.048.Suche in Google Scholar

[27] H.P. Ma, Y.B. Yang, Jacobi spectral collocation method for the time variable-order fractional mobile-immobile advection-dispersion solute transport model. East Asian J. Appl. Math. 6 (2016), 337–352; 10.4208/eajam.141115.060616a.Suche in Google Scholar

[28] W. McLean, K. Mustapha, A second-order accurate numerical method for a fractional wave equation. Numer. Math. 105 (2007), 481–510; 10.1007/s00211-006-0045-y.Suche in Google Scholar

[29] W. McLean, Exponential Sum Approximations for t^−β. In: J. Dick, F. Kuo, H. Woźniakowski (Eds). Contemporary Computational Mathematics-A Celebration of the 80th Birthday of Ian Sloan. Springer, Cham (2018).Suche in Google Scholar

[30] T.B. Nguyen, B. Jang, A high-order predictor-corrector method for solving nonlinear differential equations of fractional order. Fract. Calc. Appl. Anal. 20, No 2 (2017), 447–476; 10.1515/fca-2017-0023; https://www.degruyter.com/view/journals/fca/20/2/fca.20.issue-2.xml.Suche in Google Scholar

[31] I. Podlubny, Fractional Differential Equations. Academic Press, San Diego (1999).Suche in Google Scholar

[32] A. Schädle, M. López-Fernández, C. Lubich, Fast and oblivious convolution quadrature; SIAM J. Sci. Comput. 28 (2006), 421–438; 10.1137/050623139.Suche in Google Scholar

[33] R. Schumer, D.A. Benson, M.M. Meerschaert, B. Baeumer, Fractal mobile/immobile solute transport. Water Resour. Res. 39, (2003), 1296-1307; 10.1029/2003WR002141.Suche in Google Scholar

[34] J. Shen, T. Tang, L.L. Wang, Spectral Methods: Algorithms, Analysis and Applications. Springer Ser. Comput. Math. Springer, Heidelberg (2011).10.1007/978-3-540-71041-7Suche in Google Scholar

[35] M. Stynes, E. O'Riordan, J. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal. 55 (2017), 1057–1079; 10.1137/16m1082329.Suche in Google Scholar

[36] L.L. Wang, Analysis of spectral approximations using prolate spheroidal wave functions. Math. Comput. 79 (2010), 807–827; 10.1090/S0025-5718-09-02268-6.Suche in Google Scholar

[37] F. Zeng, F. Liu, C. Li, K. Burrage, I. Turner, V. Anh, Crank-Nicolson ADI spectral method for the two-dimensional Riesz space fractional nonlinear reaction-diffusion equation. SIAM J. Numer. Anal. 52 (2014), 2599–2622; 10.1137/130934192.Suche in Google Scholar

[38] F. Zeng, Z. Zhang, G.E. Karniadakis, Second-order numerical methods for multi-term fractional differential equations: Smooth and non-smooth solutions. Comput. Methods Appl. Mech. Engrg. 327 (2017), 478–502; 10.1016/j.cma.2017.08.029.Suche in Google Scholar

[39] F. Zeng, I. Turner, K. Burrage, A stable fast time-stepping method for fractional integral and derivative operators. J. Sci. Comput. 77 (2018), 283–307; 10.1007/s10915-018-0707-9.Suche in Google Scholar

[40] H. Zhang, F. Liu, M.S. Phanikumar, M.M. Meerschaert, A novel numerical method for the time variable fractional order mobile-immobile advection-dispersion model. Comput. Math. Appl. 66 (2013), 693–701; 10.1016/j.camwa.2013.01.031.Suche in Google Scholar

[41] J. Zhang, C.J. Xu, Finite difference/spectral approximations to a water wave model with a nonlocal viscous term. Appl. Math. Model. 38 (2014), 4912–4925; 10.1016/j.apm.2014.03.051.Suche in Google Scholar

[42] Y. Zhang, D.A. Benson, D.M. Reeves, Time and space nonlocalities underlying fractional-derivative models: Distinction and literature review of field applications. Adv. Water. Resources32 (2009), 561–581; 10.1016/j.advwatres.2009.01.008.Suche in Google Scholar

Received: 2019-10-27

Revised: 2020-11-13

Published Online: 2021-01-29

Published in Print: 2021-02-23

Artikel in diesem Heft

https://doi.org/10.1515/fca-2021-0009

Schlagwörter für diesen Artikel

two-dimensional nonlinear time fractional mobile/immobile advection-dispersion equation; correction method; stability and convergence analysis; fast method