Numerical algorithm for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term

# Numerical algorithm for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term

Jiahui Hu Jungang Wang Zhanbin Yuan Zongze Yang Yufeng Nie Research Center for Computational Science, Northwestern Polytechnical University, Xi’an 710129, China College of Science, Henan University of Technology, Zhengzhou 450001, China
###### Abstract

In this paper, an alternating direction implicit (ADI) difference scheme for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term is presented. The unique solvability of the difference solution is discussed, and the unconditional stability and convergence order of the numerical scheme are analysed. Finally, numerical experiments are carried out to verify the effectiveness and accuracy of the algorithm.

###### keywords:
Two-dimensional time-fractional wave equation of distributed-order, ADI scheme, Nonlinear source term, Stability, Convergence
###### Msc:
[2010] 35R11, 65M06, 65M12
\newdefinition

definitionDefinition[section] \newdefinitionexampleExample[section]

## 1 Introduction

The idea of distributed-order differential equation was first introduced by Caputo in his work for modeling the stress-strain behavior of an anelastic medium in 1960s Caputo (). Being different from the differential equations with the single-order fractional derivative and the ones with sums of fractional derivatives, i.e., multi-term fractional differential equations (FDEs), the distributed-order differential equations are derived by integrating the order of differentiation over a certain range. It can be regarded as a generalization of the aforementioned two classes of FDEs. A typical application of this kind of FDEs is in the retarding sub-diffusion process, where a plume of particles spreads at a logarithmic rate, which leads to ultraslow diffusion (see Sinai1983 ()ChechkinKlafterSokolov2003 ()Kochubei2008 ()). Another example is the fractional Langevin equation of distributed-order, which was proposed to model the kinetics of retarding sub-diffusion whose scaling exponent decreases in time, and then was applied to simulate the strongly anomalous ultraslow diffusion with the mean square displacement growing as a power of logarithm of time EabLim2011 (). The distributed-order FDEs were also found playing important role in other various research fields, such as control and signal processing JiaoChenPodlubny2012 (), modelling dielectric induction and diffusion Caputo2001 (), identification of systems Hartley1999 (), and so on.

Till now, there have been many important progresses for the research on analytical solutions of distributed-order FDEs. For the kinetic description of anomalous diffusion and relaxation phenomena, A. V. Chechkin et al. presented the diffusion-like equation with time fractional derivative of distributed-order in Chechkin2003 (), where the positivity of the solutions of the proposed equation was proved and the relation to the continuous-time random walk theory was established. T. M. Atanackovic et al. analysed a Cauchy problem for a time distributed-order diffusion-wave equation by means of the theory of an abstract Volterra equation atanackovic2009time (). In Gorenflo2013 (), for the one-dimensional distributed-order diffusion-wave equation, R. Gorenflo et al. gave the interpretation of the fundamental solution of the Cauchy problem as a probability density function of the space variable evolving in time in the transform domain by employing the technique of the Fourier and Laplace transforms. Using the Laplace transform method, Z. Li et al investigated the asymptotic behavior of solutions to the initial-boundary-value problem for the distributed-order time-fractional diffusion equations LiLuchkoYamamoto2014 ().

In most instances, the analytical solutions of distributed-order differential equations are not easy to available, thus it stimulates researchers to develop numerical algorithms for approximate solutions. To our knowledge, the research on numerically solving the distributed-order differential equations are still in its infancy. The literatures DiethelmFord2009 ()Podlubny2013 ()Katsikadelis2014 () concerned on developing numerical methods for solving distributed-order ordinary differential equations. In terms of the distributed-order partial differential equations, most of the work are about the one-dimensional time distributed-order differential equations, and the integrating range of the order of time derivative is the interval , which is named as time distributed-order diffusion equation. N. J. Ford et al. developed an implicit finite difference method for the solution of the diffusion equation with distributed order in time Ford2015 (). By using the Grünwald-Letnikov formula, Gao et al. proposed two difference schemes to solve the one-dimensional distributed-order differential equations, and the extrapolation method was applied to improve the approximate accuracy GaoSun2016111 (). In GaoSunSun2015 (), the authors handled the same distributed-order differential equations by employing a weighted and shifted Grünwald-Letnikov formula to derive several second-order convergent difference schemes. When the order of the time derivative is distributed over the interval , it is called the time distributed-order wave equation. The study of the numerical solution of this kind of equation is rather more limited. Ye et al. derived and analysed a compact difference scheme for a distributed-order time-fractional wave equation in YeLiuAnh2015 ().

When considering the high-dimensional models, Gao et al. investigated ADI schemes for two-dimensional distributed-order diffusion equations GaoSun2015 ()GaoSun2016 (), and they also developed two ADI difference schemes for solving the two-dimensional time distributed-order wave equations GaoSun2016a (). Due to the widespread use of the nonlinear models RidaEl-SayedArafa2010 ()WazwazGorguis2004 (), M. L. Morgado et al. developed an implicit difference scheme for one-dimensional time distributed-order diffusion equation with a nonlinear source term MorgadoRebelo2015 (). For further discussion on the numerical approaches for solving the high-dimensional distributed-order partial differential equations, this paper is devoted to develop effective numerical algorithm for two-dimensional time-fractional wave equation of distributed-order with a nonlinear source term

 ∫21p(β)C0Dβtu(x,y,t)dβ=∂2u(x,y,t)∂x2+∂2u(x,y,t)∂x2+f(x,y,t,u(x,y,t)), (x,y)∈Ω,t∈(0,T], (1.1) u(x,y,t)=ϕ(x,y,t),(x,y)∈∂Ω,0≤t

where , and is the boundary of . The fractional derivative in (1) is given in the Caputo sense

 C0Dβtv(t)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂v(t)∂t−∂v(0)∂t, β=1,1Γ(2−β)∫t0(t−ξ)1−β∂2v(ξ)∂ξ2dξ, 1<β<2,∂2v(t)∂t2, β=2,

and the function is served as weight for the order of differentiation such that and . We assume that , , , and are continuous, and the nonlinear source term satisfies a Lipschitz condition of the form

 |f(x,y,t,u1)−f(x,y,t,u2)|≤Lf|u1−u2|, (1.4)

where is a positive constant.

The main procedure of developing numerical scheme for solving problem (1)(1.3) is as follows. Firstly a suitable numerical quadrature formula is adopted to discrete the integral in (1), and a multi-term time fractional wave equation is left whereafter. Then we develop an ADI finite difference scheme which is uniquely solvable for the multi-term time fractional wave equation. By using the discrete energy method, we prove the derived numerical scheme is unconditionally stable and convergent.

The rest of this paper is organized in the following way. In Section 2, the ADI finite difference scheme is constructed and described detailedly. In Section 3, we give analysis on solvability, stability and convergence for the derived difference scheme. Numerical results are illustrated in Section 4 to confirm the effectiveness and accuracy of our method, and some conclusions are drawn in the last section.

## 2 The derivation of the ADI scheme

This section focuses on deriving the ADI scheme for the problems (1)(1.3).

Let , and be positive integers, and , and be the uniform sizes of spatial grid and time step, respectively. Then a spatial and temporal partition can be defined as for , for and for . Denote and , then the domain is covered by . Let be a grid function on . We introduce the following notations:

 un−12ij=12(unij+un−1ij),δtun−12ij=1τ(unij−un−1ij),
 δxuni−12,j=1h1(unij−uni−1,j),δ2xunij=1h1(δxuni+12,j−δxuni−12,j),
 δyuni,j−12=1h2(unij−uni,j−1),δ2yunij=1h2(δxuni,j+12−δxuni,j−12),

and

 Δhuij=δ2xuij+δ2yuij.

Consider Eq. (1) at the point , and we write it as

 ∫21p(β)C0Dβtu(xi,yj,tn)dβ (2.1) = ∂2u(xi,yj,tn)∂x2+∂2u(xi,yj,tn)∂y2+f(xi,yj,tn,u(xi,yj,tn)).

Take an average of Eq. (2.1) on time level and , then we have

 12(∫21p(β)C0Dβtu(xi,yj,tn)dβ+∫21p(β)C0Dβtu(xi,yj,tn−1)dβ) (2.2) = +12[f(xi,yj,tn,u(xi,yj,tn))+f(xi,yj,tn−1,u(xi,yj,tn−1))].

Denote by the grid functions on with , , . Eq. (2.2) can be expressed as

 ∫21p(β)C0DβtUn−12ijdβ= ∂2∂x2Un−12ij+∂2∂y2Un−12ij (2.3) +12[f(xi,yj,tn,Unij)+f(xi,yj,tn−1,Un−1ij)]

Firstly we discretize the integral term in (2.3). Suppose , and . Let be a positive integer, and be the uniform step size. Take , , then the mid-point quadrature rule is used for approximating the integral in (2.3)

 ΔβK∑l=1p(βl)C0DβltUn−12ij+R1=∂2∂x2Un−12ij+∂2∂y2Un−12ij (2.4) +12[f(xi,yj,tn,U(xi,yj,tn))+f(xi,yj,tn−1,U(xi,yj,tn−1))],

where .

Next, we solve the multi-term time fractional wave equation (2.4) with the initial and boundary conditions (1.3) and (1.2). Suppose . According to Theorem 8.2.5 in Sun2009 (), the Caputo derivative , have the fully discrete difference scheme

 C0DβltUn−12ij (2.5) = τ1−βlΓ(3−βl)[a(βl)0δtUn−12ij−n−1∑k=1(a(βl)n−k−1−a(βl)n−k)δtUk−12ij−a(βl)n−1ψ2(xi,yj)]+Rl2,

where

 a(βl)k=(k+1)2−βl−k2−βl,k=0,1,2,⋯,

and

 ∣Rl2∣≤ 1Γ(3−βl)[2−βl12+23−βl3−βl−(1+21−βl)+112]⋅ (2.6) max0≤t≤tn∣∂3u(xi,yj,t)∂t3∣τ3−βl,l=1,2,⋯,K.

In the meantime, using the second order finite difference

 ∂2g(xi)∂x2=g(xi+1)−2g(xi)+g(xi−1)(Δx)2−(Δx)212∂4g(ξi)∂x4,ξi∈(xi−1,xi+1)

to approximate the second order derivatives in (2.4), it is obtained

 ΔβK∑l=1p(βl)τ1−βlΓ(3−βl)[a(βl)0δtUn−12ij−n−1∑k=1(a(βl)n−k−1−a(βl)n−k)δtUk−12ij (2.7) −a(βl)n−1ψ2(xi,yj)]+K∑l=1Δβp(βl)Rl2+R1 = δ2xUn−12ij+δ2yUn−12ij+12(f(xi,yj,tn−1,Un−1ij)+f(xi,yj,tn,Unij))+R3,

where . Subsequently, the nonlinear source term is dealt with in the following manner to avoid a system of nonlinear equations when computing:

 f(xi,yj,tn,Unij)=f(xi,yj,tn−1,Un−1ij)+O(τ). (2.8)

Substituting (2.8) in (2.7), we are left with

 ΔβK∑l=1p(βl)τ1−βlΓ(3−βl)[a(βl)0δtUn−12ij−n−1∑k=1(a(βl)n−k−1−a(βl)n−k)δtUk−12ij−a(βl)n−1ψ2(xi,yj)] (2.9) =δ2xUn−12ij+δ2yUn−12ij+f(xi,yj,tn−1,Un−1ij)+Rn−12ij+˜Rn−12ij,

where

 Rn−12ij=−K∑l=1Δβp(βl)Rl2+O(h21+h22)+O(Δβ2)

and

 ˜Rn−12ij=O(τ).

From (2.6), we can deduce that there exists a positive constant such that

 ∣∣∣−K∑l=1Δβp(βl)Rl2∣∣∣≤C1τ1+12ΔβK∑l=1Δβp(βl).

Since

 K∑l=1Δβp(βl)∼∫21p(β)dβ=c0,

we get

 K∑l=1Δβp(βl)≤C2,

where is a positive constant. Thus there exists a positive constant such that

 ∣∣Rn−12ij∣∣≤C3(τ1+12Δβ+h21+h22+Δβ2).

Besides, it is obvious that

 ∣∣˜Rn−12ij∣∣≤C4τ,

where is a positive constant.

Denote

 μ=ΔβK∑l=1p(βl)1τβlΓ(3−βl).

Since

 ΔβK∑l=1p(βl)1τβlΓ(3−βl) ∼ ∫21p(β)1τβΓ(3−β)dβ = p(β∗)Γ(3−β∗)∫211τβdβ = p(β∗)Γ(3−β∗)1−ττ2∣lnτ∣,

it can be concluded that

 μ=1O(τ2|lnτ|).

In addition, for any positive and small when is sufficiently small, thus the term is almost the same as when is sufficiently small. Adding the high order term

 τ4μδ2xδ2yUnij−Un−1ijτ

on both sides of (2.9), we derive

 ΔβK∑l=1p(βl)τ1−βlΓ(3−βl)[a(βl)0δtUn−12ij−n−1∑k=1(a(βl)n−k−1−a(βl)n−k)δtUk−12ij−a(βl)n−1ψ2(xi,yj)] (2.10) +τ4μδ2xδ2yUnij−Un−1ijτ = δ2xUn−12ij+δ2yUn−12ij+f(xi,yj,tn−1,Un−1ij)+Rn−12ij+˜Rn−12ij+ˆRn−12ij,

where

 ˆRn−12ij=τ4μδ2xδ2yUnij−Un−1ijτ,

and it is clear that

 ∣∣ˆRn−12ij∣∣≤C5τ3|lnτ|.

Also, for the initial and boundary value conditions, we have

 U0ij=ψ1(xi,yj), (xi,yj)∈Ω, (2.11)
 Unij=ϕ(xi,yj,tn), (xi,yj)∈∂Ω, 0≤n≤N. (2.12)

Let be the numerical approximation to . Neglecting the small term , and in (2.10), and using instead of in (2.10)(2.12), we construct the difference scheme for (1)(1.3) as follows:

 ΔβK∑l=1p(βl)τ1−βlΓ(3−βl)[a(βl)0δtun−12ij−n−1∑k=1(a(βl)n−k−1−a(βl)n−k)δtuk−12ij −a(βl)n−1(ψ2)ij]+τ4μδ2xδ2yunij−un−1ijτ = δ2xun−12ij+δ2yun−12ij+f(xi,yj,tn−1,un−1ij), 1≤i≤M1−1, 1≤j≤M2−1, 1≤n≤N, (2.13) u0ij=(ψ1)ij, 1≤i≤M1−1, 1≤j≤M2−1, (2.14) unij=ϕnij, (i,j)∈γ={(i,j) | (xi,yj)∈∂Ω}, 0≤n≤N, (2.15)

where

 (ψ1)ij=ψ1(xi,yj), (ψ2)ij=ψ2(xi,yj), 1≤i≤M1−1, 1≤j≤M2−1,

and

 ϕnij=ϕ(xi,yj,tn),(i,j)∈γ,0≤n≤N.

Notice , then Eq. (2) can be rewritten as:

 ΔβK∑l=1p(βl)1τβlΓ(3−βl)unij−12δ2xunij−12δ2yunij+14μδ2xδ2yunij = ΔβK∑l=1p(βl)1τβlΓ(3−βl)[un−1ij+n−1∑k=1(a(βl)n−k−1−a(βl)n−k)(ukij−uk−1ij) +τa(βl)n−1(ψ2)ij]+12δ2xun−1ij+12δ2yun−1ij+14μδ2xδ2yun−1ij+f(xi,yj,tn−1,un−1ij),

or

 = [n−1∑k=1(a(βl)n−k−1−a(βl)n−k)(ukij−uk−1ij)+τa(βl)n−1(ψ2)ij]+f(xi,yj,tn−1,un−1ij),

where denotes the identity operator.

Let

 u∗ij=(√μI−12√μδ2y)unij.

Together with (2.14) and (2.15) the ADI difference scheme is derived, and the procedure can be executed as follows:

On each time level , firstly, for all fixed , solving a set of equations at the mesh points to get the intermediate solution :

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(√μI−12√μδ2x)u∗ij=(√μI+12√μδ2x)(√μI+12√μδ2y)un−1ij+ΔβK∑l=1p(βl)1τβlΓ(3−βl)[n−1∑k=1(a(βl)n−k−1−a(βl)n−k)(ukij−uk−1ij)+τa(βl)n−1(ψ2)ij]+f(xi,yj,tn−1,un−1ij),1≤i≤M1−1,u∗0j=(√μI−12√μδ2y)un0j,u∗M1j=(√μI−12√μδ2y)unM1j; (2.16)

afterwards, for all fixed , by computing a set of equations at the mesh points , the solution can be obtained:

 ⎧⎪ ⎪⎨⎪ ⎪⎩(√μI−12√μδ2y)unij=u∗ij,1≤j≤M2−1,uni0=ϕ(xi,y0,tn),uniM2=ϕ(xi,yM2,tn). (2.17)

## 3 Analysis of the ADI difference scheme

### 3.1 Solvability

It is clear that the ADI scheme (2.16)(2.17) is a linear tridiagonal system in unknowns, and the coefficient matrices are strictly diagonally dominant. Thus the scheme (2.16)(2.17) has a unique solution. This result can be written as following.

###### Theorem 3.1.

The ADI difference scheme (2.16)(2.17) is uniquely solvable.

### 3.2 Stability

In this subsection we prove the unconditional stability and the convergence of the difference scheme (2.16)(2.17). We start with some auxiliary definitions and useful results.

Denote the space of grid functions on

 Vh={v∣v={vij∣(xi,yj)∈¯Ωh} and vij=0 if (xi,yj)∈∂Ωh}.

For any grid function , the following discrete norms and Sobolev seminorm are introduced:

 ∥v∥= ⎷h1h2M1−1∑i=1M2−1∑j=1|vij|2,∥δxδyv∥= ⎷h1h2M1∑i=1M2∑j=1|δxδyvi−12,j−12|2,
 ∥δxv∥= ⎷h1h2M1∑i=1M2−1∑j=1|δxvi−12,j|2,∥δyv∥= ⎷h1h2M1−1∑i=1M2∑j=1|δyvi,j−12|2,
 ∥Δhv∥= ⎷h1h2M1−1∑i=1M2−1∑j=1|Δhvij|2,|v|1=√∥δxv∥2+∥δ