# Numerical Treatment of the Modified Time Fractional Fokker-Planck Equation

Authors: Yuxin Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A numerical method for the modified time fractional Fokker-Planck equation is proposed. Stability and convergence of the method are
rigorously discussed by means of the Fourier method. We prove that the difference
scheme is unconditionally stable, and convergence order isO(τ+h4), where τ and h are the temporal and spatial step sizes, respectively. Finally, numerical results are given to confirm the theoretical analysis.


---

## Body

## 1. Introduction

Fractional differential equations have attracted considerable attention due to their frequent appearance in various applications in fluid mechanics, biology, physics, and engineering [1, 2]. Usually, fractional differential equations do not have analytic solutions and can only be solved by some semianalytical and numeric techniques. Recently, several semianalytical methods, such as variational iteration method, homotopy perturbation method, Adomian decomposition method, homotopy analysis method, and collocation method, have been used to solve fractional differential equations [3–7]. Meanwhile, some effective numerical techniques are developed; see [8–17].In the present paper, we are motivated to study the following modified time fractional Fokker-Planck equation [18]:
(1)∂u(x,t)∂t=[κα∂2∂x2-να∂∂x]D0,t1-αRLu(x,t)+f(x,t),hhhhhhhhhhhhhhhhhhhhhhi0<x<L,0≤t≤T,(2)u(x,0)=ϕ(x),0<x<L,u(0,t)=φ1(t),0≤t≤T,u(L,t)=φ2(t),0≤t≤T,
where κα>0 is a fractional diffusion coefficient, and να>0 is a fractional friction coefficient. D0,t1-αRLu(x,t) denotes the temporal Riemann-Liouville derivative operator defined as [2]
(3)D0,t1-αRLu(x,t)=1Γ(α)∂∂t∫0tu(x,s)(t-s)1-αds,
where Γ(·) is the gamma function.The outline of the paper is as follows. In Section2, an effective numerical method for solving the modified time fractional Fokker-Planck equation is proposed. The solvability, stability, and convergence of the numerical method are discussed in Sections 3 and 4, respectively. In Section 5, we give some numerical results demonstrating the convergence orders of the numerical method. Also a conclusion is given in Section 6.


## 2. The Construction of Numerical Method

Letxj=jh(j=0,1,…,M) and tk=kτ(k=0,1,…,N), where h=L/M and τ=T/N are the uniform spatial and temporal mesh sizes, respectively, and M, N are two positive integers.Lemma 1.
Suppose that(4)[κα∂2∂x2-να∂∂x]u(xj,t)=g(xj,t);
then a fourth-order difference scheme for the above equation is given as follows:
(5)𝒜u(xj,t)=ℒg(xj,t)+𝒪(h4),
where 𝒜 and ℒ are two difference operators and are defined by
(6)𝒜=καcosh⁡(6ναh6κα)δx2-ναμxδx,ℒ=I+h212(δx2-νακαμxδx),
in which I is an unit operator and μxδx and δx2 are average central and second central difference operators with respect to x and are defined by
(7)μxδxu(xj,t)=u(xj+1,t)-u(xj-1,t)2h,δx2u(xj,t)=u(xj+1,t)-2u(xj,t)+u(xj-1,t)h2.Proof.
In view of Taylor expansion, we can obtain(8)Ej(t)=𝒜u(xj,t)-ℒg(xj,t)=κα∑k=0∞‍1(2k)!(6ναh6κα)2k×(∂2u(xj,t)∂x2+h212∂4u(xj,t)∂x4+O(h4))-να(∂u(xj,t)∂x+h26∂3u(xj,t)∂x3+O(h4))-(g(xj,t)+h212∂2g(xj,t)∂x2+O(h4))+ναh212κα(∂g(xj,t)∂x+O(h2))=[(κα∂2u(xj,t)∂x2-να∂u(xj,t)∂x)-g(xj,t)]-ναh212κα[(κα∂3u(xj,t)∂x3-να∂2u(xj,t)∂x2)hhhhhhhi-∂g(xj,t)∂x]+h212[(κα∂4u(xj,t)∂x4-να∂3u(xj,t)∂x3)hhhhhi-∂2g(xj,t)∂x2]+O(h4).
Noting (4), we easily obtain
(9)Ej(t)=O(h4).
This completes the proof.Combing (1), (4), and (5), we obtain
(10)ℒ·∂u(xj,tk)∂t=𝒜·D0,t1-αRLu(xj,tk)+ℒ·f(xj,tk)+𝒪(h4).Using the relation of the Riemann-Liouville fractional derivativeD0,t1-αRLu(x,t) and Grünwald-Letnikov fractional derivative D0,t1-γGLu(x,t), we can approximate the Riemann-Liouville fractional derivative D0,t1-αRLu(x,t) by [2]
(11)RLD0,t1-αu(xj,tk)=1τ1-α∑m=0k‍ϖm(1-α)u(xj,tk-m)+O(τ),
where ϖm(1-α)=(-1)m(1-αm).For first-order derivative∂u(xj,tk)/∂t, we apply the following backward difference scheme:
(12)∂u(xj,tk)∂t=u(xj,tk)-u(xj,tk-1)τ+O(τ).Letuik be the numerical approximation of u(xi,tk); substituting (11) and (12) into (10) and omitting the error term 𝒪(τ+h4), we can obtain the following difference scheme for solving (1):
(13)ℒ·(ujk-ujk-1)=τα𝒜·∑m=0k‍ϖm(1-α)ujk-m+τℒ·fjk.The discrete form of above system is(14)(112-c1-c2+c3)uj+1k+(56+2c2)ujk+(112+c1-c2-c3)uj-1k=(112-c1+c2ϖ1(1-α)-c3ϖ1(1-α))uj+1k-1+(56-2c2ϖ1(1-α))ujk-1+(112+c1+c2ϖ1(1-α)+c3ϖ1(1-α))uj-1k-1-c3∑m=2k‍ϖm(1-α)(uj+1k-m-uj-1k-m)+c2∑m=2k‍ϖm(1-α)(uj+1k-m-2ujk-m+uj-1k-m)+τ(112-c1)fj+1k+56τfjk+τ(112+c1)fj-1k,
where c1=ναh/24κα, c2=τακαθ/h2, c3=τανα/2h, and θ=cosh⁡(6ναh/6κα).The initial and boundary conditions can be discretized as(15)uj0=ϕ(jh),j=0,1,…,M,u0k=φ1(kτ),uMk=φ2(kτ),k=0,1,…,N.Obviously, the local truncation error of difference scheme (14) is Rjk=O(τ+h4).


## 3. The Solvability of the Difference Scheme

Firstly, let us denote(16)U0=(ϕ(x1),ϕ(x2),…,ϕ(xM-1))T,Uk=(u1k,u2k,…,uM-1k)T,Fk=(f1k,f2k,…,fM-1k)T,hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhk=0,1,…,N.Then we can give the compact form of difference scheme (14) as follows:
(17)AUk=BUk-1+∑m=2k‍CmUk-m+DFk+G,k=1,…,N,
where
(18)A=(s2s3s1s2s3⋱⋱⋱s1s2s3s1s2),B=(p2p3p1p2p3⋱⋱⋱p1p2p3p1p2),Cm=(q2q3q1q2q3⋱⋱⋱q1q2q3q1q2),D=(r2r3r1r2r3⋱⋱⋱r1r2r3r1r2),G=(g1,0,…,0︸M-3,gM-1)T,s1=112+c1-c2-c3,s2=56+2c2,s3=112-c1-c2+c3,p1=112+c1+c2ϖ1(1-α)+c3ϖ1(1-α),p2=56-2c2ϖ1(1-α),p3=112-c1+c2ϖ1(1-α)-c3ϖ1(1-α),q1=(c2+c3)ϖm(1-α),q2=-2c2ϖm(1-α),q3=(c2-c3)ϖm(1-α),r1=τ(112+c1),r2=56τ,r3=τ(112-c1),(19)g1=-s1φ1(tk)+p1φ1(tk-1)+(c2+c3)×∑m=2k‍ϖm(1-α)φ1(tk-m)+r1f(x0,tk),gM-1=-s3φ2(tk)+p3φ2(tk-1)+(c2+c3)×∑m=2k‍ϖm(1-α)φ2(tk-m)+r3f(xM,tk).Theorem 2.
Difference equation (17) is uniquely solvable.Proof.
It is well known that the eigenvalues of the matrixA are
(20)λk=s2+2s1s3cos⁡(kπM),=56+2τακαθh2+2(112-τακαθh2)2-(ναh24κα-τανα2h)2cos⁡(kπM),
where k=1,2,…,M-1.
Note thatθ≥1 and κα>0; if
(21)(112-τακαθh2)2-(ναh24κα-τανα2h)2<0,
then we easily know that λk≠0.
If(22)(112-τακαθh2)2-(ναh24κα-τανα2h)2≥0,
then
(23)λk≥56+2τακαθh2-2|112-τακαθh2|>0.
At the moment, we obtaindet⁡(A)≠0; that is to say, the matrix A is invertible. Hence, difference equation (17) has a unique solution.


## 4. Stability and Convergence Analysis

In this section, we analyze the stability and convergence of difference scheme (17) by the Fourier method [8]. Firstly, we give the stability analysis.Lemma 3.
The coefficientsϖm(1-α)(m=0,1,…) satisfy [8] as follows:(1)
ϖ0(1-α)=1, ϖ1(1-α)=α-1, ϖm(1-γ)<0, j=1,2,…;(2)
∑m=0∞ϖm(1-γ)=0.LetUjk be the approximate solution of (14) and define
(24)ρjk=ujk-Ujk,j=1,2,…,M-1,k=0,1,…,N,ρk=(ρ1k,ρ2k,…,ρM-1k)T,k=0,1,…,N,
respectively.So, we can easily obtain the following roundoff error equation:(25)(112-c1-c2+c3)ρj+1k+(56+2c2)ρjkhhhh+(112+c1-c2-c3)ρj-1khh=(112-c1+c2ϖ1(1-α)-c3ϖ1(1-α))ρj+1k-1hhhh+(56-2c2ϖ1(1-α))ρjk-1hhhh+(112+c1+c2ϖ1(1-α)+c3ϖ1(1-α))ρj-1k-1hhhh-c3∑m=2k‍ϖm(1-α)(ρj+1k-m-ρj-1k-m)hhhh+c2∑m=2k‍ϖm(1-α)(ρj+1k-m-2ρjk-m+ρj-1k-m),hhhhhhhhhhhhhhhhhhhhhhhj=1,2,…,M-1,ρ0k=ρMk=0,k=0,1,…,N.
Now, we define the grid functions
(26)ρk(x)={ρjk,whenxj-h2<x≤xj+h2,hhhhhhhj=1,2,…,M-1,0,when0≤x≤h2orL-h2<x≤L;
then ρk(x) can be expanded in a Fourier series:
(27)ρk(x)=∑l=-∞∞‍ξk(l)exp⁡(2πlxLi),
where
(28)ξk(l)=1L∫0L‍ρk(x)exp⁡(-2πlxLi)dx.We introduce the following norm:(29)∥ρk∥2=(∑j=1M-1‍h|ρjk|2)1/2=[∫0L‍|ρk(x)|2dx]1/2,
and according to the Parseval equality
(30)∫0L‍|ρk(x)|2dx=∑l=-∞∞‍|ξk(l)|2,
we obtain
(31)∥ρk∥22=∑l=-∞∞‍|ξk(l)|2.
Through the above analysis, we can suppose that the solution of (25) has the following form:
(32)ρjk=ξkexp⁡(iβjh),
where β=2πl/L.Substituting the above expression into (25) one gets
(33)[1-(13-4c2)sin2(12βh)-i·2(c1-c3)sin(βh)]ξk=[1-(13+4c2ϖ1(1-α))sin2(12βh)-i·2(c1+c3ϖ1(1-α))sin(βh)(13+4c2ϖ1(1-α))]ξk-1-[4c2sin2(12βh)+i·2c3sin(βh)]×∑m=2k‍ϖm(1-α)ξk-m,k=0,1,…,N.Lemma 4.
The following relation holds:(34)|QS|≤1,
where
(35)Q=[1-(13+4c2ϖ1(1-α))sin2(12βh)-i·2(c1+c3ϖ1(1-α))sin(βh)(12βh)],S=[1-(13-4c2)sin2(12βh)-i·2(c1-c3)sin(βh)(12βh)].Proof.
Because of(36)8c2-32c1c3=8τακαh2cosh⁡(6ναh6κα)-2τανα23κα=8τακαh2[1+∑k=2∞‍1(2k)!(6ναh6κα)2k]>0,
we obtain
(37){8αc2[1-13sin2(12βh)]-32αc1c3cos⁡2(12βh)}+16c22(1-ϖ1(1-α)2)sin2(12βh)+16c32(1-ϖ1(1-α)2)cos⁡2(12βh)>0.
Furthermore, we can rewrite the above inequality as(38)[1-13sin2(12βh)-4c2ϖ1(1-α)sin2(12βh)]2+4(c1+c3ϖ1(1-α))2sin2(βh)≤[1-13sin2(12βh)+4c2sin2(12βh)]2+4(c1-c3)2sin2(βh);
that is,
(39)|(12βh)-1(1-(13+4c2ϖ1(1-α))sin2(12βh)-i·2(c1+c3ϖ1(1-α))hhh×sin(βh)(12βh))×(1-(13-4c2)sin2(12βh)hhhhhhhhhhhhhhhhhhh-i·2(c1-c3)sin(βh)(12βh))-1|≤1.
This completes the proof of Lemma 4.Lemma 5.
Supposing thatξk(k=1,…,N) is the solution of (33), then we have
(40)|ξk|≤exp⁡(M(k-1)τ)|ξ0|,k=1,…,N.Proof.
Fork=0, from (33), we get
(41)|ξ1|=|QS||ξ0|.
In the light of Lemma4, it is clear that
(42)|ξ1|≤|ξ0|=exp⁡(M·0τ)|ξ0|.
Now, we suppose that(43)|ξℓ|≤exp⁡(M(ℓ-1)τ)|ξ0|,(ℓ=1,2,…,k-1).
Fork>0, from (33) with Lemmas 3 and 4, we have
(44)|ξk|=1|S||∑m=2k‍Qξk-1-[4c2sin2(12βh)+i·2c3sin(βh)]×∑m=2k‍ϖm(1-α)ξk-m|=1|S||Qξk-1-[4c2sin2(12βh)+i·2c3sin(βh)]×∑m=2k-1‍ϖm(1-α)ξk-m-[4c2sin2(12βh)+i·2c3sin(βh)]ξ0|≤1|S|{|Q||ξk-1|-|4c2sin2(12βh)+i·2c3sin(βh)|×∑m=2k-1‍ϖm(1-α)|ξk-m|-|4c2sin2(12βh)+i·2c3sin(βh)|ϖk(1-α)|ξ0|}≤1|S|{(12βh)|Q|exp⁡(M(k-2)τ)|ξ0|-|4c2sin2(12βh)+i·2c3sin(βh)|×∑m=2k-1‍ϖm(1-α)exp⁡(M(k-m-1)τ)|ξ0|-|4c2sin2(12βh)+i·2c3sin(βh)|ϖk(1-α)|ξ0|}≤1|S|{|Q|-|4c2sin2(12βh)+i·2c3sin(βh)|∑m=1∞‍ϖm(1-α)}×exp⁡(M(k-2)τ)|ξ0|=1|S|{|Q|+α|4c2sin2(12βh)+i·2c3sin(βh)|}×exp⁡(M(k-2)τ)|ξ0|≤(1+Mτ)exp⁡(M(k-2)τ)|ξ0|≤exp⁡(Mτ)exp⁡(M(k-2)τ)|ξ0|=exp⁡(M(k-1)τ)|ξ0|=K|ξ0|.
This finishes the proof of Lemma 5.Lemma 6.
Difference scheme (14) is unconditionally stable.Proof.
According to Lemma5, we obtain
(45)∥ρk∥2=(∑i=1M-1‍h|ρik|2)1/2=(∑i=1M-1‍h|ξkexp⁡(iβjh)|2)1/2=(∑i=1M-1‍h|ξk|2)1/2≤K(∑i=1M-1‍h|ξ0|2)1/2=K(∑i=1M-1‍h|ξ0exp⁡(iβjh)|2)1/2=K(∑i=1M-1‍h|ρi0|2)1/2=∥ρ0∥2,k=1,2,…,N,
which means that difference scheme (14) is unconditionally stable.
Next, we give the convergence analysis. Suppose(46)ℰik=u(xi,tk)-uik,i=1,…,M-1,k=1,…,N,
and denote
(47)ℰk=(ℰ1k,ℰ2k,…,ℰM-1k)T,Rk=(R1k,R2k,…,RM-1k)T,hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhk=1,…,N.
From (14), we obtain
(48)(112-c1-c2+c3)ℰj+1k+(56+2c2)ℰjk+(112+c1-c2-c3)ℰj-1k=(112-c1+c2ϖ1(1-α)-c3ϖ1(1-α))ℰj+1k-1+(56-2c2ϖ1(1-α))ℰjk-1+(112+c1+c2ϖ1(1-α)+c3ϖ1(1-α))ℰj-1k-1-c3∑m=2k‍ϖm(1-α)(ℰj+1k-m-ℰj-1k-m)+c2∑m=2k‍ϖm(1-α)(ℰj+1k-m-2ℰjk-m+ℰj-1k-m)+τ(112-c1)Rj+1k+56τRjk+τ(112+c1)Rj-1k,j=1,2,…,M-1,k=0,1,…,N.
Similar to the stability analysis method, we define the grid functions(49)ℰk(x)={ℰjk,whenxj-h2<x≤xj+h2,hhhhhhhhj=1,2,…,M-1,0,when0≤x≤h2orL-h2<x≤L,Rk(x)={Rjk,whenxj-h2<x≤xj+h2,hhhhhhhhj=1,2,…,M-1,0,when0≤x≤h2orL-h2<x≤L;
then ℰk(x) and Rk(x) can be expanded to the following Fourier series, respectively:
(50)ℰk(x)=∑l=-∞∞‍ζk(l)exp⁡(2πlxLi),Rk(x)=∑l=-∞∞‍ηk(l)exp⁡(2πlxLi),
where
(51)ζk(l)=1L∫0L‍ℰk(x)exp⁡(-2πlxLi)dx,ηk(l)=1L∫0L‍Rk(x)exp⁡(-2πlxLi)dx.
The same as before, we also have
(52)∥ℰk∥2=(∑j=1M-1‍h|ℰjk|2)1/2=(∑l=-∞∞‍|ζk(l)|2)1/2,(53)∥Rk∥2=(∑j=1M-1‍h|Rjk|2)1/2=(∑l=-∞∞‍|ηk(l)|2)1/2.
Based on the above analysis, we can assume that ℰik and Rik have the following forms:
(54)ℰjk=ζkexp⁡(iβjh),Rjk=ηkexp⁡(iβjh),
respectively. Substituting the above two expressions into (48) yields
(55)Sζk=Qζk-1-[4c2sin2(12βh)+i·2c3sin(βh)]×∑m=2k‍ϖm(1-α)ζk-m+τ{[1-13sin2(12βh)]-i·2c1sin(βh)}ηk,iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiik=1,2,…,N.Lemma 7.
Letζk(k=0,1,…,N) be the solution of (55); then there exists a positive constant C3, so that
(56)|ζk|≤C3kτexp⁡(M(k-1)τ)|η1|,k=0,1,…,N.Proof.
Fromℰ0=0, we have
(57)ζ0=ζ0(l)=0.
In view of the convergence of the series of (53), there is a positive constant C1, such that
(58)|ηk|=|ηk(l)|≤C1|η1|=C1|η1(l)|.
Fork=1, from (55), we have
(59)ζ1=QSζ0+τ{[1-(1/3)sin2((1/2)βh)]-i·2c1sin(βh)}Sη1=τ{[1-(1/3)sin2((1/2)βh)]-i·2c1sin(βh)}Sη1.
Noticing (58), then
(60)|ζ1|≤τ|η1|≤C1τexp⁡(M·0τ)|η1|.
Now, we suppose that(61)|ζℓ|≤C1ℓτexp⁡(M(ℓ-1)τ)|η1|,ℓ=1,…,N-1.
Then when k>1, we obtain
(62)|ζk|=1|S||Qζk-1-[4c2sin2(12βh)+i·2c3sin(βh)]hhhh×∑m=2k‍ϖm(1-α)ζk-mhhhh+τ{[1-13sin2(12βh)]-i·2c1sin(βh)}ηk|≤1|S|{|Q||ζk-1|-|4c2sin2(12βh)+i·2c3sin(βh)|hhhh×∑m=2k‍ϖm(1-α)|ζk-m|hhhh+τ|[1-13sin2(12βh)]-i·2c1sin(βh)||ηk|}≤1|S|{(12βh)|Q|(k-1)exp⁡(M(k-2)τ)hhhh-|4c2sin2(12βh)+i·2c3sin(βh)|(k-m)hhhh×∑m=2k‍ϖm(1-α)exp⁡(M(k-m-1)τ)hhhh+|[1-13sin2(12βh)]-i·2c1sin(βh)|}C1τ|η1|≤1|S|{(12βh)|Q|(k-1)exp⁡(M(k-2)τ)hhhh+α|4c2sin2(12βh)+i·2c3sin(βh)|hhhh×(k-1)exp⁡(M(k-1)τ)hhhh+|[1-13sin2(12βh)]-i·2c1sin(βh)|}C1τ|η1|≤1|S|{(k-1)[|Q|+α|4c2sin2(12βh)+i·2c3sin(βh)|]hhhh+|[1-13sin2(12βh)]-i·2c1sin(βh)|}×C1τexp⁡(M(k-2)τ)|η1|≤{|[1-(1/3)sin2((1/2)βh)]-i·2c1sin(βh)||S|(k-1)[|Q|+α|4c2sin2(12βh)+i·2c3sin(βh)|]hhh+|[1-(1/3)sin2((1/2)βh)]-i·2c1sin(βh)||S|}×C1τexp⁡(M(k-2)τ)|η1|≤{k(1+Mτ)}C1τexp⁡(M(k-2)τ)|η1|≤C1kτexp⁡(Mτ)exp⁡(M(k-2)τ)|η1|=C1kτexp⁡(M(k-1)τ)|η1|.
This completes the proof.Theorem 8.
Difference scheme (14) is convergent, and the convergence order is O(τ+h4).Proof.
Firstly, we know that there are exist positive constantsC2, such that
(63)|Rjk|≤C2(τ+h4),(64)∥Rjk∥≤C2(M-1)h(τ+h4)≤C2L(τ+h4).
Using (52) and (64) with Lemma 7, we get
(65)∥ℰk∥2≤C1kτexp⁡(M(k-1)τ)∥R1∥2≤C1C2Lkτexp⁡(MT)(τ+h4)=C(τ+h4).
Due to k≤N,
(66)kτ≤T,
so,
(67)∥ℰk∥2≤C(τ+h4),
where C=C1C2TLexp⁡(MT). This ends the proof.Remark 9.
From above discussion, we know that difference scheme is an implicit scheme and it is unconditionally stable and convergent. If we take(68)∂u(xj,tk)∂t=u(xj,tk+1)-u(xj,tk)τ+O(τ)
in (1), then we can obtain an explicit scheme and it is conditionally stable and convergent.


## 5. Numerical Example

In this section, a numerical example is presented to confirm our theoretical results.Example 10.
Consider the following equation:(69)∂u(x,t)∂t=[∂2∂x2-∂∂x]D0,t1-αRLu(x,t)+f(x,t),hhhhhhhhhhhhhhhhhhh0<x<2,0≤t≤1,
with the initial and boundary conditions
(70)u(x,0)=0,0≤x≤2,u(0,t)=0,u(2,t)=0,0<t≤1,
where f(x,t)=2tsin(x)sin(2-x)+(2t1+α/Γ(2+α))[2cos⁡2(x-1)+sin2(1-x)]. The analytical solution of this equation is u(x,t)=t2sin(x)sin(2-x).
The maximum error, temporal, and spatial convergence orders by difference scheme (14) for various α are listed in Table 1. From the obtained results, we can draw the following conclusions: the experimental convergence orders are approximately 1 and 4 in temporal and spatial directions, respectively. Figures 1 and 2 show the comparison of the numerical solution with the analytical solution at t=0.5 and α=0.8 for different temporal and spatial mesh sizes.
By Table1 and Figures 1 and 2, it can be seen that the numerical solution is in excellent agreement with the analytical solution. These results confirm our theoretical analysis.Table 1
The maximum error, temporal, and spatial convergence orders by difference scheme (14) for different α.
α
The maximum error
Temporal convergence order
Spatial convergence order
0.1
h
=
1
/
4
,
τ
=
1
/
4
0.1650
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0103
1.0004
4.0017
0.2
h
=
1
/
4
,
τ
=
1
/
4
0.1526
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0095
1.0014
4.0057
0.3
h
=
1
/
4
,
τ
=
1
/
4
0.1400
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0087
1.0021
4.0083
0.4
h
=
1
/
4
,
τ
=
1
/
4
0.1274
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0080
0.9983
3.9932
0.5
h
=
1
/
4
,
τ
=
1
/
4
0.1150
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0072
0.9994
3.9975
0.6
h
=
1
/
4
,
τ
=
1
/
4
0.1027
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0065
0.9955
3.9819
0.7
h
=
1
/
4
,
τ
=
1
/
4
0.0907
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0058
0.9917
3.9670
0.8
h
=
1
/
4
,
τ
=
1
/
4
0.0789
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0051
0.9879
3.9515
0.9
h
=
1
/
4
,
τ
=
1
/
4
0.0676
—
—
h
=
1
/
8
,
τ
=
1
/
64
0.0045
0.9773
3.9090Figure 1
The comparison of the numerical solution with the analytical solution att=0.5. (τ=1/400,h=1/20).Figure 2
The comparison of the numerical solution with the analytical solution att=0.5. (τ=1/200,h=1/50).


## 6. Conclusion

In this paper, a computationally effective numerical method is proposed for simulating the modified time fractional Fokker-Planck equation. It has proven the unconditional stability and solvability of proposed scheme. Also we showed that the method is convergent with orderO(τ+h4). The numerical results demonstrate the effectiveness of the proposed scheme.



---

*Source: 282190-2014-03-27.xml*