# On the Numerical Solution of Fractional Parabolic Partial Differential Equations with the Dirichlet Condition

Authors: Allaberen Ashyralyev; Zafer Cakir
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

The first and second order of accuracy stable difference schemes for the numerical solution of the mixed problem for the fractional parabolic equation are presented. Stability and almost coercive
stability estimates for the solution of these difference schemes are obtained. A procedure of modified
Gauss elimination method is used for solving these difference schemes in the case of one-dimensional
fractional parabolic partial differential equations.


---

## Body

## 1. Introduction

It is known that various problems in fluid mechanics (dynamics, elasticity) and other areas of physics lead to fractional partial differential equations. Methods of solutions of problems for fractional differential equations have been studied extensively by many researchers (see, e.g., [1–28] and the references therein).The role played by stability inequalities (well posedness) in the study of boundary value problems for parabolic partial differential equations is well known (see, e.g., [29–34]). In the present paper, the mixed boundary value problem for the fractional parabolic equation
(1.1)∂u(t,x)∂t+Dt1/2u(t,x)-∑p=1m‍(ap(x)uxp)xp=f(t,x),x=(x1,…,xm)∈Ω,0<t<T,u(t,x)=0,x∈S,u(0,x)=0,x∈Ω-
is considered. Here Dt1/2=D0+1/2 is the standard Riemann-Liouville's derivative of order 1/2 and Ω is the open cube in the m-dimensional Euclidean space
(1.2)ℝm:{x∈Ω:x=(x1,…,xm);0<xj<1,1≤j≤m}
with boundary S,Ω-=Ω∪S,ap(x)(x∈Ω) and f(t,x)(t∈(0,T),x∈Ω) are given smooth functions and ap(x)≥a>0.Thefirst and second order of accuracy in t and second orders of accuracy in space variables difference schemes for the approximate solution of problem (1.1) are presented. The stability and almost coercive stability estimates for the solution of these difference schemes are established. A procedure of modified Gauss elimination method is used for solving these difference schemes in the case of one-dimensional fractional parabolic partial differential equations.


## 2. Difference Schemes and Stability Estimates

The discretization of problem (1.1) is carried out in two steps. In the first step, let us define the grid space
(2.1)Ω-h={x=xp=(h1p1,…,hmpm),p=(p1,…,pm),0≤pj≤Mj,hjMj=1,j=1,…,m},Ωh=Ω-h∩Ω,Sh=Ω-h∩S.
We introduce the Hilbert space L2h=L2(Ω-h) of the grid function φh(x)={φ(h1j1,…,hmjm)} defined on Ω-, equipped with the norm
(2.2)∥φh∥L2(Ω-h)=(∑x∈Ω-h‍|φh(x)|2h1⋯hm)1/2.
To the differential operator Ax generated by problem (1.1), we assign the difference operator Ahx by the formula
(2.3)Ahxuh=-∑p=1m‍(ap(x)ux-ph)xp,jp
acting in the space of grid functions uh(x), satisfying the conditions uh(x)=0 for all x∈Sh. It is known that Ahx is a self-adjoint positive definite operator in L2(Ω-h). Here,
(2.4)φxp,jp=1hp(φ(h1j1,…,hj(jj+1),…,hmjm)-φ(h1j1,…,hjjj,…,hmjm)),φx¯p,jp=1hp(φ(h1j1,…,hjjj,…,hmjm)-φ(h1j1,…,hj(jj-1),…,hmjm)).
With the help of Ahx, we arrive at the initial boundary value problem
(2.5)dvh(t,x)dt+Dt1/2vh(t,x)+Ahxvh(t,x)=fh(t,x),0<t<T,x∈Ωh,vh(0,x)=0,x∈Ω-
for a finite system of ordinary fractional differential equations.In the second step, applying the first order of approximation formula(2.6)Dtk1/2uk=1π∑r=1k‍Γ(k-r+1/2)(k-r)!(ur-ur-1τ1/2)
for
(2.7)Dtk1/2u(tk)=1Γ(1/2)∫0tk(tk-s)-1/2u′(s)ds
(see [35]) and using the first order of accuracy stable difference scheme for parabolic equations, one can present the first order of accuracy difference scheme with respect to t(2.8)ukh(x)-uk-1h(x)τ+Dtk1/2ukh(x)+Ahxukh(x)=fkh(x),x∈Ω-h,fkh(x)=fh(tk,x),tk=kτ,1≤k≤N,Nτ=T,u0h(x)=0,x∈Ω-h
for the approximate solution of problem (2.5). Here
(2.9)Γ(k-r+12)=∫0∞tk-r+1/2e-tdt.Moreover, applying the second order of approximation formula(2.10)Dtk-τ/21/2uk={-223πτu0+223πτu1+2τ3πu′(0),k=1,6πτ{45u0+25u1+25u2}-6τ5πu′(0),k=2,d∑m=2k-1‍{[(k-m)b1(k-m)+b2(k-m)]um-2+[(2m-2k-1)b1(k-m)-2b2(k-m)]um-1+[(k-m+1)b1(k-m)+b2(k-m)]um}+c[-uk-2-4uk-1+5uk],3≤k≤N
for
(2.11)Dtk-τ/21/2u(tk-τ2)=1Γ(1/2)∫0tk-τ/2‍(tk-τ2-s)-1/2u′(s)ds
(see [27]) and the Crank-Nicholson difference scheme for parabolic equations, one can present the second order of accuracy difference scheme with respect to t and to x and
(2.12)ukh(x)-uk-1h(x)τ+Dtk1/2ukh(x)+12Ahx(ukh(x)+uk-1h(x))=fkh(x),x∈Ω-h,fkh(x)=f(tk-τ2,x),tk=kτ,1≤k≤N,Nτ=T,u0h(x)=0,x∈Ω-h
for the approximate solution of problem (2.5). Here and in the future
(2.13)d=2πτ,c=26πτ,b1(r)=r+12-r-12,b2(r)=-13((r+12)3/2-(r-12)3/2).Theorem 2.1.
Letτ and |h|=h12+⋯+hn2 be sufficiently small positive numbers. Then, the solutions of difference scheme (2.8) and (2.12) satisfy the following stability estimate:
(2.14)max1≤k≤N∥ukh∥L2h≤C1max1≤k≤N∥fkh∥L2h,
where C1 does not depend on τ,h and fkh, 1≤k≤N.Proof.
We consider the difference scheme (2.8). We have that
(2.15)ukh(x)=∑s=1k‍Rk-s+1Fsh(x)τ,1≤k≤N,
where
(2.16)R=(I+τAhx)-1,Fkh(x)=fkh(x)-Dtk1/2ukh(x),Dtk1/2ukh(x)=1π∑m=1k‍Γ(k-m+1/2)(k-m)!τ-1/2[-Dtm1/2umh(x)+fmh(x)].
Using formula (2.15), we can write
(2.17)ukh(x)=∑s=1k‍Rk-s+1[-Dts1/2ush(x)+fsh(x)]τ=-∑s=1k‍Rk-s+1Dts1/2ush(x)τ+∑s=1k‍Rk-s+1fsh(x)τ,1≤k≤N.
First, we will prove that
(2.18)max1≤k≤N∥Dtk1/2ukh∥L2h≤Mmax1≤k≤N∥fkh∥L2h.
Using formula (2.17), we get
(2.19)ukh(x)-uk-1h(x)τ=-Dtk1/2ukh(x)+fkh(x)-Ahxukh(x)=-Dtk1/2ukh(x)+fkh(x)+∑s=1k‍AhxRk-s+1Dts1/2ush(x)τ-∑s=1k‍AhxRk-s+1fsh(x)τ.
Using formulas (2.16) and (2.19), we obtain
(2.20)Dtk1/2ukh(x)=1π∑m=1k‍Γ(k-m+1/2)(k-m)!(umh(x)-um-1h(x)τ1/2)=1π∑m=1k‍Γ(k-m+1/2)(k-m)!τ1/2[-Dtm1/2umh(x)+fmh(x)]+1π∑s=1k‍∑m=sk‍Γ(k-m+1/2)(k-m)!τ3/2AhxRm-s+1Dts1/2ush(x)-1π∑s=1k‍∑m=sk‍Γ(k-m+1/2)(k-m)!τ3/2AhxRm-s+1fsh(x).
Now, let us estimate zk=∥Dtk1/2ukh∥L2h,1≤k≤N. Applying the triangle inequality and the estimate [34]
(2.21)∥AhxRk∥L2h→L2h≤Mkτ,∥Rk∥L2h→L2h≤M,1≤k≤N,
we get
(2.22)zk≤1π∑m=1k‍Γ(k-m+1/2)(k-m)!τ1/2[zm+∥fmh∥L2h]+1π∑s=1k‍∥∑m=sk‍Γ(k-m+1/2)(k-m)!AhxRm-s+1∥L2h→L2hzsτ3/2+1π∑s=1k‍∥∑m=sk‍Γ(k-m+1/2)(k-m)!AhxRm-s+1∥L2h→L2h∥fsh∥L2hτ3/2≤M3∑s=1k-1‍1(k-s)ττ[zs+∥fsh∥L2h]+M4[zs+∥fsh∥L2h]τ1/2
for any k=1,…,N.  Then, using the difference analogy of integral inequality, we get (2.18).
Second, applying formula (2.17), estimates (2.18) and (2.21), we obtain
(2.23)∥ukh∥L2h=∑s=1k‍∥Rk-s+1∥L2h→L2h∥Dts1/2ush∥L2hτ+∑s=1k‍∥Rk-s+1∥L2h→L2h∥fsh∥L2hτ≤C1max1≤k≤N∥fkh∥L2h.
Estimate (2.14) for the solution of (2.8) is proved. The proof of estimate (2.14) for the solution of (2.12) follows the scheme of the proof of estimate (2.14) for the solution of (2.8) and rely on the estimate
(2.24)∥AhxBkC2∥L2h→L2h≤1kτ,∥Bk∥L2h→L2h≤1,1≤k≤N.
Here,
(2.25)B=(I-τ2Ahx)(I+τ2Ahx)-1,C=(I+τ2Ahx)-1.
Theorem 2.1 is proved.Theorem 2.2.
Letτ and |h|=h12+⋯+hn2 be sufficiently small positive numbers. Then, the solutions of difference scheme (2.8) satisfy the following almost coercive stability estimate:
(2.26)max1≤k≤N∥ukh-uk-1hτ∥L2h+max1≤k≤N∑p=1m‍∥(ukh)x-pxp,jp∥L2h≤C2ln1τ+|h|max1≤k≤N∥fkh∥L2h,
where C2 is independent of τ,h and fkh, 1≤k≤N.Proof.
We will prove the estimate(2.27)max1≤k≤N∥ukh-uk-1hτ∥L2h≤Mmin{ln1τ,1+|ln∥Ahx∥L2h→L2h|}max1≤k≤N∥fkh∥L2h.
Using formula (2.19) and estimate (2.21), we obtain
(2.28)max1≤k≤N∑s=1m‍∥AhxRk-s+1fshτ∥L2h≤Mmin{ln1τ,1+|ln∥Ahx∥L2h→L2h|}max1≤k≤N∥fkh∥L2h,max1≤k≤N∑s=1m‍∥AhxRk-s+1Dts1/2ushτ∥L2h≤Mmin{ln1τ,1+|ln∥Ahx∥L2h→L2h|}max1≤k≤N∥Dtk1/2ukhτ∥L2h
and estimate (2.18), the triangle inequality and equation (2.8), we get (2.27). From that it follows:
(2.29)max1≤k≤N∥Ahxukh∥L2h≤M1min{ln1τ,1+|ln∥Ahx∥L2h→L2h|}max1≤k≤N∥fkh∥L2h.
Then, the proof of estimate (2.26) is based on estimates (2.27), (2.29), and the following theorem on coercivity inequality for the solution of the elliptic difference problem in L2h.
Theorem  2.3.For the solutions of the elliptic difference problem(2.30)Ahxuh(x)=wh(x),x∈Ωh,uh(x)=0,x∈Shthe following coercivity inequality holds (see [14, 36])(2.31)∑p=1m‍∥uxpx-p,jph∥L2h≤C∥wh∥L2h,whereCdoes not depend onhandwh.
Theorem2.2 is proved.Theorem 2.4.
Letτ and |h|=h12+⋯+hm2 be sufficiently small positive numbers. Then, the solutions of difference scheme (2.12) satisfy the following almost coercive stability estimate:
(2.32)max1≤k≤N∥ukh-uk-1hτ∥L2h+max1≤k≤N12∑p=1m‍∥(ukh+uk-1h)x-pxp,jp∥L2h≤C3ln1τ+|h|max1≤k≤N∥fkh∥L2h,
where C3 does not depend on τ,h and fkh, 1≤k≤N.The proof of Theorem2.4 follows the proof of Theorem 2.2 and on the estimate (2.24) and the self-adjointness and positive definiteness of operator Ahx in L2h and Theorem  2.3.Remark 2.5.
The stability estimates of Theorems2.1, 2.2, and 2.4 are satisfied in the case of operator
(2.33)Au=-∑k=1n‍ak(x)∂2u∂xk2+∑k=1n‍bk(x)∂u∂xk+c(x)u
with Dirichlet condition u=0 in S. In this case, A is not self-adjoint operator in H. Nevertheless, Au=A0u+Bu and A0 is a self-adjoint positive definite operator in H and BA0-1 is bounded in H. The proof of this statement is based on the abstract results of [14] and difference analogy of integral inequality.
The method of proofs of Theorems2.1, 2.2, and 2.4 enables us to obtain the estimate of convergence of difference schemes of the first and second order of accuracy for approximate solutions of the initial-boundary value problem
(2.34)∂u(t,x)∂t-∑p=1n‍ap(x)uxpxp+∑p=1n‍bp(x)uxp+Dtαu(t,x)=f(t,x;u(t,x),ux1(t,x),…,uxn(t,x)),x=(x1,…,xn)∈Ω,0<t<T,u(0,x)=0,x∈Ω-,u(t,x)=0,x∈S
for semilinear fractional parabolic partial differential equations.Note that, one has not been able to obtain a sharp estimate for the constant figuring in the stability estimates of Theorems2.1, 2.2, and 2.4. Therefore, our interest in the present paper is studying the difference schemes (2.8) and (2.12) by numerical experiments. Applying these difference schemes, the numerical methods are proposed in the following section for solving the one-dimensional fractional parabolic partial differential equation. The method is illustrated by numerical experiments.


## 3. Numerical Results

For the numerical result, the mixed problem(3.1)∂u(t,x)∂t+Dt1/2u(t,x)-∂∂x((1+x)∂u(t,x)∂x)=f(t,x),f(t,x)=(3+16t5π+π2t(1+x))t2sinπx-πt3cosπx,0<t<1,0<x<1,u(t,0)=u(t,1)=0,0≤t≤1,u(0,x)=0,0≤x≤1
for the one-dimensional fractional parabolic partial differential equation is considered. The exact solution of problem (3.1) is
(3.2)u(t,x)=t3sinπx.
First, applying difference scheme (2.8), we obtain
(3.3)unk-unk-1τ+1π∑r=1k‍Γ(k-r+1/2)(k-r)!(unr-unr-1τ1/2)-1h[(1+xn+1)un+1k-unkh-(1+xn)unk-un-1kh]=φnk,φnk=f(tk,xn),tk=kτ,xn=nh,1≤k≤N,1≤n≤M-1,u0k=uMk=0,0≤k≤N,un0=0,0≤n≤M.
We can rewrite it in the system of equations with matrix coefficients
(3.4)AUn+1+BUn+CUn-1=Dφn,1≤n≤M-1,U0=0~,UM=0~.
Here and in the future 0~ is the (N+1)×1 zero matrix and A=anD,C=cnD,
(3.5)D=[000⋯00010⋯00001⋯00⋯⋯⋯⋯⋯⋯000⋯10000⋯01](N+1)×(N+1),B=[b1100⋯00b21b220⋯00b31b32b33⋯00⋯⋯⋯⋯⋯⋯bN,1bN,2bN,3⋯bN,N0bN+1,1bN+1,2bN+1,3⋯bN+1,NbN+1,N+1](N+1)×(N+1),φn=[φn0φn1φn2⋮φnN-1φnN](N+1)×1,Uq=[uq0uq1uq2⋮uqN-1uqN](N+1)×1,q=n±1,n,an=-1+xn+1h2,cn=-1+xnh2,b11=1,b21=-1τ-1τ,b22=1τ+1τ+2+xn+1+xnh2,b31=-Γ(1+1/2)πτ,b32=Γ(1+1/2)-Γ(1/2)πτ-1τ,b33=1τ+1τ+2+xn+1+xnh2,bij={-Γ(i-2+1/2)πτ(i-2)!,j=1,Γ(i-j+1/2)πτ(i-j)!-Γ(i-j-1+1/2)πτ(i-j-1)!,2≤j≤i-2,Γ(1+1/2)-Γ(1/2)πτ-1τ,j=i-1,1τ+1τ+2+xn+1+xnh2,j=i,0,i<j≤N+1
for i=4,5,…,N+1 and
(3.6)φnk=[3+16kτ5π+π2(kτ)(1+nh)](kτ)2sin(πnh)-π(kτ)3cos(πnh).So, we have the second-order difference equation with respect ton matrix coefficients. This type system was developed by Samarskii and Nikolaev [37]. To solve this difference equation we have applied a procedure for difference equation with respect to k matrix coefficients. Hence, we seek a solution of the matrix equation in the following form:
(3.7)Uj=αj+1Uj+1+βj+1,UM=0,j=M-1,…,2,1,
where αj(j=1,2,…,M) are (N+1)×(N+1) square matrices and βj(j=1,2,…,M) are (N+1)×1 column matrices defined by
(3.8)αj+1=-(B+Cαj)-1A,(3.9)βj+1=(B+Cαj)-1(Dφj-Cβj),j=1,2,…,M-1,
where j=1,2,…,M-1, α1 is the (N+1)×(N+1) zero matrix and β1 is the (N+1)×1 zero matrix.Second, applying difference scheme (2.12), we obtain
(3.10)unk-unk-1τ+Dtk-τ/21/2unk-12[(1+xn)un+1k-2unk+un-1kh2+un+1k-un-1k2h+(1+xn)un+1k-1-2unk-1+un-1k-1h2+un+1k-1-un-1k-12h]=φnk,φnk=f(tk-τ2,xn),tk=kτ,xn=nh,1≤k≤N,1≤n≤M-1,u0k=uMk=0,0≤k≤N,un0=0,0≤n≤M,
where
(3.11)Dtk-τ/21/2unk={-223πτun0+223πτun1+2τ3πu′(0,xn),k=1,6πτ{45un0+25un1+25un2}-6τ5πu′(0,xn),k=2,d∑m=2k-1‍{[(k-m)b1(k-m)+b2(k-m)]unm-2+[(2m-2k-1)b1(k-m)-2b2(k-m)]unm-1+[(k-m+1)b1(k-m)+b2(k-m)]unm}+c[-unk-2-4unk-1+5unk],3≤k≤N
for any n, 1≤n≤M-1. We get the system of equations in the matrix form
(3.12)AUn+1+BUn+CUn-1=Dφn,1≤n≤M-1,U0=0~,UM=0~,
where A=anF,C=cnF,
(3.13)F=[000⋯00110⋯00011⋯00⋯⋯⋯⋯⋯⋯000⋯10000⋯11](N+1)×(N+1),B=[b1100⋯00b21b220⋯00b31b32b33⋯00⋯⋯⋯⋯⋯⋯bN,1bN,2bN,3⋯bN,N0bN+1,1bN+1,2bN+1,3⋯bN+1,NbN+1,N+1](N+1)×(N+1),D=[000⋯00010⋯00001⋯00⋯⋯⋯⋯⋯⋯000⋯10000⋯01](N+1)×(N+1),φn=[φn0φn1φn2⋮φnN-1φnN](N+1)×1,Uq=[Uq0Uq1Uq2⋮UqN-1UqN](N+1)×1,q=n±1,n,an=-12(1+xnh2+12h),cn=-12(1+xnh2-12h),b11=1,b21=-223πτ-1τ+1+xnh2,b22=223πτ+1τ+1+xnh2,b31=465πτ,b32=265πτ-1τ+1+xnh2,b33=265πτ+1τ+1+xnh2,b41=d[1b1(1)+b2(1)],b42=d[-3b1(1)-2b2(1)]-c,b43=d[2b1(1)+b2(1)]-4c-1τ+1+xnh2,b44=5c+1τ+1+xnh2,b51=d[2b1(2)+b2(2)],b52=d[-5b1(2)-2b2(2)+1b1(1)+b2(1)],b53=d[3b1(2)+b2(2)-3b1(1)-2b2(1)]-c,b54=d[2b1(1)+b2(1)]-4c-1τ+1+xnh2,b55=5c+1τ+1+xnh2,bij={d[(i-3)b1(i-3)+b2(i-3)],j=1,d[(5-2i)b1(i-3)-2b2(i-3)+(i-4)b1(i-4)+b2(i-4)],j=2,d[(i-j+1)b1(i-j)+b2(i-j)+(2j-2i+1)b1(i-j-1)-2b2(i-j-1)+(i-j-2)b1(i-j-2)+b2(i-j-2)],3≤j≤i-3,d[3b1(2)+b2(2)-3b1(1)-2b2(1)]-c,j=i-2,d[2b1(1)+b2(1)]-4c-1τ+1+xnh2,j=i-1,5c+1τ+1+xnh2,j=i,0,i<j≤N+1
for i=6,7,…,N+1 and
(3.14)φnk=[3+16kτ5π+π2(kτ)(1+nh)](kτ)2sin(πnh)-π(kτ)3cos(πnh).
So, we have again the second-order difference equation with respect to n matrix coefficients. Therefore, applying the same procedure of modified Gauss elimination method (3.7) and (3.8) difference equation (3.12).Finally, we give the results of the numerical analysis. The numerical solutions are recorded for different values ofN and M and unk represents the numerical solutions of these difference schemes at (tk,xn). The error is computed by the following formula:
(3.15)EMN=max1≤k≤N,1≤n≤M-1|u(tk,xn)-unk|.
Table 1 is constructed for N=M=20,40, and 80, respectively.Table 1
Error analysis.
Method
N
=
M
=
20
N
=
M
=
40
N
=
M
=
80
Difference scheme (2.8)
0.0040
0.0020
0.0010
Difference scheme (2.12)
0.0006726
0.0001678
0.00004187Thus, by using the Crank-Nicholson difference scheme, the accuracy of solution increases faster than the first order of accuracy difference scheme.


## 4. Conclusion

In this study, the first and second order of accuracy stable difference schemes for the numerical solution of the mixed problem for the fractional parabolic equation are investigated. We have obtained stability and almost coercive stability estimates for the solution of these difference schemes. The theoretical statements for the solution of these difference schemes for one-dimensional parabolic equations are supported by numerical example in computer. We showed that the second order of accuracy difference scheme is more accurate comparing with the first order of accuracy difference scheme.



---

*Source: 696179-2012-09-04.xml*