A numerical method for distributed order time fractional diffusion equation with weakly singular solutions

A numerical method for distributed order time fractional diffusion equation with weakly singular solutions

Applied Mathematics Letters 96 (2019) 159–165 Contents lists available at ScienceDirect Applied Mathematics Letters www.elsevier.com/locate/aml A n...

701KB Sizes 0 Downloads 116 Views

Applied Mathematics Letters 96 (2019) 159–165

Contents lists available at ScienceDirect

Applied Mathematics Letters www.elsevier.com/locate/aml

A numerical method for distributed order time fractional diffusion equation with weakly singular solutions Jincheng Ren a ,1 , Hu Chen b ,∗,2 a

College of Mathematics and Information Science, Henan University of Economics and Law, Zhengzhou 450045, China b School of Mathematics and Information Science, Yantai University, Yantai 264005, China

article

info

Article history: Received 28 January 2019 Received in revised form 27 April 2019 Accepted 27 April 2019 Available online 15 May 2019 Keywords: Distributed order derivative L1 scheme Weak singularity Convergence analysis

abstract A finite difference/spectral method is proposed for the numerical approximation of a distributed order time fractional diffusion equation with initial singularity on two dimensional spatial domain. The L1 scheme on graded mesh for the discretization of time Caputo fractional derivative is used to capture the weak initial singularity and Legendre spectral method is adopted for the spatial discretization. It is proved that with appropriate choice of the grading parameter, the scheme can attain order 2 − β convergence in temporal direction, where β (0 < β < 1) is the upper bound of the order of distributed order fractional derivative. Numerical results confirm the sharpness of the error analysis. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction In this paper, we consider the following distributed order time-fractional diffusion equation in the spatial domain Ω = (−1, 1) × (−1, 1) ∈ R2 : Dtw u(x, t) − ν∆u(x, t) + c(x)u(x, t) = f (x, t)

for x := (x1 , x2 ) ∈ Ω , 0 < t ≤ T,

(1.1a)

subject to the initial and boundary conditions for x ∈ Ω ,

(1.1b)

for 0 ≤ t ≤ T,

(1.1c)

u(x, 0) = u0 (x) u(x, t)|∂Ω = 0

here ν > 0 is a constant diffusion coefficient, ∆ is the two dimensional Laplacian operator, f is a given ¯ ), and Dtw u(x, t) denotes the distributed order fractional derivative of u in time t, function, c(x) ≥ 0 ∈ C(Ω ∗

Corresponding author. E-mail addresses: [email protected] (J. Ren), [email protected] (H. Chen). 1 The research of this author is supported in part by the NSF of China (No. 11601119), sponsored by Program for HASTIT (No. 18HASTIT027) and Young talents Fund of HUEL. 2 The research of this author is supported in part by the Chinese Postdoc Foundation (No. 2018M631316) and the NSF of China (No. 11801026). https://doi.org/10.1016/j.aml.2019.04.030 0893-9659/© 2019 Elsevier Ltd. All rights reserved.

160

J. Ren and H. Chen / Applied Mathematics Letters 96 (2019) 159–165

given by Dtw u(x, t) =



β

w(α)Dtα u(x, t)dα,

0 < β ≤ 1,

(1.2)

0

∫β where w(α) ≥ 0, 0 w(α)dα = c0 > 0, Dtα u(x, t) (0 < α < 1) is defined as the Caputo fractional derivative of order α with respect to t, defined by ∫ t ∂u(x, s) ds 1 , (1.3) Dtα u(x, t) = Γ (1 − α) 0 ∂s (t − s)α in which Γ (·) is Euler’s Gamma function. Due to the excellent modeling capability of Eqs. (1.1a)–(1.1c), its accurate numerical treatment has been the subject of numerous studies; see [1,2] and the references therein. However, most of the previous published work are concerned with globally smooth solutions. Recently, it is well-known that for time fractional derivative problems, the solutions usually exhibit some weak singularity at t = 0, i.e, although the solution is continuous at t = 0, its temporal first and higher order derivatives blow up as t → 0+ . For single term time fractional diffusion equation, L1 scheme on graded meshes is used to tackle this weak initial singularity, see [3,4]. At present, it still leaves the question how to extend the L1 scheme to distributed order time fractional diffusion equation with non-smooth solution. In this work, we develop finite difference/spectral method for problem (1.1a)–(1.1c) whose solution has some weak singularity at initial time t = 0, and establish rigorous error estimates for the fully discrete scheme. Our error analysis is based on the ideas of [5] and is performed on graded meshes. Numerical results show the sharpness of the result. Throughout the paper, we use A ≲ B (A ≳ B) to mean that A ≤ cB (A ≥ cB), where c is a generic positive constant, and A ≃ B to mean that A ≲ B and A ≳ B. 2. The L1 scheme for time semidiscrete system We introduce the following composite midpoint formula for numerical integration, this formula can be found in [6, (5.1.19)]. Lemma 2.1. Let s(α) ∈ C 2 [0, β]. Divide the interval [0, β] into J-subintervals with hα = β/J and αl = (l − 1/2)hα , l = 1, . . . , J. Then we have ∫

β

s(α)dα = hα 0

J ∑

s(αl ) +

l=1

βh2α ′′ s (ξ), 24

ξ ∈ [0, β].

∑J αl 2 Denote α := (α1 , . . . , αJ ), and Dα t u := hα l=1 w(αl )Dt u. Thus, under the conditions w(α) ∈ C [0, β], Dtα u(·, ·) ∈ C 2 [0, β], the initial equation (1.1a) is equivalent to 2 Dα t u(x, t) − ∆u(x, t) + c(x)u(x, t) = f (x, t) + O(hα ).

(2.1)

For a positive integer M , considering the graded mesh {tj = T (j/M )r }M j=0 with some r ≥ 1. For k ≥ 1, denote ∫ tj τj−1 (α) dk,j = (tk − η)−α dη for j = 1, . . . , k. (2.2) Γ (1 − α) tj−1 (α)

(α)

It is easy to see that dk,j+1 > dk,j > 0 for all admissible k and j. For any mesh function {V j }M j=0 , define (α)

δtα V k := dk,k V k −

k−1 ∑ j=0

(α)

(α)

(dk,j+1 − dk,j )V j

for k = 1, . . . , M ,

(2.3)

J. Ren and H. Chen / Applied Mathematics Letters 96 (2019) 159–165

(α)

where dk,0 := 0. With dα k,j := hα k Lα t V := hα

J ∑

161

(α )

∑J

w(αl )dk,jl , we then define

l=1 α

k w(αl )δt l V k = dα k,k V −

k−1 ∑

α j (dα k,j+1 − dk,j )V

for k = 1, . . . , M .

(2.4)

j=0

l=1

α We also have dα k,j > dk,j−1 > 0 for all admissible k and j. M For any mesh function {V j }M j=0 on {tj }j=0 , one has ⎧( ⎫ )−1 J ⎨ ⎬ ∑ −α j |V k | ≤ |V 0 | + Γ (1 − αJ ) max hα w(αl )tj l Lα for k = 1, . . . , M. t |V | ⎭ j=1,...,k ⎩

Lemma 2.2.

l=1

Proof . Suppose maxj=1,...,k |V j | = |V n | for some n ∈ {1, . . . , k}. Set W j = |V j | − |V 0 | for all j. Then W j = 0, and ∑ ( j ) ( n ) n−1 α 0 n α 0 (dα Lα |V | = d W + |V | − n,j+1 − dn,j ) W + |V | t n,n j=0 n = dα n,n W −

n−1 ∑

⎛ α j ⎝dα (dα n,j+1 − dn,j )W + n,n −

j=1 n ≥ dα n,n W −

n−1 ∑

n−1 ∑

⎞ α ⎠ 0 (dα n,j+1 − dn,j ) |V |

j=0

( n ) α 0 (dα n,j+1 − dn,j ) |V | − |V |

j=1

( n ) 0 = dα n,1 |V | − |V | . −1 α Hence |V n | − |V 0 | ≤ (dα Lt |V n |. The result now follows from dα n,1 ) n,1 ≥ hα

Lemma 2.3.

∑J

l=1

−α

w(αl )tn l /Γ (1 − αl ). □

For j = 1, . . . , M , suppose hα is sufficiently small, we have hα

J ∑

−αl

w(αl )tj

−αJ

≳ min{tj−α1 , tj

}.

l=1

Proof . First, we have hα ∫β w(α)dα = c0 and 0

∑J

l=1

−αl

w(αl )tj



J ∑

−αJ

1 ≥ min{t−α , tj j

∫ w(αl ) =

∑J

l=1

β

w(α)dα + O(h2α ) ≥

0

l=1

providing hα is sufficiently small.

}hα

w(αl ). Then the results follow from c0 2



From Kopteva’s paper [5, Lemma 2.3], we have Lemma 2.4 ([7, Lemma 2.2]). Suppose u(t) ∈ C 2 (0, T ], then (α)

|δtα u(tk ) − Dtα u(tk )| ≲ t−α k { max ψj } j=1,...,k

for k = 1, . . . , M,

(2.5)

where (α)

ψ1

(α)

ψj

= τ1α sup

(

) η 1−α |δt u(t1 ) − u′ (η)| ,

(2.6a)

η∈(0,t1 )

= τj2−α tα j

sup η∈(tj−1 ,tj )

|u′′ (η)| for j = 2, . . . , M.

(2.6b)

162

J. Ren and H. Chen / Applied Mathematics Letters 96 (2019) 159–165

We then have the following lemma. Lemma 2.5.

For any function u(x, ·) ∈ C 2 (0, T ], one has α |Lα t u(tk ) − Dt u(tk )| ≲ hα

J ∑ l=1

Proof . By direct triangular inequality.

−αl

w(αl )tk

(αl )

{ max ψj j=1,...,k

}.

(2.7)



Lemma 2.6 ([7, Lemma 2.3]). Let σ ∈ (0, 1). Assume that |u(l) (t)| ≲ 1 + tσ−l for l = 0, 1, 2, and t ∈ (0, T ]. (α) Recall the definitions of ψj from (2.6). Then (α)

ψj

≲ M − min{rσ,2−α} for j = 1, . . . , M.

3. The fully discrete method Set Λ := (−1, 1). Let N be a positive integer. Denote by PN (Λ) the space of all polynomials of degree less than or equal to N . Set P0N := {ϕ ∈ PN (Λ) : ϕ(±1) = 0} and H01 (Ω ) := {v ∈ H 1 (Ω ) : v|∂Ω = 0}. 1,0 Let πN be the H01 -orthogonal projection operator from H01 (Ω ) into (P0N )2 ; that is, for each v ∈ H01 (Ω ), ( ) 1,0 ∇πN v, ∇vN = (∇v, ∇vN ) , ∀vN ∈ (P0N )2 . (3.1) For this projection operator the following approximation result is valid: Lemma 3.1 ([8, eq.(5.8.15)]). For any v ∈ H01 (Ω ) ∩ H s (Ω ), one has 1,0 v∥k ≤ CN k−s ∥v∥s for 0 ≤ k ≤ 1 and s ≥ 1, ∥v − πN

where C is a positive constant independent of N and v. We discretize in space using a Legendre spectral method. Let ujN ∈ (P0N )2 for j = 0, . . . , M be the approximation of the exact solution u at time t = tj . Then the weak formulation of the fully discrete method 1,0 for problem (1.1) is: find ukN ∈ (P0N )2 , with u0N = πN u0 , such that k k k k 0 2 (Lα t uN , vN ) + ν(∇uN , ∇vN ) + (cuN , vN ) = (f , vN ), ∀vN ∈ (PN ) ,

(3.2)

for k = 1, . . . , M . It is well known that the well-posedness of the discrete problem (3.2) is guaranteed by the Lax–Milgram lemma. 3.1. Error analysis of the fully discrete scheme Next we prove the stability of the method (3.2). Lemma 3.2.

Let {ujN }M j=0 be the solution of (3.2). Then ∥ukN ∥ ≲ ∥u0N ∥ + max ∥f j ∥ 1≤j≤M

for k = 1, . . . , M.

J. Ren and H. Chen / Applied Mathematics Letters 96 (2019) 159–165

163

Proof . For k = 1, . . . , M , taking vN = ukN in (3.2) gives k k k 2 k k k k (Lα t uN , uN ) + ν∥∇uN ∥ + (cuN , uN ) = (f , uN ).

Using the definition (2.4), c ≥ 0 and Cauchy–Schwarz inequality, we get k 2 k 2 dα k,k ∥uN ∥ + ν∥∇uN ∥ ≤

k−1 ∑

j α k k k (dα k,j+1 − dk,j )∥uN ∥ · ∥uN ∥ + ∥f ∥ · ∥uN ∥.

(3.3)

j=0

Then omitting the positive term ν∥∇ukN ∥2 , dividing both sides by ∥ukN ∥, and according to definition of Lα t , k k one gets Lα ∥u ∥ ≤ ∥f ∥. Thus the results follow from Lemma 2.2, and t N ( )−1 J ∑ −αl αJ 1 hα ≲ max{tα w(αl )tj j , tj } ≲ 1 l=1

according to Lemma 2.3. □ For the error analysis of the fully discrete scheme (3.2), we have Theorem 3.1. Let u be the solution of (1.1) and let {ukN }M k=0 be the solution of the fully discrete method (3.2) ∞ s on graded mesh. Assume that u ∈ L∞ (0, T ; H01 (Ω ) ∩ H s (Ω )), that Dα t u ∈ L (0, T ; H (Ω )) for some s ≥ 1, σ−l l for l = 0, 1, 2, σ ∈ (0, 1). Then and that ∥∂t u(t)∥2 ≲ 1 + t ∥u(tk ) − ukN ∥ ≲ M − min{rσ,

2−β+ h2α }

+ N −s + h2α ,

k = 0, . . . , M.

1,0 1,0 Proof . For j = 0, . . . , M , let ejN = u(tj ) − ujN , e˜jN = πN u(tj ) − ujN and eˆjN = u(tj ) − πN u(tj ). Then 1,0 j j j 0 0 0 eN = e˜N + eˆN , and in particular eN = u0 − πN u0 = eˆN and e˜N = 0. Subtracting the differential equation (1.1a) and the fully discrete scheme (3.2) gives the following error equation for all vN ∈ (P0N )2 :

(Lα ˜m em em t e N , vN ) + ν(∇˜ N , ∇vN ) + (c˜ N , vN ) ( α ) α α m w = Lt u(tm ) − Dt u(tm ), vN + (−cˆ em ˆN , vN ) + (Dα N − Lt e t u(tm ) − Dt u(tm ), vN ) =: (rm , vN ) + (P m , vN ) + (Rm , vN ). Taking vN = e˜m N here and using the same argument in Lemma 3.2, one gets m m m Lα em t ∥˜ N ∥ ≤ ∥r ∥ + ∥P ∥ + ∥R ∥.

Then, according to Lemma 2.2, we have ⎧( ⎫ )−1 J ⎨ ⎬ ∑ ( ) l ∥˜ ekN ∥ ≲ max hα w(αl )t−α ∥rm ∥ + ∥P m ∥ + ∥Rm ∥ . m ⎭ 1≤m≤k ⎩ l=1

With the help of Lemma 3.1 gives α ∥P m ∥ = ∥Lα ˆm ˆN (tm ) + Dα ˆN (tm ) + cˆ em t e N − Dt e t e N∥  α ( ) 1,0 α α   + ∥Dα ≤ Lt u(tm ) − Dt u(tm ) − π Lt u(tm ) − Dα ˆN (tm )∥ + ∥cˆ em t u(tm ) t e N∥ N

−s ≲ ∥rm ∥1 + N −s ∥Dα ∥u(tm )∥s . t u∥L∞ (H s ) + N

By virtue of Lemmas 2.5 and 2.6, we have ⎧( ⎫ )−1 {( )} J ⎨ ⎬ ∑ hα (αi ) −αl m max hα w(αl )tm ∥r ∥1 ≲ max max ∥ψj ∥1 ≲ M − min{rσ,2−β+ 2 } . j=1,...,m ⎭ 1≤m≤k 1≤m≤k ⎩ l=1

(3.4)

J. Ren and H. Chen / Applied Mathematics Letters 96 (2019) 159–165

164

Table 4.1 Maximum L2 (Ω ) error and convergence rate with r = (2 − β)/σ. M

64 128 256 512 1024 2048

σ = β = 0.4

σ = β = 0.6

σ = β = 0.8

Error

Rate

Error

Rate

Error

Rate

5.3963e−04 2.0061e−04 7.1851e−05 2.5101e−05 8.6165e−06 2.9210e−06

1.4276 1.4813 1.5173 1.5426 1.5607 *

9.2655e−04 3.9184e−04 1.6224e−04 6.5550e−05 2.6057e−05 1.0237e−05

1.2416 1.2721 1.3074 1.3309 1.3479 *

1.2008e−03 5.9097e−04 2.8328e−04 1.3336e−04 6.1879e−05 2.8395e−05

1.0228 1.0609 1.0869 1.1078 1.1238 *

Table 4.2 Maximum L2 (Ω ) error and convergence rate with r = (2 − β)/σ. M

64 128 256 512 1024 2048

σ = 0.3, β = 0.6

σ = 0.5, β = 0.6

σ = 0.7, β = 0.6

Error

Rate

Error

Rate

Error

Rate

3.0912e−03 1.2454e−03 4.9193e−04 1.9210e−04 7.4363e−05 2.8625e−05

1.3115 1.3401 1.3566 1.3692 1.3773 *

1.4025e−03 5.8849e−04 2.3959e−04 9.5637e−05 3.7658e−05 1.4697e−05

1.2530 1.2964 1.3249 1.3446 1.3574 *

5.7960e−04 2.5009e−04 1.0402e−04 4.2501e−05 1.7059e−05 6.7526e−06

1.2126 1.2656 1.2913 1.3169 1.3370 *

Combining (2.1), from (3.4) one gets ∥˜ ekN ∥ ≲ M − min{rσ,2−β+

hα 2 }

+ N −s + h2α .

ekN ∥ and Lemma 3.1, we get the desired ekN ∥ + ∥ˆ Finally, using the triangle inequality ∥u(tk ) − ukN ∥ ≤ ∥˜ result. □ 4. Numerical results Example 4.1. Consider the problem (1.1) with ν = 0.1, T = 1, c(x) = x21 +x22 , w(α) = Γ (σ+1−α)/Γ (1+σ) and with an exact solution that displays a typical weak singularity at t = 0: u(x, t) = (1 + tσ )(1 − x21 )(1 − x22 ) exp(x1 + x2 ). The right-hand side f (x, t) can be computed from (1.1a). We first check the temporal accuracy of our scheme. Choose N = 20, J = 100 to eliminate the contamination of the spatial error and the error caused by numerical quadrature in distributed variable. Tables 4.1–4.2 present the maximum L2 (Ω ) error and convergence rate for different values of M , σ and β, it can be seen the temporal accuracy can recover the optimal order 2 − β, which conforms with our theoretical prediction. Next we take M = 3000, N = 20 to test the convergence order in distributed variable. From Table 4.3, we can see the convergence rate in the distributed variable is two order, thus justifies our theoretical analysis. At last we check the spatial accuracy with respect to N . We choose M = 3000, J = 100, r = (2 − β)/σ and σ = β = 0.6 as an example. Fig. 1, whose vertical axis is on a log scale, shows that one obtains spectral accuracy in space for this problem whose spatial component is very smooth. Acknowledgments We would like to thank the anonymous referees for many constructive comments and suggestions which led to an improved presentation of this paper.

J. Ren and H. Chen / Applied Mathematics Letters 96 (2019) 159–165

165

Table 4.3 Maximum L2 (Ω ) error and convergence rate with r = (2 − β)/σ and σ = β = 0.6 in the distributed-order variable. hα

Error

Rate

β/2 β/4 β/8 β/16 β/32 β/64

2.5913e−02 6.8078e−03 1.7266e−03 4.3192e−04 1.0663e−04 2.5234e−05

1.9284 1.9793 1.9991 2.0181 2.0792 *

Fig. 1. Error function for Example 4.1 in the case σ = β = 0.6.

References [1] H. Chen, S. L¨ u, W. Chen, Finite difference/spectral approximations for the distributed order time fractional reaction– diffusion equation on an unbounded domain, J. Comput. Phys. 315 (2016) 84–97, http://dx.doi.org/10.1016/j.jcp.2016. 03.044. [2] B.P. Moghaddama, J.A. Tenreiro Machadob, M.L. Morgado, Numerical approach for a class of distributed order time fractional partial differential equations, Appl. Numer. Math. 136 (2019) 152–162, http://dx.doi.org/10.1016/j.apnum.2018. 09.019. [3] 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 (2) (2017) 1057–1079, https://doi.org/10.1137/16M1082329, 3639581. [4] H.-L. Liao, D. Li, J. Zhang, Sharp error estimate of the nonuniform L1 formula for linear reaction–subdiffusion equations, SIAM J. Numer. Anal. 56 (2) (2018) 1112–1133, http://dx.doi.org/10.1137/17M1131829. [5] N. Kopteva, Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions, Math. Comp. (2019) http://dx.doi.org/10.1090/mcom/3410, in press. [6] G. Dahlquist, ˚ A. Bj¨ orck, Numerical Methods in Scientific Computing: Vol. 1, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008. [7] H. Chen, X. Hu, J. Ren, T. Sun, Y. Tang, L1 scheme on graded mesh for the linearized time fractional KdV equation with initial singularity, Int. J. Mod. Sim. Sci. Comp. 10 (2019) 1941006, http://dx.doi.org/10.1142/S179396231941006X. [8] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods: Fundamentals in Single Domains, in: Scientific Computation, Springer-Verlag, Berlin, 2006, p. xxii+563.