# Numerical Method to Modify the Fractional-Order Diffusion Equation

Authors: Yunkun Chen
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

Time or space or time-space fractional-order diffusion equations (FODEs) are widely used to describe anomalous diffusion processes in many physical and biological systems. In recent years, many authors have proposed different numerical methods to solve the modified fractional-order diffusion equations, and some achievements have been obtained. However, to our knowledge of the literature, up to date, all the proposed numerical methods to modify FODE have achieved at most a second-order time accuracy. In this study, we focus mainly on the numerical methods based on numerical integration in order to modify the fractional-order diffusion equation:1+1/12δx2pjk−1+1/12δx2pjk−1=μα∑l=0kλαlδx2pjk−l+μβ∑l=0kλβlδx2pjk−l+τ/21+1/12δx2fjk−1+fjk, k=1,2,…,K;j=1,2,…,j−1,pj0=ωj,j=0,1,…,J, p0k=φtk,pjk=ψtk,k=0,1,…,K,fjl=fxj,tl,ωj=ωxj. Accordingly, numerical methods can be built to modify FODE with second-order time accuracy and fourth-order spatial accuracy in ∂px,t/∂t=∂1−α/∂t1−α+B∂1−β/∂t1−β∂2px,t/∂x2+fx,t,0<t≤1,0<x<1,px,0=0,0≤x≤1. p0,t=t2,p1,t=et2,0≤t≤1. Our suggested method can improve the time precision with a certain value.


---

## Body

## 1. Introduction

In recent years, fractional-order diffusion equations (FODEs) have been widely employed to describe anomalous diffusion processes in many physical and biological systems. By acting upon the diffusion factor with two levels of fractional-order time derivatives, models of special diffusion phenomena have been proposed [1–3].(1)∂px,t∂t=A∂1−α∂t1−α+B∂1−β∂t1−β∂2px,t∂t2.Some scholars have constructed different numerical methods to solve the modified FODE, and accordingly, some useful achievements have been obtained. For one-dimensional positive FODE, Langlands [4] proposed an interpretation with an infinite series form of the Fox function over an infinite region. Liu et al. [5, 6] discussed numerical methods and analytical techniques to develop a finite element approximation with first-order time accuracy and mth order spatial accuracy, where m is the number of segmented polynomials. Mohebbi et al. [7] applied a fourth-order compact formula for the second-order spatial partial derivatives and the discretization of Riemann-Liouville fractional-order time derivatives to provide a higher-order and absolutely stable format with first-order time accuracy and fourth-order spatial accuracy. Bhrawy [8] studied the compact subdiffusion scheme including second-order time accuracy and fourth-order spatial accuracy. For the multirepair positive FODE, Zhang et al. [9] developed a finite difference/finite element method with (1 + min{α, β}) order time accuracy and mth order spatial accuracy, where m represents the number of segmented polynomials. Mohebbi et al. [10] proposed a fourth-order compact solution method with first-order time accuracy and fourth-order spatial accuracy. Abbaszadeh and Mohebbi [11] discussed a solution obtained by a radial basis function (RBF) meshless method with min{α, β} order time accuracy. Wang and Wang [12] analyzed the tight LOD method and its extrapolation method including 2 min{α, β} order time accuracy and fourth-order spatial accuracy. For two-dimensional variable-order modified FODE, Chen and Liu [13] examined a numerical method with first-order time accuracy and fourth-order spatial accuracy and further developed the numerical method to improve the time accuracy.According to the literature, the existing numerical methods for the modified FODE can achieve at most second-order time accuracy. However, this study develops a new numerical method with not only second-order time accuracy but also fourth-order spatial accuracy based on the numerical integration of the modified FODE.


### 1.1. Basic Concepts and Properties

Definition 1.
(Grünwald–Letnikov fractional stratification number). Letα be a positive real number, and n−1≤α<n, n is a positive integer, and let the function f (x) be defined on the interval [a, b] as(2)Dαxαfx=limh⟶0∑j=0x−a/h−1jαjfx−jh.
It is theα-order Grünwald–Letnikov (G-L) fractional-order derivative of function f (x), [z] is the largest integer that does not exceeds z, and(3)αj=αα−1…α−j+1j!.Definition 2 (Riemann-Liouville fractional stratification number).
Letα be a positive real number, and n−1≤α<n, n is a positive integer, and let the function fx be defined on the interval [a, b], respectively. Here,(4)Dαxαfx=1Гn−αdndxn∫αxfτdτx−τα−n−1,Dαbαfx=−1nГn−αdndxn∫αxfτdττ−xα−n−1,represent the left and right α-order Riemann-Liouville fractional derivatives.Property 1.
Letα be a positive real number, and n − 1 ≤ α < n, n is a positive integer defined in the interval of [a, b], the function fx has up to n − 1 continuous functions, and fnx is integrable in [a, b]. Then, the Riemann-Liouville fractional derivative is equivalent to the Grünwald–Letnikov (GL) fractional derivative.


### 1.2. Construction of Numerical Methods

In this research, the following numerical method is developed to modify the initial and boundary values of the fractional diffusion equation (MFDE):(5)∂px,t∂t=∂1−α∂t1−α+B∂1−β∂t1−β∂2px,t∂x2+fx,t,(6)px,0=ωx,0≤x≤L,(7)p0,t=φt,pL,t=ψt,0≤t≤T,where 0<α<β<1,D0t1−γPx,t is the fractional large derivative of 1−γ.(8)D0t1−γPx,t=1Гγ∂∂t∫0tpx,st−s1−γds.Supposepx,t∈PΩ,px,t is the exact solution of problems (1)–(6) and ∂2fx,t/∂t2∈CΩ, where(9)Ω=x,t|0≤x≤L,0≤t≤T,PΩ=px,t|∂6px,t∂x6,∂4px,t∂x2∂t2∈Ω.Suppose that(10)xj=jh,j=0,1,…,J,tk=kτ,k=0,1,…,K,where h=L/j and τ=T/K are the space step and time step, respectively.We define(11)δx2νji=νj−1i−2vji+vj+1i,and integrating the two sides of equation (5) with respect to t on the interval tk−1,tk, we get(12)Pxj,tk−Pxj,tk−1=1Гα∫0tk∂2pxj,s∂x2dstk−s1−α−∫0tk−1∂2pxj,s∂x2dstk−1−s1−α+1Гβ∫0tk∂2pxj,s∂x2dstk−s1−β−∫0tk−1∂2pxj,s∂x2dstk−1−s1−β+∫tk−1tkfxj,sds.Hence, we get(13)1+112δx2pxj,tk−1+112δx2pxj,tk−1=1+112δx21Гα∫0tk∂2pxj,s∂x2dstk−s1−α−∫0tk−1∂2pxj,s∂x2dstk−1−s1−α+1+112δx21Гβ∫0tk∂2pxj,s∂x2dstk−s1−β−∫0tk−1∂2pxj,s∂x2dstk−1−s1−β+1+112δx2∫tk−1tkfxj,sds.Sincepx,t∈PΩ, then(14)1+112δx2∂2pxj,ti∂x2=δx2pxj,tih2+Oh4,∂2pxj,s∂x2=tl−sτ∂2pxj,tl−1∂x2+s−tl−1τ∂2pxj,tl∂x2+Oτ2,tl−1≤s≤tl,hence,(15)1+112δx21Гα∫0tk∂2pxj,s∂x2dstk−s1−α−∫0tk−1∂2pxj,s∂x2dstk−1−s1−α=1+112δx21Гα∑l=1k∫tl−1tltl−sτ∂2pxj,tl−1∂x2+s−tl−1τ∂2pxj,tl∂x2+Oτ2dstk−s1−α−1+112δx21Гα∑l=1k−1∫tl−1tltl−sτ∂2pxj,tl−1∂x2+s−tl−1τ∂2pxj,tl∂x2+Oτ2dstk−1−s1−α=1Гα∑l=1k∫tl−1tltl−sτ1+112δx2∂2pxj,tl−1∂x2+s−tl−1τ1+112δx2∂2pxj,tl∂x2+Oτ2dstk−s1−α−1Гα∑l=1k−1∫tl−1tltl−sτ1+112δx2∂2pxj,tl−1∂x2+s−tl−1τ1+112δx2∂2pxj,tl∂x2+Oτ2dstk−1−s1−α=1Гα∑l=1k∫tl−1tltl−sτδx2pxj,tl−1h2+Oh4+Oτ2dstk−1−s1−α−1Гα∑l=1k∫tl−1tltl−sτδx2pxj,tl−1h2+Oh4+s−tl−1τδx2pxj,tlh2+Oh4+Oτ2dstk−1−s1−α=μα∑l=1kλαlδx2pxj,tk−l+R¯jk,where(16)μα=ταΓ2+αh2,R¯jk=Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−α−∑l=1k−1∫tl−1tldstk−1−s1−α,λl=λ=0,21+ϒ−3,λ=1,−l−21+ϒ+3l−21+ϒ−l−21+ϒ+l+11+ϒ,λ=2,…,k−1,1+Υkϒ−k−1ϒ−k−21+ϒ+2k−11+ϒ−k1+ϒ,λ=k.Similarly, then, we have(17)1+112δx21Гβ∫0tk∂2pxj,s∂x2dstk−s1−β−∫0tk−1∂2pxj,s∂x2dstk−1−s1−β=μβ∑l=0kλβlδx2pxj,tk−l+R˜jk,where(18)μβ=τβΓ2+βh2,R˜jk=Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−β−∑l=1k−1∫tl−1tldstk−s1−β.Because∂2fx,t/∂t2∈CΩ, the following trapezoidal formula holds:(19)∫tk−1tkfxj,sds=τ2fxj,tk−1+fxj,tk+Oτ3,and according to the above analysis, we have(20)1+112δx2pxj,tk−1+112δx2pxj,tk−1=μα∑l=0kλαlδx2pxj,tk−l+μβ∑l=0kλβlδx2pxj,tk−l+τ21+112δx2fxj,tk−1+fxj,tk+Rjk,where(21)Rjk=R¯jk+R˜jk+Oτ3=Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−α−∑l=1k−1∫tl−1tldstk−1−s1−α+Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−β−∑l=1k−1∫tl−1tldstk−1−s1−β+Oτ3.Please note that(22)∑l=1k∫tl−1tldstk−s1−α=∫0tkdstk−s1−α=kταα,∑l=1k−1∫tl−1tldstk−1−s1−α=∫0tk−1dstk−1−s1−α=k−1ταα,∑l=1k∫tl−1tldstk−s1−β=∫0tkdstk−s1−β=kτββ,∑l=1k−1∫tl−1tldstk−1−s1−β=∫0tk−1dstk−1−s1−β=k−1τββ,with kτ≤T; then,(23)Rjk=Oτ2+h4.Now, the numerical methods should be taken into account for solving problems (5)–(7):(24)1+112δx2pjk−1+112δx2pjk−1=μα∑l=0kλαlδx2pjk−l+μβ∑l=0kλβlδx2pjk−l+τ21+112δx2fjk−1+fjk,(25)k=1,2,…,K;j=1,2,…,j−1,pj0=ωj,j=0,1,…,J,(26)p0k=φtk,pjk=ψtk,k=0,1,…,K,where fjl=fxj,tl,ωj=ωxj, and pjl is the approximation of the exact solution pxj,tl.


### 1.3. Convergence of Numerical Methods

In this section, we discuss the convergence of numerical equations (24)–(26). Subtract (24) from (20) and then obtain the following error equation:(27)1+112δx2Ejk−1+112δx2Ejk−1=μα∑l=0kλαlδx2Ejk−l+μβ∑l=0kλβlδx2Ejk−l+Rjk,k=1,2,…,K;j=1,2,…,j−1,where(28)Ejk=pxj,tk−pjk.Fork=0,1,…,K, the following grid functions are defined, respectively:(29)Ekx=Ejk,x∈xj−1/2,xj+1/2,j=1,2,…,j−1,0,x∈0,h2∪L−H2,L,Rkx=Rjk,x∈xj−1/2,xj+1/2,j=1,2,…,j−1,0,x∈0,h2∪L−H2,L,then, Ekx and Rkx can be expanded with the following Fourier series, respectively:(30)Ekx=∑−∞∞ξklei2πlx/L,k=0,1,…,K,Rkx=∑−∞∞ηklei2πlx/L,k=0,1,…,K,where(31)ξkl=∫0LEkle−i2πlx/Ldx,ηkl=∫0LRkle−i2πlx/Ldx.Let(32)Ekx=E1k,E2k,…,Ej−1kT,Rkx=R1k,R2k,…,Rj−1kT.The following Parseval equation can be derived:(33)∫0LEkx2dx=∫0LEkx⋅Ekx¯dx=∫0LEkx⋅∑−∞∞ξklei2πlx/L¯dx=∫0LEkx⋅∑−∞∞ξkl¯e−i2πlx/Ldx=∑−∞∞ξkl¯∫0LEkxe−i2πlx/Ldx=∑−∞∞ξkl¯Lξkl=L∑−∞∞ξkl¯ξkl=L∑−∞∞ξkl2.Similarly, the following Parseval equation can also be derived:(34)∫0LRkx2dx=L∑−∞∞ηkl2.Notice that(35)∫0LEkx2dx=∑j=1j−1hEjk21/2,∫0LRkx2dx=∑j=1j−1hRjk21/2.So far, we obtain(36)Ek2≡∑j=1j−1hEjk21/2=L∑l=−∞∞ξkl21/2,k=0,1,…,K,(37)Rk2≡∑j=1j−1hRjk21/2=L∑l=−∞∞ηkl21/2,k=0,1,…,K.Assume, respectively,(38)Ejk=ξkeiσjh,Rjk=ηkeiσjh,where σ=2π/L. Substituting the above expression into (27), we get(39)1−13sin2σh2ξk−1−13sin2σh2ξk−1=−4μαsin2σh2∑l=0kλαlξk−l−4μβsin2σh2∑l=0kλβlξk−l+ηk,k=1,2,…,K.After rearrangement,(40)ξk=1ΦΨξk−l−4μαsin2σh2∑l=0kλαlξk−l4μβsin2σh2∑l=0kλβlξk−l+ηk,k=1,2,…,K,where(41)Φ=1−13sin2σh2+4μα+μβsin2σh2,ψ=1−13sin2σh2−4μαλα1+μβλβ1sin2σh2.Similar to reference [14], we also obtain the following lemma:Lemma 1.
If0<ϒ<1, the coefficient λϒll=0,1,…,k has the following properties:(1)
λϒ0=1,λϒ1=21+ϒ−3, λϒl<0, l=2,…,k;(2)
∑l=0kλϒl=1+ϒkϒ−k−1>0;(3)
−1<−λϒl<1, −∑l=2kλϒl<1+λϒ1.Lemma 2.
Set 0 <α < β < 1, and the time and space steps satisfy τ2=h4; then,(42)Φ≥1.Proof.
From0<α<β<1 and τ2=h4, the following equation is built.(43)μα+μβ=ταΓ2+αh2+τβΓ2+βh2=h2α−1Γ2+α+h2β−1Γ2+β≥1Γ3+1Γ3=1.Lemma 3.
If0<α<β<1,μαλα1+μβλβ1≤0, or μαλα1+μβλβ1≤1/6, then(44)Φ≥0.
From (23), therefore, constant C1 exists and makes(45)Rjk≤C1τ2+h4,k=1,2,…,K.
Again, paying attention to the first equation in (37), we have(46)RK≤C1Lτ2+h4,k=1,2,…,K.
By the convergence of the series at the left end of (37), there is a positive constant C2, so that(47)ηjk≡ηKl≤C2τη1l≡C2τη1,k=1,2,…,K.Theorem 1.
Assuming thatpx,t∈PΩ is the exact solution of problems (5)–(7), ∂2fx,t/∂t2∈CΩ, and time and space steps satisfy τ2=h4, if μαλα1+μβλβ1≤0, or μαλα1+μβλβ1≤1/6, then the numerical equations (24)–(26) converge, and the convergence order is Oτ2+h4.Proof.
First, we prove the following conclusion, where there is a positive constantC, such that(48)ξk≤Ckη1,k=1,2,…,K,where ξkk=1,2,…,K is the solution of (40).
Fork = 1, since E0=0, get ξ0=0, and further applying Lemma 2 and (47), we can obtain from equation (40),(49)ξ1=1Φη1≤η1≤C2τη1.
Assume(50)ξk≤C2nτη1,k=1,2,…,K−1.
According to Lemmas1–3, we have through (40),(51)ξ1=1ΦΨξk−1−4μαsin2σh2∑l=2kλαlξk−l−4μβsin2σh2∑l=2kλβlξk−l+ηk≤1ΦΨξk−1−4μαsin2σh2∑l=2kλαlξk−l−4μβsin2σh2∑l=2kλβlξk−l+ηk≤1ΦΨk−1+4μαsin2σh2∑l=2kλαlk−1+4μβsin2σh2∑l=2kλβlk−1+1C2τηk≤1ΦΨ+4μαsin2σh2∑l=2kλαl+4μβsin2σh2∑l=2kλβlK−1+1C2τηk=1ΦΨ−4μαsin2σh2∑l=2kλαl−4μβsin2σh2∑l=2kλβlK−1+1C2τηk≤1ΦΨ−4μα1+λαlsin2σh2+4μβ1+λβlsin2σh2K−1+1C2τηk=K−1+1ΦC2τηk≤C2kτη1.
Conclusion (48) is proven by induction. From conclusions (36), (37), (45), and (48), we obtain the following one:(52)Ek2≤C2kτR12≤Cτ2+h2,where C=C1C2LT.
We know that because there are many high-precision numerical integration formulas (45) and (46), it is not difficult to continue to build up a high time precision numerical method for solving the modified fractional diffusion equation (5). However, the problem is that the qualitative analysis of convergence and stability can be very difficult.


## 1.1. Basic Concepts and Properties

Definition 1.
(Grünwald–Letnikov fractional stratification number). Letα be a positive real number, and n−1≤α<n, n is a positive integer, and let the function f (x) be defined on the interval [a, b] as(2)Dαxαfx=limh⟶0∑j=0x−a/h−1jαjfx−jh.
It is theα-order Grünwald–Letnikov (G-L) fractional-order derivative of function f (x), [z] is the largest integer that does not exceeds z, and(3)αj=αα−1…α−j+1j!.Definition 2 (Riemann-Liouville fractional stratification number).
Letα be a positive real number, and n−1≤α<n, n is a positive integer, and let the function fx be defined on the interval [a, b], respectively. Here,(4)Dαxαfx=1Гn−αdndxn∫αxfτdτx−τα−n−1,Dαbαfx=−1nГn−αdndxn∫αxfτdττ−xα−n−1,represent the left and right α-order Riemann-Liouville fractional derivatives.Property 1.
Letα be a positive real number, and n − 1 ≤ α < n, n is a positive integer defined in the interval of [a, b], the function fx has up to n − 1 continuous functions, and fnx is integrable in [a, b]. Then, the Riemann-Liouville fractional derivative is equivalent to the Grünwald–Letnikov (GL) fractional derivative.


## 1.2. Construction of Numerical Methods

In this research, the following numerical method is developed to modify the initial and boundary values of the fractional diffusion equation (MFDE):(5)∂px,t∂t=∂1−α∂t1−α+B∂1−β∂t1−β∂2px,t∂x2+fx,t,(6)px,0=ωx,0≤x≤L,(7)p0,t=φt,pL,t=ψt,0≤t≤T,where 0<α<β<1,D0t1−γPx,t is the fractional large derivative of 1−γ.(8)D0t1−γPx,t=1Гγ∂∂t∫0tpx,st−s1−γds.Supposepx,t∈PΩ,px,t is the exact solution of problems (1)–(6) and ∂2fx,t/∂t2∈CΩ, where(9)Ω=x,t|0≤x≤L,0≤t≤T,PΩ=px,t|∂6px,t∂x6,∂4px,t∂x2∂t2∈Ω.Suppose that(10)xj=jh,j=0,1,…,J,tk=kτ,k=0,1,…,K,where h=L/j and τ=T/K are the space step and time step, respectively.We define(11)δx2νji=νj−1i−2vji+vj+1i,and integrating the two sides of equation (5) with respect to t on the interval tk−1,tk, we get(12)Pxj,tk−Pxj,tk−1=1Гα∫0tk∂2pxj,s∂x2dstk−s1−α−∫0tk−1∂2pxj,s∂x2dstk−1−s1−α+1Гβ∫0tk∂2pxj,s∂x2dstk−s1−β−∫0tk−1∂2pxj,s∂x2dstk−1−s1−β+∫tk−1tkfxj,sds.Hence, we get(13)1+112δx2pxj,tk−1+112δx2pxj,tk−1=1+112δx21Гα∫0tk∂2pxj,s∂x2dstk−s1−α−∫0tk−1∂2pxj,s∂x2dstk−1−s1−α+1+112δx21Гβ∫0tk∂2pxj,s∂x2dstk−s1−β−∫0tk−1∂2pxj,s∂x2dstk−1−s1−β+1+112δx2∫tk−1tkfxj,sds.Sincepx,t∈PΩ, then(14)1+112δx2∂2pxj,ti∂x2=δx2pxj,tih2+Oh4,∂2pxj,s∂x2=tl−sτ∂2pxj,tl−1∂x2+s−tl−1τ∂2pxj,tl∂x2+Oτ2,tl−1≤s≤tl,hence,(15)1+112δx21Гα∫0tk∂2pxj,s∂x2dstk−s1−α−∫0tk−1∂2pxj,s∂x2dstk−1−s1−α=1+112δx21Гα∑l=1k∫tl−1tltl−sτ∂2pxj,tl−1∂x2+s−tl−1τ∂2pxj,tl∂x2+Oτ2dstk−s1−α−1+112δx21Гα∑l=1k−1∫tl−1tltl−sτ∂2pxj,tl−1∂x2+s−tl−1τ∂2pxj,tl∂x2+Oτ2dstk−1−s1−α=1Гα∑l=1k∫tl−1tltl−sτ1+112δx2∂2pxj,tl−1∂x2+s−tl−1τ1+112δx2∂2pxj,tl∂x2+Oτ2dstk−s1−α−1Гα∑l=1k−1∫tl−1tltl−sτ1+112δx2∂2pxj,tl−1∂x2+s−tl−1τ1+112δx2∂2pxj,tl∂x2+Oτ2dstk−1−s1−α=1Гα∑l=1k∫tl−1tltl−sτδx2pxj,tl−1h2+Oh4+Oτ2dstk−1−s1−α−1Гα∑l=1k∫tl−1tltl−sτδx2pxj,tl−1h2+Oh4+s−tl−1τδx2pxj,tlh2+Oh4+Oτ2dstk−1−s1−α=μα∑l=1kλαlδx2pxj,tk−l+R¯jk,where(16)μα=ταΓ2+αh2,R¯jk=Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−α−∑l=1k−1∫tl−1tldstk−1−s1−α,λl=λ=0,21+ϒ−3,λ=1,−l−21+ϒ+3l−21+ϒ−l−21+ϒ+l+11+ϒ,λ=2,…,k−1,1+Υkϒ−k−1ϒ−k−21+ϒ+2k−11+ϒ−k1+ϒ,λ=k.Similarly, then, we have(17)1+112δx21Гβ∫0tk∂2pxj,s∂x2dstk−s1−β−∫0tk−1∂2pxj,s∂x2dstk−1−s1−β=μβ∑l=0kλβlδx2pxj,tk−l+R˜jk,where(18)μβ=τβΓ2+βh2,R˜jk=Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−β−∑l=1k−1∫tl−1tldstk−s1−β.Because∂2fx,t/∂t2∈CΩ, the following trapezoidal formula holds:(19)∫tk−1tkfxj,sds=τ2fxj,tk−1+fxj,tk+Oτ3,and according to the above analysis, we have(20)1+112δx2pxj,tk−1+112δx2pxj,tk−1=μα∑l=0kλαlδx2pxj,tk−l+μβ∑l=0kλβlδx2pxj,tk−l+τ21+112δx2fxj,tk−1+fxj,tk+Rjk,where(21)Rjk=R¯jk+R˜jk+Oτ3=Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−α−∑l=1k−1∫tl−1tldstk−1−s1−α+Oτ2+Oh4∑l=1k∫tl−1tldstk−s1−β−∑l=1k−1∫tl−1tldstk−1−s1−β+Oτ3.Please note that(22)∑l=1k∫tl−1tldstk−s1−α=∫0tkdstk−s1−α=kταα,∑l=1k−1∫tl−1tldstk−1−s1−α=∫0tk−1dstk−1−s1−α=k−1ταα,∑l=1k∫tl−1tldstk−s1−β=∫0tkdstk−s1−β=kτββ,∑l=1k−1∫tl−1tldstk−1−s1−β=∫0tk−1dstk−1−s1−β=k−1τββ,with kτ≤T; then,(23)Rjk=Oτ2+h4.Now, the numerical methods should be taken into account for solving problems (5)–(7):(24)1+112δx2pjk−1+112δx2pjk−1=μα∑l=0kλαlδx2pjk−l+μβ∑l=0kλβlδx2pjk−l+τ21+112δx2fjk−1+fjk,(25)k=1,2,…,K;j=1,2,…,j−1,pj0=ωj,j=0,1,…,J,(26)p0k=φtk,pjk=ψtk,k=0,1,…,K,where fjl=fxj,tl,ωj=ωxj, and pjl is the approximation of the exact solution pxj,tl.


## 1.3. Convergence of Numerical Methods

In this section, we discuss the convergence of numerical equations (24)–(26). Subtract (24) from (20) and then obtain the following error equation:(27)1+112δx2Ejk−1+112δx2Ejk−1=μα∑l=0kλαlδx2Ejk−l+μβ∑l=0kλβlδx2Ejk−l+Rjk,k=1,2,…,K;j=1,2,…,j−1,where(28)Ejk=pxj,tk−pjk.Fork=0,1,…,K, the following grid functions are defined, respectively:(29)Ekx=Ejk,x∈xj−1/2,xj+1/2,j=1,2,…,j−1,0,x∈0,h2∪L−H2,L,Rkx=Rjk,x∈xj−1/2,xj+1/2,j=1,2,…,j−1,0,x∈0,h2∪L−H2,L,then, Ekx and Rkx can be expanded with the following Fourier series, respectively:(30)Ekx=∑−∞∞ξklei2πlx/L,k=0,1,…,K,Rkx=∑−∞∞ηklei2πlx/L,k=0,1,…,K,where(31)ξkl=∫0LEkle−i2πlx/Ldx,ηkl=∫0LRkle−i2πlx/Ldx.Let(32)Ekx=E1k,E2k,…,Ej−1kT,Rkx=R1k,R2k,…,Rj−1kT.The following Parseval equation can be derived:(33)∫0LEkx2dx=∫0LEkx⋅Ekx¯dx=∫0LEkx⋅∑−∞∞ξklei2πlx/L¯dx=∫0LEkx⋅∑−∞∞ξkl¯e−i2πlx/Ldx=∑−∞∞ξkl¯∫0LEkxe−i2πlx/Ldx=∑−∞∞ξkl¯Lξkl=L∑−∞∞ξkl¯ξkl=L∑−∞∞ξkl2.Similarly, the following Parseval equation can also be derived:(34)∫0LRkx2dx=L∑−∞∞ηkl2.Notice that(35)∫0LEkx2dx=∑j=1j−1hEjk21/2,∫0LRkx2dx=∑j=1j−1hRjk21/2.So far, we obtain(36)Ek2≡∑j=1j−1hEjk21/2=L∑l=−∞∞ξkl21/2,k=0,1,…,K,(37)Rk2≡∑j=1j−1hRjk21/2=L∑l=−∞∞ηkl21/2,k=0,1,…,K.Assume, respectively,(38)Ejk=ξkeiσjh,Rjk=ηkeiσjh,where σ=2π/L. Substituting the above expression into (27), we get(39)1−13sin2σh2ξk−1−13sin2σh2ξk−1=−4μαsin2σh2∑l=0kλαlξk−l−4μβsin2σh2∑l=0kλβlξk−l+ηk,k=1,2,…,K.After rearrangement,(40)ξk=1ΦΨξk−l−4μαsin2σh2∑l=0kλαlξk−l4μβsin2σh2∑l=0kλβlξk−l+ηk,k=1,2,…,K,where(41)Φ=1−13sin2σh2+4μα+μβsin2σh2,ψ=1−13sin2σh2−4μαλα1+μβλβ1sin2σh2.Similar to reference [14], we also obtain the following lemma:Lemma 1.
If0<ϒ<1, the coefficient λϒll=0,1,…,k has the following properties:(1)
λϒ0=1,λϒ1=21+ϒ−3, λϒl<0, l=2,…,k;(2)
∑l=0kλϒl=1+ϒkϒ−k−1>0;(3)
−1<−λϒl<1, −∑l=2kλϒl<1+λϒ1.Lemma 2.
Set 0 <α < β < 1, and the time and space steps satisfy τ2=h4; then,(42)Φ≥1.Proof.
From0<α<β<1 and τ2=h4, the following equation is built.(43)μα+μβ=ταΓ2+αh2+τβΓ2+βh2=h2α−1Γ2+α+h2β−1Γ2+β≥1Γ3+1Γ3=1.Lemma 3.
If0<α<β<1,μαλα1+μβλβ1≤0, or μαλα1+μβλβ1≤1/6, then(44)Φ≥0.
From (23), therefore, constant C1 exists and makes(45)Rjk≤C1τ2+h4,k=1,2,…,K.
Again, paying attention to the first equation in (37), we have(46)RK≤C1Lτ2+h4,k=1,2,…,K.
By the convergence of the series at the left end of (37), there is a positive constant C2, so that(47)ηjk≡ηKl≤C2τη1l≡C2τη1,k=1,2,…,K.Theorem 1.
Assuming thatpx,t∈PΩ is the exact solution of problems (5)–(7), ∂2fx,t/∂t2∈CΩ, and time and space steps satisfy τ2=h4, if μαλα1+μβλβ1≤0, or μαλα1+μβλβ1≤1/6, then the numerical equations (24)–(26) converge, and the convergence order is Oτ2+h4.Proof.
First, we prove the following conclusion, where there is a positive constantC, such that(48)ξk≤Ckη1,k=1,2,…,K,where ξkk=1,2,…,K is the solution of (40).
Fork = 1, since E0=0, get ξ0=0, and further applying Lemma 2 and (47), we can obtain from equation (40),(49)ξ1=1Φη1≤η1≤C2τη1.
Assume(50)ξk≤C2nτη1,k=1,2,…,K−1.
According to Lemmas1–3, we have through (40),(51)ξ1=1ΦΨξk−1−4μαsin2σh2∑l=2kλαlξk−l−4μβsin2σh2∑l=2kλβlξk−l+ηk≤1ΦΨξk−1−4μαsin2σh2∑l=2kλαlξk−l−4μβsin2σh2∑l=2kλβlξk−l+ηk≤1ΦΨk−1+4μαsin2σh2∑l=2kλαlk−1+4μβsin2σh2∑l=2kλβlk−1+1C2τηk≤1ΦΨ+4μαsin2σh2∑l=2kλαl+4μβsin2σh2∑l=2kλβlK−1+1C2τηk=1ΦΨ−4μαsin2σh2∑l=2kλαl−4μβsin2σh2∑l=2kλβlK−1+1C2τηk≤1ΦΨ−4μα1+λαlsin2σh2+4μβ1+λβlsin2σh2K−1+1C2τηk=K−1+1ΦC2τηk≤C2kτη1.
Conclusion (48) is proven by induction. From conclusions (36), (37), (45), and (48), we obtain the following one:(52)Ek2≤C2kτR12≤Cτ2+h2,where C=C1C2LT.
We know that because there are many high-precision numerical integration formulas (45) and (46), it is not difficult to continue to build up a high time precision numerical method for solving the modified fractional diffusion equation (5). However, the problem is that the qualitative analysis of convergence and stability can be very difficult.


## 2. Stability of the Numerical Methods

The stability of numerical format was discussed in (24)–(26); consider the following difference equation:(53)1+112δx2ρjk−1+112δx2ρjk−1=μα∑l=0kλαlδx2ρjk−l+μβ∑l=0kλβlδx2ρjk−l+Rjkk=1,2,…,K;j=1,2,…,j−1,where ρjk=pjk−p˜jk, p˜jk is the approximation of pjk.Fork = 0, 1, …, K, respectively, define the following grid function:(54)ρkx=ρjk,x∈xj−1/2,xj+1/2,j=1,2,…,j−1,0,x∈0,h2∪L−h2,L.Then,ρkx can be expanded by Fourier series:(55)ρkx=∑−∞∞ζklei2πlx/L,k=0,1,…,K,where(56)ζkl=1L∫0Lρkle−i2πlx/Ldx.Let(57)ρkx=ρ1k,ρ2k,…,ρj−1kT.Similar to the derivation of (36), the following Parseval equation can also be obtained:(58)ρk2≡∑j=1j−1hρjk21/2=L∑l=−∞∞ζkl21/2,k=0,1,…,K.Assume(59)Ejk=ζkeiσjh,where σ=2π/L.Substituting the above expression into (53), we have(60)1−13sin2σh2ζk−1−13sin2σh2ζk−1=−4μαsin2σh2∑l=0kλαlζk−l−4μβsin2σh2∑l=0kλβlζk−l,k=1,2,…,K.After rearranging, we get(61)ζk=1ΦΨζk−l−4μαsin2σh2∑l=0kλαlζk−l4μβsin2σh2∑l=0kλβlζk−l,k=1,2,…,K.Theorem 2.
Assume that the time and space steps satisfyτ2=h4, if μαλα1+μβλβ1≤0, or 0<μαλα1+μβλβ1≤1/6, then numerical equations (24)–(26) are stable.Proof.
First, the following conclusion is proved:(62)ζk≤ζ0,k=1,2,…,K,where ζkk=1,2,…,K is the solution of (61).
Whenk=1, from Lemmas 2 and 3, we get from equation (61)(63)ζ1=1ΦΨζ0≤ζ0,and assume(64)ζn≤ζ0,n=1,2,…,K−1.
According to Lemmas1–3, the following equation is obtained from (61):(65)ζ1=1ΦΨζk−1−4μαsin2σh2∑l=2kλαlζk−l−4μβsin2σh2∑l=2kλβlζk−l≤1ΦΨζk−1+4μαsin2σh2∑l=2kλαlζk−l+4μβsin2σh2∑l=2kλβlζk−l≤1ΦΨ+4μαsin2σh2∑l=2kλαl+4μβsin2σh2∑l=2kλβlζ0=1ΦΨ−4μαsin2σh2∑l=2kλαl−4μβsin2σh2∑l=2kλβlζ0≤1ΦΨ+4μα1+λαlsin2σh2+4μβ1+λβlsin2σh2ζ0=ζ0.
Conclusion (62) can be achieved by induction. From conclusion (58) and (62), we obtain the following one:(66)ρk2≤ρ02,k=1,2,…,K.


## 3. Solvability of Numerical Methods

This section discusses the solvability of numerical equations (24)–(26). The corresponding chi-square linear algebraic equations with numerical equations (24)–(26) are(67)1+112δx2ρjk−1+112δx2ρjk−1,(68)=μα∑l=0kλαlδx2ρjk−l+μβ∑l=0kλβlδx2ρjk−l,(69)k=1,2,…,K;j=1,2,…,j−1,pj0=0,j=0,1,…,J,(70)p0k=0,pjk=0,k=0,1,…,K.Similar to the proof of Theorem2, it can also be shown that the solution pk of the chi-square linear equations (67)–(71) satisfies(71)ρk2≤ρ02,k=1,2,…,K,from p0=0, we have(72)p0k=0,pjk=0,k=0,1,…,K.This means that the numerical equations (24)–(26) have only zero solutions for the corresponding linear algebraic equations of the same order. Thus, the following conclusion is drawn.Theorem 3.
Suppose that the time and space steps satisfyτ2=h4, if μαλα1+μβλβ1≤0, or 0<μαλα1+μβλβ1≤1/6, then the numerical equations (24)–(26) are uniquely solvable.


## 4. Numerical Tests

The numerical equations (24)–(26) are used to solve the following modified fractional-order diffusion equation for the initial value problem:(73)∂px,t∂t=∂1−α∂t1−α+B∂1−β∂t1−β∂2px,t∂x2+fx,t,(74)0<t≤1,0<x<1,px,0=0,0≤x≤1,(75)p0,t=t2,p1,t=et2,0≤t≤1,fx,t=2ext−t1+∝Γ2+α−t1+βΓ2+β,where fx,t=2ext−t1+∝/Γ2+α−t1+β/Γ2+β, and the exact solution of problems (73)–(75) is(76)px,t=ext2.Let the maximum error of the numerical solution be(77)E∞=max0≤k≤KEk2=Oτa+hb.Assume that the maximum errors of spatial step, time step, and numerical solution arehf, τf, and E∞f in the previous experiment and the maximum error of the spatial step, time step, and numerical solution are hs, τs, and E∞s in the later experiment. Then, the time accuracy a and spatial accuracy b are calculated by the following equation:(78)a=logτs/τfE∞sE∞f,b=loghs/hfE∞sE∞f.Table1 provides that the numerical equations (24)–(26) for problems (73)–(75) are obtained for different α, β, andτ2 = h2 with maximum error E∞, time accuracy α, and spatial accuracy β, respectively.Table 1
Maximum errorE∞, time accuracy, and spatial accuracy b of numerical equations (24)–(26) for solving problems (73)–(75).
αβτ2=h2=1/16τ2=h4=1/81abτ2=h4=1/256ab0.10.73.5759 × 10−57.4432 × 10−61.943.872.4953 × 10−61.903.800.20.83.5591 × 10−57.4432 × 10−61.933.862.4953 × 10−61.903.800.30.93.5422 × 10−57.3461 × 10−61.943.882.2508 × 10−62.064.110.40.53.5675 × 10−57.4432 × 10−61.933.872.4272 × 10−61.953.900.50.63.5505 × 10−57.3461 × 10−61.943.892.4272 × 10−61.923.850.60.73.5253 × 10−57.2489 × 10−61.953.902.4272 × 10−61.903.800.70.713.5169 × 10−57.2489 × 10−61.953.902.2508 × 10−62.034.070.80.823.4916 × 10−57.3461 × 10−61.933.842.4272 × 10−61.923.850.90.933.4748 × 10−57.1571 × 10−61.953.902.0066 × 10−62.214.42


## 5. Conclusion

As given in Table1, our theoretical analysis results are strongly supported by the numerical experiments.Figure1 shows that the numerical methods (24)–(26) are used to compare the numerical solutions (N.S.) of problems (73)–(75) with the exact solutions (E.S.) of problems (73)–(75) in t=1/4,2/4,3/4,1, when α=0.5, β=0.6, τ2=h4=1/256.Figure 1
Whenα=0.5, β=0.6, and τ2=h4=1/256, numerical methods (24)–(26) are utilized to compare the numerical solutions (N.S.) of problems (73)–(75) with the exact solutions (E.S.) of problems (73)–(75) in t=1/4,2/4,3/4,1.Figure2 shows that the numerical methods (24)–(26) are used to compare the numerical solutions (N.S.) of problems (73)–(75) with the exact solutions (E.S.) of problems (73)–(75) in t=1/4,2/4,3/4,1, when α=0.5, β=0.6, τ2=h4=1/256.Figure 2
Whenα=0.5, β=0.6, and τ2=h4=1/256, numerical methods (24)–(26) are utilized to compare the numerical solutions (N.S.) of problems (73)–(75) with the exact solutions (E.S.) of problems (73)–(75) in t=1/4,2/4,3/4,1.As shown in Figures1 and 2, the numerical solutions of problems (73)–(75) and the exact solutions of problems (73)–(75) obtained by numerical methods (24)–(26) have a good approximation effect. This finding revealed that our theoretical analysis results are reliable.In conclusion, one-dimensional fractional diffusion equations were studied in the present study, with the purpose that the construction techniques of numerical methods for solving one-dimensional fractional diffusion equations and the corresponding numerical analysis can be extended to multidimensional fractional diffusion equations.



---

*Source: 4846747-2022-03-07.xml*