# Numerical Solution of Time-Fractional Diffusion-Wave Equations via Chebyshev Wavelets Collocation Method

Authors: Fengying Zhou; Xiaoyong Xu
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

The second-kind Chebyshev wavelets collocation method is applied for solving a class of time-fractional diffusion-wave equation. Fractional integral formula of a single Chebyshev wavelet in the Riemann-Liouville sense is derived by means of shifted Chebyshev polynomials of the second kind. Moreover, convergence and accuracy estimation of the second-kind Chebyshev wavelets expansion of two dimensions are given. During the process of establishing the expression of the solution, all the initial and boundary conditions are taken into account automatically, which is very convenient for solving the problem under consideration. Based on the collocation technique, the second-kind Chebyshev wavelets are used to reduce the problem to the solution of a system of linear algebraic equations. Several examples are provided to confirm the reliability and effectiveness of the proposed method.


---

## Body

## 1. Introduction

Many phenomena in various fields of the science and engineering can be modeled by fractional differential equations, in which time-fractional diffusion-wave equation is a mathematical model of a wide class of important physical phenomena. It is obtained from the diffusion-wave equation by replacing the second-order time derivative term by a fractional derivative order1<α≤2. In this paper, our study focuses on the following time-fractional diffusion-wave ([1]) with the Caputo fractional derivative (1)∂αμx,t∂tα+λ∂μx,t∂t=∂2μx,t∂x2+qx,t,0≤x≤1,0≤t≤1, with the initial condition (2)μx,0=f0x,∂μx,0∂t=f1x,0≤x≤1, and the boundary conditions (3)μ0,t=g0t,μ1,t=g1t,0≤t≤1, where f0(x),f1(x),g0(t),g1(t) are given functions with second-order continuous derivatives and λ>0 is a constant and q(x,t) is a given known function.It is noted that most fractional diffusion-wave equations do not have closed form solutions. Many researchers have proposed various methods to solve the time-fractional diffusion-wave equations from the perspective of analytical solution and numerical solution. The method of separation of variables in [1], Sumudu transform method in [2], and decomposition method in [3] were used to construct analytical approximate solutions of fractional diffusion-wave equations, respectively. Finite difference schemes in [4–7] were widely used to solve the numerical solutions of the fractional diffusion-wave equations. The authors of [8] employed radial point interpolation method for solving the fractional diffusion-wave equations. B-spline collocation method was proposed to solve the fractional diffusion-wave equations in [9]. In [10, 11], Sinc-finite difference method and Sinc-Chebyshev method were employed for solving the fractional diffusion-wave equations respectively. Recently methods based on operational matrix of Jacobi and Chebyshev polynomials were proposed to deal with the fractional diffusion-wave equations ([12–14]). In [15], the authors applied fractional order Legendre functions method depending on the choices of two parameters to solve the fractional diffusion-wave equations. Two-dimensional Bernoulli wavelets with satisfier function in the Ritz-Galerkin method were proposed for the time-fractional diffusion-wave equation in [16]. The author of [17] proposed a numerical method based on the Legendre wavelets with their operational matrix of fractional integral to solve the time-fractional diffusion-wave equations.Since fractional derivative is a nonlocal operator, it is natural to consider a global scheme such as the collocation method for its numerical solution. Spectral methods are widely used in seeking numerical solutions of fractional order differential equations, due to their excellent error properties and exponential rates of convergence for smooth problems. Collocation methods, one of the three most common spectral schemes, have been applied successfully to numerical simulations of many problems in science and engineering.Wavelets, as another basis set and very well-localized functions, are considerably useful for solving differential and integral equations. Particularly, orthogonal wavelets are widely used in approximating numerical solutions of various types of fractional order differential equations in the relevant literatures; see [18–22]. Among them, the second-kind Chebyshev wavelets have gained much attention due to their useful properties ([23–26]) and can handle different types of differential problems. It is observed that most papers using these wavelets methods to approximate numerical solutions of fractional order differential equations are based on the operational matrix of fractional integral or fractional derivatives. It is inevitable to produce approximation error during the process of constructing the operational matrix. Regarding this point, analysis in [27] shows some disadvantages of using the operational matrix of Legendre and Chebyshev wavelets.Inspired and motivated by the work mentioned above, our main purpose of this paper is to extend the second-kind Chebyshev wavelets for solving time-fractional diffusion-wave equations (1)–(3) and to show that it is not necessary to establish the operational matrix of fractional integrals and fractional derivatives when applying wavelets to solve various types of fractional partial differential equations. To reduce the approximation error at most during the calculation process, fractional integral formula of a single Chebyshev wavelet in the Riemann-Liouville sense is derived by means of the shifted Chebyshev polynomials of the second kind. By utilizing the collocation method and some properties of the second-kind Chebyshev wavelets, the problem under consideration is reduced to the solution to a system of linear algebraic equations. The proposed method is very convenient for solving such problems, since the initial and boundary conditions are taken into account automatically.The rest of the paper is organized as follows. Section2 describes some necessary definitions and preliminaries of calculus. Section 3 gives the convergence and accuracy estimation of the second-kind Chebyshev wavelets expansion of two-dimension. Section 4 is devoted to deriving the fractional integral formula of a single Chebyshev wavelet in the Riemann-Liouville sense by means of shifted Chebyshev polynomials of the second kind. The proposed method is described for solving time-fractional diffusion-wave equations in Section 5. In Section 6, numerical results of some test problems are presented. Finally, a brief conclusion is given in Section 7.


## 2. Definitions and Preliminaries

In this section, we present some necessary definitions and preliminaries of the fractional calculus theory which will be used later. The widely used definitions of fractional integral and fractional derivative are the Riemann-Liouville definition and the Caputo definition, respectively.Definition 1.
A real functionf(x), x>0, is said to be in the space Cσ, σ∈R, if there is a real number ρ with ρ>σ such that f(x)=xρf0(x), where f0(x)∈C[0,∞), and f(x)∈Cσn if f(n)(x)∈Cσ, n∈N.Definition 2 (see [28]).
The Riemann-Liouville fractional integral operatorIα of order α(α≥0) for a function f(x)∈Cσ(σ≥-1) is defined as (4)Iαfx=fx,α=0,1Γα∫0xx-tα-1ftdt,α>0.Definition 3 (see [28]).
The Caputo fractional derivative operatorDα of order α(α≥0) for a function f(x)∈C1n is defined as (5)Dαfx=fnx,α=n∈N,1Γn-α∫0xfntx-tα+1-ndt,n-1<α<n.Some important properties of the operatorsIα and Dα are needed in this paper; we only mention the following properties:(1)
Iα1Iα2f(x)=Iα1+α2f(x) for α1,α2>0.(2)
DαIαf(x)=f(x), DβIαf(x)=Iα-βf(x), α>β.(3)
IαDαf(x)=f(x)-∑k=0⌈α⌉-1‍f(k)(0+)xk/k!, x>0, where the ceiling function ⌈α⌉ denotes the integer smaller than or equal to α and N0={0,1,2,…}. One can see more details about fractional calculus in [28].


## 3. The Second-Kind Chebyshev Wavelets and Their Properties

The second-kind Chebyshev waveletsψn,m(t)=ψ(k,n,m,t) have four arguments: k can assume any positive integer, n=1,2,3,…,2k-1, m, is the degree of the second-kind Chebyshev polynomials and t is the normalized time. They are defined on the interval [0,1) as (6)ψn,mt=2k/2U~m2kt-2n+1,n-12k-1≤t<n2k-1,0,otherwise,where(7)U~mt=2πUmt,m=0,1,2,…. Here Um(t) are the second-kind Chebyshev polynomials of degree m which are orthogonal with respect to the weight function ω(t)=1-t2 on the interval [-1,1] and satisfy the following recursive formula: (8)U0t=1,U1t=2t,Um+1t=2tUmt-Um-1t,m=1,2,3,…. Note that when dealing with the second-kind Chebyshev wavelets the weight function has to be dilated and translated as (9)ωnt=ω2kt-2n+1.A functionf(x)∈L2(R) defined over [0,1) may be expanded by the second-kind Chebyshev wavelets as (10)fx=∑n=1∞∑m=0∞cn,mψn,mx, where (11)cn,m=fx,ψn,mxLω20,1=∫01fxψn,mxωnxdx,in which 〈·,·〉Lω2[0,1) denotes the inner product in Lω2[0,1). If the infinite series in (10) is truncated, then it can be written as (12)fx≅∑n=12k-1∑m=0M-1cn,mψn,mx=CTΨx, where C and Ψ(x) are 2k-1M×1 matrices given by (13)C=c1,0c1,1⋯c1,M-1c2,0c2,1⋯c2,M-1⋯c2k-1,0c2k-1,1⋯c2k-1,M-1T,Ψx=ψ1,0ψ1,1⋯ψ1,M-1ψ2,0ψ2,1⋯ψ2,M-1⋯ψ2k-1,0ψ2k-1,1⋯ψ2k-1,M-1T.Observe that{ψn,m(x)ψn′,m′(y):n,n′=1,2,3,…,2k-1,m,m′=0,1,2,…} is an orthonormal set over [0,1)×[0,1). A function μ(x,y)∈L2(R2) defined over [0,1)×[0,1) may be expanded as (14)μx,y=∑n=12k-1∑n′=12k-1∑m=0∞∑m′=0∞dn,n′,m,m′ψn,mxψn′,m′y, where (15)dn,n′,m,m′=μx,y,ψn,mxψn′,m′yLω20,1×0,1=∫01∫01μx,yψn,mxψn′,m′yωnxωn′ydxdy.If the infinite series in (14) is truncated, then (14) can be written as(16)μx,y≅∑n=12k-1∑n′=12k-1∑m=0M-1∑m′=0M-1dn,n′,m,m′ψn,mxψn′,m′y=∑n=12k-1∑m=0M-1ψn,mx∑n′=12k-1∑m′=0M-1dn,n′,m,m′ψn′,m′y=∑n=12k-1∑m=0M-1ψn,mxDn,mTΨy,where Dn,m and Ψ(y) are 2k-1M×1 matrices given by (17)Dn,m=dn,1,m,0dn,1,m,1⋯dn,1,m,M-1⋯dn,2k-1,m,0dn,2k-1,m,1⋯dn,2k-1,m,M-1T,Ψy=ψ1,0yψ1,1y⋯ψ1,M-1y⋯ψ2k-1,0yψ2k-1,1y⋯ψ2k-1,M-1yT.Write Ψn(x)=ψn,0(x)ψn,1(x)⋯ψn,M-1(x) and Dn=Dn,0Dn,1⋯Dn,M-1T, n=1,2,3,…,2k-1; then μ(x,y) can be written as (18)μx,y≅∑n=12k-1∑m=0M-1ψn,mxDn,mTΨy=∑n=12k-1ΨnxDnΨy=Ψ1x⋯Ψ2k-1xD1D2⋮D2k-1Ψy=ΨTxDΨy, where D is a 2k-1M×2k-1M matrix given by D=D1D2⋮D2k-1.The following theorems give the convergence and accuracy estimation of the second-kind Chebyshev wavelets expansion.Theorem 4 (see [29]).
Letf(x) be a second-order derivative square-integrable function defined on [0,1) with bounded second-order derivative; say f′′x≤B for some constant B; then
(i)f(x) can be expanded as an infinite sum of the second-kind Chebyshev wavelets and the series converges to f(x) uniformly; that is, (19)fx=∑n=1∞∑m=0∞cn,mψn,mx, where cn,m=〈f(x),ψn,m(x)〉Lω2[0,1);
(ii)(20)σf,k,M<πB23∑n=2k-1+1∞1n5∑m=M∞1m-141/2,where σf,k,M=(∫01‍f(x)-∑n=12k-1‍∑m=0M-1‍cn,mψn,m(x)2ωn(x)dx)1/2.Theorem 5.
Supposing thatμ(x,y)∈L2(R2) is a continuous function defined on [0,1)×[0,1), ∂2μ(x,y)/∂x2, ∂2μ(x,y)/∂y2, and ∂4μ(x,y)/∂x2∂y2 are bounded with some positive constant B. Then, for any positive integer k,
(i) the series(21)∑n=12k-1∑n′=12k-1∑m=0∞∑m′=0∞dn,m,n′,m′ψn,mxψn′,m′yconverges to μ(x,y) uniformly in L2(R2); that is, (22)μx,y=∑n=12k-1∑n′=12k-1∑m=0∞∑m′=0∞dn,m,n′,m′ψn,mxψn′,m′y,where dn,m,n′,m′=〈μ(x,y),ψn,m(x)ψn′,m′(y)〉Lω2([0,1)×[0,1));
(ii)(23)σμ,k,M2<Bπ226∑n=12k-11n3∑n′=12k-11n′3∑m=M∞1m2m-12+Bπ2215∑n=12k-11n5∑n′=12k-11n′5∑m=2∞1m2m-12∑m′=M∞1m′2m′-12,where (24)σμ,k,M=∫01∫01μx,y-∑n=12k-1∑n′=12k-1∑m=0M-1∑m′=0M-1dn,n′,m,m′ψn,mxψn′,m′y2ωnxωn′ydxdy1/2.Proof.
(i)(25)dn,m,n′,m′=μx,y,ψn,mxψn′,m′yLω20,1×0,1=∫01∫01μx,yψn,mxψn′,m′yωnxωn′ydxdy=∫01ψn′,m′yωn′y∫01μx,yψn,mxωnxdxdy=∫01ψn′,m′yωn′yAn,mydy, where An,m(y)=∫01‍μ(x,y)ψn,m(x)ωn(x)dx. Similar to the proof of Theorem  1 in [29], we get, for m>1, (26)An,my=2-5k/222π1mm-1m+1∫0π∂2μcos⁡θ+2n-1/2k,y∂θ2τmθdθ, where τm(θ)=m+1sin⁡m-1θsin⁡θ-m-1sin⁡m+1θsin⁡θ. So, for m>1 and m′>1, (27)dn,m,n′,m′=∫01ψn′,m′yωn′yAn,mydy=2-5k/222π1mm-1m+1∫0πτmθ∫01∂2μcos⁡θ+2n-1/2k,y∂θ2·ψn′,m′yωn′ydydθ=2-5k/222π1mm-1m+1∫0πτmθ·2-5k/222π1m′m′-1m′+1·∫0π∂4μcos⁡θ+2n-1/2k,cos⁡η+2n′-1/2k∂θ2∂η2τm′ηdηdθ=2-5k23π1mm-1m+11m′m′-1m′+1·∫0π∫0π∂4μcos⁡θ+2n-1/2k,cos⁡η+2n′-1/2k∂θ2∂η2τmθτm′ηdθdη. Note that, for m>1, (28)∫0πτmθdθ2≤π∫0πτmθ2dθ≤π∫0πm+1sin⁡m-1θ-m-1sin⁡m+1θ2dθ≤π∫0πm+12sin2⁡m-1θdθ+∫0πm-12sin2⁡m+1θdθ=π22m+12+m-12≤π2m+12. Hence, for m>1 and m′>1, (29)dn,m,n′,m′≤Bπ28·1n5/2·1n′5/2·1mm-1·1m′m′-1 since n≤2k-1 and n′≤2k-1.
Form>1 and m′=0,1, (30)dn,m,n′,m′=∫01ψn′,m′yωn′yAn,mydy=2-5k/222π1mm-1m+1∫0πτmθ∫01∂2μcos⁡θ+2n-1/2k,y∂θ2·ψn′,m′yωn′ydydθ=2-5k/222π1mm-1m+1∫0πτmθ∫n′-1/2k-1n′/2k-1∂2μcos⁡θ+2n-1/2k,y∂θ2·2k/22πUm′2ky-2n′+1ω2ky-2n′+1dydθ=2-3k2π1mm-1m+1∫0πτmθ∫0π∂2μcos⁡θ+2n-1/2k,cos⁡η+2n′-1/2k∂θ2·sin⁡m′+1ηsin⁡ηdηdθ, where we use the variable substitution 2ky-2n′+1=cos⁡η and the definition of the second-kind Chebyshev polynomials. Therefore, for m>1 and m′=0,1, we obtain (31)dn,m,n′,m′≤Bπ24·1n3/2·1n′3/2·1mm-1 by n≤2k-1 and n′≤2k-1. Similarly, for m′>1 and m=0,1, we also have (32)dn,m,n′,m′≤Bπ24·1n3/2·1n′3/2·1m′m′-1. Relations (29), (31), and (32) show that the series (33)∑n=12k-1∑n′=12k-1∑m=0∞∑m′=0∞dn,m,n′,m′ψn,mxψn′,m′yconverges to μ(x,y) uniformly in L2(R2).
(ii)(34)σμ,k,M2=∫01∫01μx,y-∑n=12k-1∑n′=12k-1∑m=0M-1∑m′=0M-1dn,n′,m,m′ψn,mxψn′,m′y2·ωnxωn′ydxdy=∫01∫01∑n=12k-1∑n′=12k-1∑m′=0M-1∑m=M∞​+∑m=0M-1∑m′=M∞​+∑m′=M∞∑m=M∞​dn,n′,m,m′·ψn,mxψn′,m′y2ωnxωn′ydxdy. According to the orthonormality of ψn,m(x)ψn′,m′(y), we have (35)σμ,k,M2=∑n=12k-1∑n′=12k-1∑m′=0M-1∑m=M∞​+∑m=0M-1∑m′=M∞​+∑m′=M∞∑m=M∞​dn,n′,m,m′2=∑n=12k-1∑n′=12k-1∑m′=01∑m=M∞​+∑m′=2M-1∑m=M∞​+∑m=01∑m′=M∞​+∑m=2M-1∑m′=M∞​+∑m′=M∞∑m=M∞​·dn,n′,m,m′2. By relations (29), (31), and (32), it gives (36)σμ,k,M2≤∑n=12k-1∑n′=12k-1∑m′=01∑m=M∞Bπ2281n31n′31m2m-12+∑m′=2M-1∑m=M∞Bπ22161n51n′51m2m-121m′2m′-12+∑m=01∑m′=M∞Bπ2281n31n′31m′2m′-12+∑m=2M-1∑m′=M∞Bπ22161n51n′51m2m-121m′2m′-12+∑m=M∞∑m′=M∞Bπ22161n51n′51m2m-121m′2m′-12=∑n=12k-1∑n′=12k-1∑m=M∞Bπ2261n31n′31m2m-12+∑m=2M-1∑m′=M∞Bπ22151n51n′51m2m-121m′2m′-12+∑m=M∞∑m′=M∞Bπ22161n51n′51m2m-121m′2m′-12<Bπ226∑n=12k-11n3∑n′=12k-11n′3∑m=M∞1m2m-12+Bπ2215∑n=12k-11n5∑n′=12k-11n′5∑m=2∞1m2m-12∑m′=M∞1m′2m′-12.The proof is completed.


## 4. The Fractional Integral of a Single Second-Kind Chebyshev Wavelet

In this section, fractional integral formula of a single Chebyshev wavelet in the Riemann-Liouville sense is derived by means of the shifted second-kind Chebyshev polynomialsUm∗, which plays an important role in dealing with the time-fractional diffusion-wave equations.Theorem 6.
The fractional integral of a Chebyshev wavelet defined on the interval[0,1] with compact support [n-1/2k-1,n/2k-1] is given by (37)Iαψn,mx=0,x<n-12k-1,1Γα2k/22π∑r=0m∑i=rm∑j=0rTi,i-rm,n,k-1jj+αCrjxr-jx-n-12k-1j+α,n-12k-1≤x≤n2k-1,1Γα2k/22π∑r=0m∑i=rm∑j=0rTi,i-rm,n,k-1jj+αCrjxr-jx-n-12k-1j+α-x-n2k-1j+α,x>n2k-1,where Ti,i-rm,n,k=(-1)m-r22i2r(k-1)(n-1)i-rΓ(m+i+2)/Γ(m-i+1)Γ(2i+2)i!/(i-r)!r!,Crj=r!/j!(r-j)!.Proof.
The analytical form of the shifted Chebyshev polynomials [30] of the second-kind Um∗ of degree m is given by(38)Um∗x=∑i=0m-1i22m-2iΓ2m-i+2Γi+1Γ2m-2i+2xm-i=∑i=0m-1m-i22iΓm+i+2Γm-i+1Γ2i+2xi.
According to the relation betweenUm∗(x) and Um(x), it gives that(39)Um2kx-2n+1=∑i=0m-1m-i22iΓm+i+2Γm-i+1Γ2i+22k-1x-n-1i=∑i=0m-1m-i22iΓm+i+2Γm-i+1Γ2i+22ik-1x-n-12k-1i=∑i=0m-1m-i22iΓm+i+2Γm-i+1Γ2i+22ik-1∑r=0iCirxi-r-n-12k-1r=∑i=0m∑r=0i-1m-i+r22i2i-rk-1n-1rΓm+i+2Γm-i+1Γ2i+2i!r!i-r!xi-r.By interchanging the summation and substituting i-r with r, Um(2kx-2n+1) can be written as(40)Um2kx-2n+1=∑r=0m∑i=rm-1m-r22i2rk-1n-1i-rΓm+i+2Γm-i+1Γ2i+2i!i-r!r!xr.
LetTi,i-rm,n,k=(-1)m-r22i2r(k-1)(n-1)i-rΓ(m+i+2)/Γ(m-i+1)Γ(2i+2)i!/(i-r)!r!; then(41)Um2kx-2n+1=∑r=0m∑i=rmTi,i-rm,n,kxr. Therefore (42)ψn,mx=2k/22π∑r=0m∑i=rmTi,i-rm,n,kxr,n-12k-1≤x<n2k-1.0,otherwise. So, when n-1/2k-1≤x≤n/2k-1, let u=x-t; then (43)Iαxr=1Γα∫n-1/2k-1xx-tα-1trdt=1Γα∫0x-n-1/2k-1uα-1x-urdu=1Γα∫0x-n-1/2k-1uα-1∑j=0rCrjxr-j-ujdu=1Γα∑j=0r-1jj+αCrjxr-jx-n-12k-1j+α. Similarly, when x>n/2k-1, we have (44)Iαxr=1Γα∫n-1/2k-1n/2k-1x-tα-1trdt=1Γα∫x-n/2k-1x-n-1/2k-1uα-1x-urdu=1Γα∫x-n/2k-1x-n-1/2k-1uα-1∑j=0rCrjxr-j-ujdu=1Γα∑j=0r-1jj+αCrjxr-jx-n-12k-1j+α-x-n2k-1j+α.Applying the Riemann-Liouville fractional integral of order α with respect to x on ψn,m(x), we obtain (45)Iαψn,mx=0,x<n-12k-1,1Γα∫n-1/2k-1xx-tα-1ψn,mtdt,n-12k-1≤x≤n2k-1,1Γα∫n-1/2k-1n/2k-1x-tα-1ψn,mtdt,x>n2k-1,=0,x<n-12k-1,1Γα2k/22π∑r=0m∑i=rmTi,i-rm,n,k∫n-1/2k-1xx-tα-1trdt,n-12k-1≤x≤n2k-1,1Γα2k/22π∑r=0m∑i=rmTi,i-rm,n,k∫n-1/2k-1n/2k-1x-tα-1trdt,x>n2k-1,=0,x<n-12k-1,1Γα2k/22π∑r=0m∑i=rmTi,i-rm,n,k∑j=0r-1jj+αCrjxr-jx-n-12k-1j+α,n-12k-1≤x≤n2k-1,1Γα2k/22π∑r=0m∑i=rmTi,i-rm,n,k∑j=0r-1jj+αCrjxr-jx-n-12k-1j+α-x-n2k-1j+α,x>n2k-1. Thus, we have (46)Iαψn,mx=0,x<n-12k-1,1Γα2k/22π∑r=0m∑i=rm∑j=0rTi,i-rm,n,k-1jj+αCrjx-n-12k-1j+α,n-12k-1≤x≤n2k-1,1Γα2k/22π∑r=0m∑i=rm∑j=0rTi,i-rm,n,k-1jj+αCrjxr-jx-n-12k-1j+α-x-n2k-1j+α,x>n2k-1.The proof is completed.For example, in the case ofk=2,M=3,x=0.65,α=1.7, we obtain (47)IαΨ6×10.65=0.455620526461909-0.137175775876020.1444590990338910.041066248466599-0.06388083094804290.0620309999360338,where(48)Ψ6×1x=ψ1,0xψ1,1xψ1,2xψ2,0xψ2,1xψ2,2xT.


## 5. Description of the Proposed Method

Consider the time-fractional diffusion-wave equation with the following form:(49)∂αμx,t∂tα+∂μx,t∂t=∂μ2x,t∂x2+qx,t,0≤x≤1,0≤t≤1,1<α≤2,with initial condition (50)μx,0=f0x,μtx,0=f1x0≤x≤1, and boundary conditions (51)μ0,t=g0t,μ1,t=g1t,0≤t≤1, where f0(·),f1(·),g0(·),g1(·) are given functions with second-order continuous derivatives in L2[0,1) and q(·,·) is a given function in L2([0,1)×[0,1)).To solve this problem, we assume(52)∂4μx,t∂x2∂t2≈ΨTx·U·Ψt, where U=ui,j2k-1M×2k-1M is an unknown matrix which should be determined and Ψ(·) is as in (13). By integrating two times with respect to t on both sides of (52) and together with (2), we have (53)∂2μx,t∂x2≈f0′′x+tf1′′x+ΨTx·U·I2Ψt, Also by integrating (53) two times with respect to x, we obtain(54)∂μx,t∂x=∂μx,t∂xx=0+f0′x-f0′0+tf1′x-f1′0+IΨxT·U·I2Ψt,(55)μx,t≈μ0,t+x∂μx,t∂xx=0+f0x-f00-xf0′0+tf1x-f10-xf1′0+I2ΨxT·U·I2Ψt.Putting x=1 in (55) and considering the boundary conditions (51), we get(56)μ1,t≈μ0,t+∂μx,t∂xx=0+f01-f00-f0′0+tf11-f10-f1′0+I2Ψ1T·U·I2Ψt.Thus, we have(57)∂μx,t∂xx=0≈g1t-g0t-f01-f00-f0′0-tf11-f10-f1′0-I2Ψ1T·U·I2Ψt.Write(58)Ht=g1t-g0t-f01-f00-f0′0-tf11-f10-f1′0-I2Ψ1T·U·I2Ψt.So(59)μx,t≈g0t+xHt+f0x-f00-xf0′0+tf1x-f10-xf1′0+I2ΨxT·U·I2Ψt.Applying fractional differentiation of order α and order one on both sides of (59) with respect to t, we get (60)∂μαx,t∂tα≈Dαg0t+xDαHt+I2ΨxT·U·I2-αΨt,∂μx,t∂t≈g0′t+xH′t+f1x-f10-xf1′0+I2ΨxT·U·IΨt,where(61)DαHt=Dαg1t-Dαg0t-I2Ψ1T·U·I2-αΨt,H′t=g1′t-g0′t-f11-f10-f1′0-I2Ψ1T·U·IΨt.Now by substituting (53), (59), and (60) into (49), replacing ≈ by =, and taking collocation points xi=2i-1/2kM, tj=2j-1/2kM, i,j=1,2,…,2k-1M, we obtain the following linear system of algebraic equations: (62)Dαg0tj+xiDαHtj+I2ΨxiT·U·I2-αΨtj+g0′tj+xiH′tj+f1xi-f10-xif1′0+I2ΨxiT·U·Ψtj-f0′′xi+tjf1′′xi+ΨTxi·U·I2Ψtj=qxi,tj,i,j=1,2,…,2k-1M. By solving this system to determine the unknown matrix U, we can achieve an approximate solution for the problem by substituting U into (59).


## 6. Numerical Examples

In this section, we give some numerical examples to demonstrate the efficiency and reliability of the proposed method. In all the examples, the package of Matlab 2016a has been used.Example 1.
Consider the following time-fractional diffusion-wave equation ([17]) of order α(1<α≤2)(63)∂αμx,t∂tα+∂μx,t∂t=∂2μx,t∂x2+qx,t,0<x<1,0<t≤1, with the initial conditions (64)μx,0=0,μtx,0=0,0≤x≤1,and the boundary conditions (65)μ0,t=0,μ1,t=0,0<t≤1,where q(x,t)=(2t2-α/Γ(3-α)+2t)(x-x2)+2t2. The exact solution of this problem is μ(x,t)=t2x(1-x).
The problem is solved by the proposed method fork=2,M=3 and the CPU time is 3.818 seconds. Figure 1 shows the approximate solution and the absolute error of this problem for α=1.1. Table 1 gives the absolute errors for different values of α at different points. To make a comparison, in Table 2, we list the absolute error obtained by the Legendre wavelets method in [17] based on fractional operational matrix of integral with k=3,M=3. From Figure 1 and Table 1, it can be seen that the presented method is very efficient and accurate in solving this problem.Table 1
The absolute error for Example1 for some different values 1<α≤2 at some different points.
(
x
,
t
)
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
(
0.1,0.1
)
3.25261
e
-
19
4.33681
e
-
19
3.25261
e
-
19
8.67362
e
-
19
1.19262
e
-
18
(
0.2,0.2
)
5.20417
e
-
18
3.46945
e
-
18
6.07153
e
-
18
2.60209
e
-
18
7.80626
e
-
18
(
0.3,0.3
)
1.38778
e
-
17
6.93889
e
-
18
2.77556
e
-
17
6.93889
e
-
18
1.38778
e
-
17
(
0.4,0.4
)
2.77556
e
-
17
1.38778
e
-
17
5.55112
e
-
17
3.46945
e
-
17
2.77556
e
-
17
(
0.5,0.5
)
2.77556
e
-
17
0
8.32667
e
-
17
3.46945
e
-
17
4.16334
e
-
17
(
0.6,0.6
)
6.93889
e
-
17
4.16334
e
-
17
2.35922
e
-
16
1.38778
e
-
17
6.93889
e
-
17
(
0.7,0.7
)
1.66533
e
-
16
1.38778
e
-
16
3.60822
e
-
16
5.55112
e
-
17
2.77556
e
-
17
(
0.8,0.8
)
1.38778
e
-
16
1.94289
e
-
16
4.16334
e
-
16
2.77556
e
-
17
1.38778
e
-
16
(
0.9,0.9
)
1.38778
e
-
17
1.24900
e
-
16
3.19189
e
-
16
1.24900
e
-
16
1.24900
e
-
16Table 2
The absolute error for Example1 for some different values 1<α≤2 at some different points in [17].
(
x
,
t
)
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
(
0.1,0.1
)
4.3939
e
-
6
1.3694
e
-
5
2.2112
e
-
5
2.3491
e
-
5
1.2825
e
-
5
(
0.2,0.2
)
7.5711
e
-
6
2.2302
e
-
5
3.7623
e
-
5
4.6694
e
-
5
3.1949
e
-
5
(
0.3,0.3
)
8.3240
e
-
6
2.6323
e
-
5
4.9656
e
-
5
6.9397
e
-
5
5.3169
e
-
5
(
0.4,0.4
)
7.5553
e
-
6
2.4508
e
-
5
4.9748
e
-
5
7.7937
e
-
5
6.7115
e
-
5
(
0.5,0.5
)
5.9012
e
-
6
1.8821
e
-
5
3.9629
e
-
5
6.8327
e
-
5
6.7208
e
-
5
(
0.6,0.6
)
4.4231
e
-
6
1.3440
e
-
5
2.7980
e
-
5
5.1193
e
-
5
5.5604
e
-
5
(
0.7,0.7
)
2.9153
e
-
6
8.1172
e
-
6
1.5667
e
-
5
2.9601
e
-
5
3.6814
e
-
5
(
0.8,0.8
)
1.6695
e
-
6
4.1365
e
-
6
5.1427
e
-
6
8.8490
e
-
6
1.7911
e
-
5
(
0.9,0.9
)
6.9982
e
-
7
3.1435
e
-
6
2.0508
e
-
4
2.8695
e
-
4
1.2030
e
-
5Figure 1
Approximate solution (a) and absolute error (b) for Example1 with α=1.1, k=2, and M=3.
(a)
(b)Example 2.
Consider the following time-fractional diffusion-wave equation ([17]) of order α(1<α≤2)(66)∂αμx,t∂tα+∂μx,t∂t=∂2μx,t∂x2+qx,t,0<x<1,0<t≤1, with the initial conditions (67)μx,0=0,μtx,0=0,0≤x≤1,and the boundary conditions (68)μ0,t=t3,μ1,t=et3,0<t≤1,where q(x,t)=(6/Γ(4-α)t3-α+3t2-t3)ex. The exact solution of this problem is μ(x,t)=t3ex.
The problem is solved by the proposed method fork=2,M=6 and the CPU time is 28.79 seconds. Figure 2 shows the approximate solution and the absolute error of this problem for α=1.3. Table 3 gives the absolute error for different values of α at different points with k=2,M=6. In order to compare our results with obtained results in [17], we list the absolute error obtained by the Legendre wavelets method in [17] with k=2,M=5 in Table 4. Table 5 lists the maximum absolute errors obtained by the proposed method for different choices of M and α at the points (xi,tj), where xi=i/40,tj=j/40,i,j=0,1,2,…,40. It is obvious that the results provided by the proposed method are more accurate than those given in [17].Table 3
The absolute errors for Example2 for some different values 1<α≤2 at some different points.
(
x
,
t
)
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
(
0.1,0.1
)
3.94481
e
-
13
3.86294
e
-
13
3.59444
e
-
13
3.55895
e
-
13
3.40724
e
-
13
(
0.2,0.2
)
2.10685
e
-
12
1.30045
e
-
12
8.35330
e
-
13
7.96873
e
-
13
9.32525
e
-
13
(
0.3,0.3
)
1.80894
e
-
11
1.19129
e
-
11
6.13486
e
-
12
1.57591
e
-
12
1.36735
e
-
12
(
0.4,0.4
)
1.32073
e
-
10
1.12757
e
-
10
9.26274
e
-
11
7.36479
e
-
11
5.68607
e
-
11
(
0.5,0.5
)
3.59640
e
-
10
3.20785
e
-
10
2.78375
e
-
10
2.35879
e
-
10
1.96482
e
-
10
(
0.6,0.6
)
2.96547
e
-
10
2.39021
e
-
10
1.74491
e
-
10
1.07139
e
-
10
4.29922
e
-
11
(
0.7,0.7
)
5.52939
e
-
10
4.83775
e
-
10
4.04305
e
-
10
3.17635
e
-
10
2.30672
e
-
10
(
0.8,0.8
)
5.01937
e
-
10
4.34392
e
-
10
3.55367
e
-
10
2.65621
e
-
10
1.69385
e
-
10
(
0.9,0.9
)
9.61414
e
-
10
9.16178
e
-
10
8.62620
e
-
10
7.99952
e
-
10
7.29467
e
-
10Table 4
The absolute errors for Example2 for some different values 1<α≤2 at some different points in [17].
(
x
,
t
)
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
(
0.1,0.1
)
6.7028
e
-
5
6.2270
e
-
5
6.0407
e
-
5
4.9516
e
-
5
2.0243
e
-
5
(
0.2,0.2
)
1.8718
e
-
4
1.6817
e
-
5
1.5683
e
-
4
1.3074
e
-
4
5.9155
e
-
5
(
0.3,0.3
)
3.0913
e
-
4
2.8256
e
-
4
2.7985
e
-
4
2.3269
e
-
4
1.0947
e
-
4
(
0.4,0.4
)
4.0221
e
-
4
3.6211
e
-
4
3.7045
e
-
4
3.4739
e
-
4
1.6790
e
-
4
(
0.5,0.5
)
4.5801
e
-
4
3.8782
e
-
4
3.8089
e
-
4
3.2553
e
-
4
1.6277
e
-
4
(
0.6,0.6
)
4.5260
e
-
4
3.7198
e
-
4
3.6309
e
-
4
3.1089
e
-
4
1.9284
e
-
4
(
0.7,0.7
)
4.0597
e
-
4
3.0859
e
-
4
3.2603
e
-
4
4.6656
e
-
5
6.2825
e
-
5
(
0.8,0.8
)
3.1039
e
-
4
7.9174
e
-
4
6.5594
e
-
4
1.2388
e
-
4
1.0181
e
-
5
(
0.9,0.9
)
1.7283
e
-
4
2.1787
e
-
4
7.1269
e
-
3
1.6600
e
-
3
2.0918
e
-
5Table 5
Maximum absolute error for Example2 with various choices of M and α.
M
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
3
4.06786
e
-
6
3.90708
e
-
6
3.71627
e
-
6
3.48835
e
-
6
3.22168
e
-
6
4
6.03876
e
-
7
5.76233
e
-
7
5.43586
e
-
7
5.04721
e
-
7
4.59432
e
-
7
5
4.17452
e
-
9
3.98846
e
-
9
3.76854
e
-
9
3.50696
e
-
9
3.20222
e
-
9
6
4.74514
e
-
10
4.51960
e
-
10
4.25304
e
-
10
3.93658
e
-
10
3.56805
e
-
10
7
2.38765
e
-
12
2.22533
e
-
12
2.14806
e
-
12
1.98641
e
-
12
1.77769
e
-
12
8
2.42473
e
-
13
2.90434
e
-
13
3.26850
e
-
13
1.74305
e
-
13
1.15907
e
-
13Figure 2
Approximate solution (a) and absolute error (b) for Example2 with α=1.3, k=2, and M=6.
(a)
(b)Example 3.
Consider the following time-fractional diffusion-wave equation of orderα(1<α≤2): (69)∂αμx,t∂tα+∂μx,t∂t=∂2μx,t∂x2+qx,t,0<x<1,0<t≤1, with the initial conditions (70)μx,0=0,μtx,0=0,0≤x≤1,and the boundary conditions (71)μ0,t=0,μ1,t=t3sin2⁡1,0<t≤1,where q(x,t)=Γ(4)/Γ(4-α)t3-αsin2⁡x+3t2sin2⁡x-2t3cos⁡(2x). The exact solution of this problem is μ(x,t)=t3sin2⁡x.
The problem is solved by the proposed method fork=2, M=6 and the CPU time is 28.42 seconds. Figure 3 shows the approximate solution and the absolute error of this problem in the case of α=1.7,k=2, and M=6. Table 6 gives the absolute error of the approximate solutions for different values of α at different points with k=2,M=6. Table 7 lists the maximum absolute error obtained by the proposed method for different choices of M and α at the points (xi,tj), where xi=i/40,tj=j/40,i,j=0,1,2,…,40.Table 6
The absolute error for Example3 for some different values 1<α≤2 at some different points.
(
x
,
t
)
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
(
0.1,0.1
)
5.49562
e
-
13
5.47621
e
-
13
4.98133
e
-
13
4.68255
e
-
13
5.26687
e
-
13
(
0.2,0.2
)
5.77200
e
-
12
2.35711
e
-
12
4.05830
e
-
13
8.65065
e
-
14
1.51648
e
-
14
(
0.3,0.3
)
8.08760
e
-
11
5.48063
e
-
11
3.10533
e
-
11
1.33666
e
-
11
3.69026
e
-
12
(
0.4,0.4
)
4.20186
e
-
10
3.35630
e
-
10
2.47962
e
-
10
1.65984
e
-
10
9.57750
e
-
11
(
0.5,0.5
)
1.64121
e
-
9
1.46671
e
-
9
1.27599
e
-
9
1.08434
e
-
9
9.04889
e
-
10
(
0.6,0.6
)
2.05951
e
-
9
1.79826
e
-
9
1.50447
e
-
9
1.19672
e
-
9
9.03980
e
-
10
(
0.7,0.7
)
2.38806
e
-
9
2.07520
e
-
9
1.71522
e
-
9
1.32122
e
-
9
9.24245
e
-
10
(
0.8,0.8
)
2.39668
e
-
9
2.09691
e
-
9
1.74583
e
-
9
1.34653
e
-
9
9.14931
e
-
10
(
0.9,0.9
)
1.75917
e
-
9
1.56064
e
-
9
1.32547
e
-
9
1.04978
e
-
9
7.37741
e
-
10Table 7
Maximum absolute error for Example3 with various choices of M and α.
M
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
3
2.40504
e
-
5
2.26983
e
-
5
2.11041
e
-
5
1.92403
e
-
5
1.71638
e
-
5
4
3.99132
e
-
6
3.77465
e
-
6
3.51924
e
-
6
3.22175
e
-
6
2.89080
e
-
6
5
1.08639
e
-
7
1.02446
e
-
7
9.51434
e
-
8
8.66286
e
-
8
7.71522
e
-
8
6
1.23281
e
-
8
1.16368
e
-
8
1.08219
e
-
8
9.87176
e
-
9
8.81455
e
-
9
7
2.64246
e
-
10
2.49148
e
-
10
2.31377
e
-
10
2.10629
e
-
10
1.87484
e
-
10
8
3.21829
e
-
11
3.06864
e
-
11
2.88753
e
-
11
2.66064
e
-
11
2.40237
e
-
11Figure 3
Approximate solution (a) and absolute error (b) for Example3 with α=1.7, k=2, and M=6.
(a)
(b)Example 4.
Consider the following time-fractional diffusion-wave equation of orderα(1<α≤2)(72)∂αμx,t∂tα+∂μx,t∂t=∂2μx,t∂x2+fx,t,0<x<1,0<t≤1, with the initial conditions (73)μx,0=0,μtx,0=ex,0≤x≤1,and the boundary conditions (74)μ0,t=tα+3+t,μ1,t=tα+3+te,0<t≤1,where q(x,t)=(Γ(4+α)/Γ(4)t3+(α+3)tα+2+1-tα+3-t)ex. The exact solution of this problem is μ(x,t)=(tα+3+t)ex.
The space-time graph of the approximate solution and the absolute error forα=1.9,k=2, and M=6 is presented in Figure 4. The absolute error for different values of α at different points with k=2 and M=6 is tabulated in Table 8. Table 9 lists the maximum absolute error obtained by the proposed method for different choices of M and α at the points (xi,tj), where xi=i/40,tj=j/40,i,j=0,1,2,…,40. The CPU time for this problem is 28.32 seconds.Table 8
The absolute errors for Example4 for some different values 1<α≤2 at some different points.
(
x
,
t
)
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
(
0.1,0.1
)
1.06076
e
-
9
2.86893
e
-
9
1.38103
e
-
8
2.42944
e
-
8
1.75130
e
-
8
(
0.2,0.2
)
1.40209
e
-
9
6.00404
e
-
9
2.27258
e
-
8
5.24253
e
-
8
4.98649
e
-
8
(
0.3,0.3
)
3.37278
e
-
9
1.98379
e
-
9
2.09526
e
-
8
7.25016
e
-
8
8.37632
e
-
8
(
0.4,0.4
)
2.36158
e
-
9
9.26811
e
-
9
2.79830
e
-
8
8.15413
e
-
8
1.06968
e
-
7
(
0.5,0.5
)
8.93198
e
-
8
9.49760
e
-
8
4.89761
e
-
8
6.65251
e
-
8
1.02985
e
-
7
(
0.6,0.6
)
1.44456
e
-
7
2.91469
e
-
7
2.53681
e
-
7
1.67284
e
-
7
1.10555
e
-
7
(
0.7,0.7
)
8.49482
e
-
8
1.98030
e
-
7
2.18826
e
-
7
1.66823
e
-
7
9.46295
e
-
8
(
0.8,0.8
)
3.27487
e
-
8
9.88976
e
-
8
1.30364
e
-
7
1.11109
e
-
7
6.11390
e
-
8
(
0.9,0.9
)
2.73139
e
-
9
3.18928
e
-
8
5.20473
e
-
8
4.77641
e
-
8
2.47031
e
-
8Table 9
Maximum absolute error for Example4 with various choices of M and α.
M
α
=
1.1
α
=
1.3
α
=
1.5
α
=
1.7
α
=
1.9
3
2.27899
e
-
5
6.43156
e
-
5
9.16715
e
-
5
1.32018
e
-
4
1.02247
e
-
4
4
2.34479
e
-
6
5.42695
e
-
6
5.49280
e
-
6
4.91082
e
-
6
3.00590
e
-
6
5
4.92849
e
-
7
1.03746
e
-
6
9.65101
e
-
7
7.05607
e
-
7
3.62955
e
-
7
6
1.49081
e
-
7
2.97338
e
-
7
2.65557
e
-
7
1.97476
e
-
7
1.16463
e
-
7
7
5.83218
e
-
8
1.08640
e
-
7
8.96853
e
-
8
5.56267
e
-
8
4.37784
e
-
8
8
2.69626
e
-
8
4.77286
e
-
8
3.86373
e
-
8
2.74826
e
-
8
2.01824
e
-
8Figure 4
Approximate solution (a) and absolute error (b) for Example4 with α=1.9, k=2, and M=6.
(a)
(b)Example 5.
Consider the following time-fractional diffusion-wave equation:(75)∂1.5μx,t∂t1.5+∂μx,t∂t=∂2μx,t∂x2+qx,t,0<x<1,0<t≤1, with the initial conditions (76)μx,0=ex2,μtx,0=ex2,0≤x≤1,and the boundary conditions (77)μ0,t=et,μ1,t=et+1,0<t≤1,where q(x,t)=(erf(t)-4x2-1)ex2+t. The exact solution of this problem is μ(x,t)=ex2+t.
This problem is solved by the proposed method withk=2, M=3 to 8. The absolute error of the approximate solutions for some different points are shown in Table 10. Figure 5 shows the approximate solution and the absolute error of this problem in the case of k=2 and M=6. The CPU time for this problem is 28.82 seconds.Table 10
The absolute errors for Example5 at some different points with k=2 and M=3 to 8.
(
x
,
t
)
M
=
3
M
=
4
M
=
5
M
=
6
M
=
7
M
=
8
(
0.1,0.1
)
6.48090
e
-
8
9.39323
e
-
9
6.41101
e
-
10
6.04636
e
-
11
2.16316
e
-
12
1.01696
e
-
13
(
0.2,0.2
)
2.27341
e
-
8
5.69677
e
-
8
7.56717
e
-
9
1.05616
e
-
9
9.40217
e
-
11
1.22968
e
-
11
(
0.3,0.3
)
7.34805
e
-
6
1.29090
e
-
6
1.27698
e
-
7
2.06609
e
-
8
1.84093
e
-
9
2.37917
e
-
10
(
0.4,0.4
)
3.45200
e
-
5
6.82472
e
-
6
7.01216
e
-
7
1.09420
e
-
7
9.84536
e
-
9
1.27289
e
-
9
(
0.5,0.5
)
1.03114
e
-
4
2.27987
e
-
5
2.20253
e
-
6
3.58064
e
-
7
3.15850
e
-
8
4.12741
e
-
9
(
0.6,0.6
)
8.61804
e
-
5
1.96315
e
-
5
2.03271
e
-
6
3.53811
e
-
7
3.28933
e
-
8
4.39198
e
-
9
(
0.7,0.7
)
6.20577
e
-
5
2.06432
e
-
5
2.61296
e
-
6
3.78036
e
-
7
3.54969
e
-
8
4.53937
e
-
9
(
0.8,0.8
)
2.11165
e
-
4
1.58365
e
-
5
2.34163
e
-
6
3.37270
e
-
7
3.52359
e
-
8
4.14875
e
-
9
(
0.9,0.9
)
1.36820
e
-
4
1.04491
e
-
6
2.87015
e
-
6
1.50481
e
-
7
2.93362
e
-
8
2.69187
e
-
9Figure 5
Approximate solution (a) and absolute error (b) for Example5 with k=2 and M=6.
(a)
(b)From Figures1–5 and Tables 1–10, it can be seen that the proposed method is very efficient and accurate in solving this problem and the obtained approximate solutions are very close to the exact ones for all chosen α with 1<α≤2. It is worth noting that more accurate result can be obtained by increasing basis functions.Remark 6.
It should be noted that the presented method can also be applied to the calculation of long timeT. We only need to extend the second-kind Chebyshev wavelet basis functions on the interval [0,1] to the general interval [a,b]. We choose Example 2 to illustrate the feasibility of this method. For example, Figure 6 shows the behavior of exact and approximation solution and corresponding error at T=10.Figure 6
The exact and approximation solution (a) and corresponding error (b) atk=2, M=6, and T=10.
(a)
(b)


## 7. Conclusion

In this paper, the collocation method based on the second-kind Chebyshev wavelets has been applied to the time-fractional diffusion-wave equations. We derived the fractional integral formula of a single Chebyshev wavelet in the Riemann-Liouville sense via the shifted Chebyshev polynomials of the second kind. Convergence and accuracy estimation of the second-kind Chebyshev wavelets expansion of two-dimension were given. The second-kind Chebyshev wavelets and their properties have been used to convert the problems under consideration into some corresponding linear systems of algebraic equations. The proposed method is very convenient for solving time-fractional diffusion-wave equations, since the boundary conditions are taken into account automatically during the process of establishing the expression of the approximate solution. Applicability and accuracy have been tested on some numerical examples. We also compared our results with the results obtained by Legendre wavelet, which showed that the proposed approach was more accurate. Furthermore, the proposed method can be expected to solve other types of fractional partial differential equations numerically.



---

*Source: 2610804-2017-09-07.xml*