Stability of Difference Schemes for Fractional Parabolic PDE with the Dirichlet-Neumann Conditions

Abstract

The stable difference schemes for the fractional parabolic equation with Dirichlet and Neumann boundary conditions are presented. Stability estimates and almost coercive stability estimates with ln (1/(𝜏+|ℎ|)) for the solution of these difference schemes are obtained. A procedure of modified Gauss elimination method is used for solving these difference schemes of one-dimensional fractional parabolic partial differential equations.

1. Introduction

Theory and applications, methods of solutions of problems for fractional differential equations have been studied extensively by many researchers [1–18]. In this study, initial-boundary-value problem for the fractional parabolic equation𝜕𝑢(𝑡,𝑥)𝜕𝑡+𝐷_𝑡^1/2𝑢(𝑡,𝑥)−^𝑚_𝑝=1𝑎_𝑝(𝑥)𝑢_{𝑥𝑝𝑥𝑝}()(),𝑥+𝜎𝑢𝑡,𝑥=𝑓𝑡,𝑥𝑥=₁,…,𝑥_𝑚𝑢()∣∈Ω,0<𝑡<𝑇,𝑡,𝑥_𝑆1()=0,𝜕𝑢𝑡,𝑥𝜕𝑛⃗_𝑆2=0,0≤𝑡≤𝑇,𝑆₁∪𝑆₂=𝑆=𝜕𝑢()Ω,0,𝑥=0,𝑥∈Ω, (1.1) with Dirichlet and Neumann conditions is considered. Here 𝐷_𝑡^1/2=𝐷^1/2₀₊ is the standard Riemann-Liouville’s derivative of order 1/2 and Ω is the open cube in the 𝑚-dimensional Euclidean space𝑅^𝑚∶𝑥∈Ω∶𝑥=𝑥₁,…,𝑥_𝑚;0<𝑥_𝑗<1,1≤𝑗≤𝑚, (1.2) with boundary 𝑆,Ω=Ω∪𝑆,𝑎_𝑝(𝑥)(𝑥∈Ω) and 𝑓(𝑡,𝑥)(𝑡∈(0,𝑇),𝑥∈Ω) are given smooth functions, 𝑎_𝑝(𝑥)≥𝑎>0, 𝜎>0, and 𝑛⃗ is the normal vector to Ω.

The first and second orders of accuracy stable difference schemes for the numerical solution of problem (1.1) are presented. Stability estimates and almost coercive stability estimates with ln(1/(𝜏+|ℎ|)) for the solution of these difference schemes are obtained. The method is illustrated by numerical examples.

2. The First and Second Orders of Accuracy Stable 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Ω_ℎ=𝑥=𝑥_𝑝=ℎ₁𝑝₁,…,ℎ_𝑚𝑝_𝑚,𝑝=𝑝₁,…,𝑝_𝑚,0≤𝑝_𝑗≤𝑀_𝑗,ℎ_𝑗𝑀_𝑗,Ω=1,𝑗=1,…,𝑚_ℎ=Ω_ℎ∩Ω,𝑆_ℎ=Ω_ℎ∩𝑆. (2.1) We introduce the Hilbert space 𝐿_2ℎ=𝐿₂(Ω_ℎ) of the grid function 𝜑^ℎ(𝑥)={𝜑(ℎ₁𝑝₁,…,ℎ_𝑚𝑝_𝑚)} defined on Ω, equipped with the norm𝜑^ℎ_𝐿2(Ωℎ)=_𝑥∈Ωℎ𝜑^ℎ(𝑥)²ℎ₁⋯ℎ_𝑚^1/2. (2.2) To the differential operator 𝐴^𝑥 generated by problem (1.1), we assign the difference operator 𝐴^𝑥_ℎ by the formula𝐴^𝑥_ℎ𝑢^ℎ=−^𝑚_𝑝=1𝑎_𝑝(𝑥)𝑢^ℎ_{𝑥𝑝𝑥𝑝,𝑗𝑝}+𝜎𝑢^ℎ_𝑥, (2.3) acting in the space of grid functions 𝑢^ℎ(𝑥), satisfying the conditions 𝑢^ℎ(𝑥)=0 for all 𝑥∈𝑆¹_ℎ and 𝐷^ℎ𝑢^ℎ(𝑥)=0 for all 𝑥∈𝑆²_ℎ. Here 𝐷^ℎ𝑢^ℎ(𝑥) is the approximation of 𝜕𝑢/𝜕𝑛⃗. It is known that 𝐴^𝑥_ℎ is a self-adjoint positive definite operator in 𝐿₂(Ω_ℎ). With the help of 𝐴^𝑥_ℎ we arrive at the initial-boundary-value problem𝑑𝑣^ℎ(𝑡,𝑥)𝑑𝑡+𝐷_𝑡^1/2𝑣^ℎ(𝑡,𝑥)+𝐴^𝑥_ℎ𝑣^ℎ(𝑡,𝑥)=𝑓^ℎ(𝑡,𝑥),0<𝑡<𝑇,𝑥∈Ω_ℎ,𝑣^ℎ(0,𝑥)=0,𝑥∈Ω_ℎ, (2.4) for a finite system of ordinary fractional differential equations.

In the second step, applying the first order of approximation formula𝐷_𝜏^1/2𝑢_𝑘=1√𝜋^𝑘_𝑟=1Γ(𝑘−𝑟+1/2)()!𝑢𝑘−𝑟^𝑟−𝑢_𝑟−1𝜏^1/2,1≤𝑘≤𝑁, (2.5) for𝐷_𝑡^1/2𝑢𝑡_𝑘=1Γ(1/2)^𝑡𝑘₀𝑡_𝑘−𝑠^−1/2𝑢^′(𝑠)𝑑𝑠, (2.6) (see [19]) 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 𝑡𝑢^ℎ_𝑘(𝑥)−𝑢^ℎ_𝑘−1(𝑥)𝜏+𝐷_𝜏^1/2𝑢^ℎ_𝑘+𝐴^𝑥_ℎ𝑢^ℎ_𝑘(𝑥)=𝑓^ℎ_𝑘(𝑥),𝑥∈Ω_ℎ,𝑓^ℎ_𝑘(𝑥)=𝑓^ℎ𝑡_𝑘,𝑥,𝑡_𝑘𝑢=𝑘𝜏,1≤𝑘≤𝑁,𝑁𝜏=𝑇,_ℎ0(𝑥)=0,𝑥∈Ω_ℎ, (2.7) for the approximate solution of problem (2.4). HereΓ1𝑘−𝑟+2=^∞₀𝑡^{𝑘−𝑟+1/2}𝑒^−𝑡𝑑𝑡. (2.8) Moreover, applying the second order of approximation formula𝐷_𝜏^1/2𝑢_𝑘=−𝑑√23𝑢₀+𝑑√23𝑢₁,√𝑘=1,2𝑑65𝑢₀+𝑑√65𝑢₁+𝑑√65𝑢₂,𝑑𝑘=2,_{𝑘−1𝑚=2}(𝑘−𝑚)𝑏₁(𝑘−𝑚)+𝑏₂(𝑘−𝑚)𝑢_𝑚−2+(2𝑚−2𝑘−1)𝑏₁(𝑘−𝑚)−2𝑏₂(𝑘−𝑚)𝑢_𝑚−1+(𝑘−𝑚+1)𝑏₁(𝑘−𝑚)+𝑏₂(𝑘−𝑚)𝑢_𝑚+𝑑6√2−𝑢_𝑘−2−4𝑢_𝑘−1+5𝑢_𝑘,3≤𝑘≤𝑁, (2.9) for𝐷_𝑡^1/2𝑢𝑡_𝑘−𝜏2=1Γ(1/2)^𝑡𝑘₀^−𝜏/2𝑡_𝑘−𝜏2−𝑠^−1/2𝑢^′(𝑠)𝑑𝑠 (2.10) and the Crank-Nicholson difference scheme for parabolic equations, one can present the second order of accuracy difference scheme with respect to 𝑡 and to 𝑥 and𝑢^ℎ_𝑘(𝑥)−𝑢^ℎ_𝑘−1(𝑥)𝜏+𝐷_𝜏^1/2𝑢^ℎ_𝑘(𝑥)+12𝐴^𝑥_ℎ𝑢^ℎ_𝑘(𝑥)+𝑢^ℎ_𝑘−1(𝑥)=𝑓^ℎ_𝑘(𝑥),𝑥∈Ω_ℎ,𝑓^ℎ_𝑘(𝑥)=𝑓𝑡_𝑘−𝜏2,𝑥,𝑡_𝑘𝑢=𝑘𝜏,1≤𝑘≤𝑁,𝑁𝜏=𝑇,_ℎ0(𝑥)=0,𝑥∈Ω_ℎ, (2.11) for the approximate solution of problem (2.4). Here2𝑑=√𝜋𝜏,𝑏₁(𝑟)=1𝑟+2−1𝑟−2,𝑏₂(𝑟)1=−31𝑟+2^3/2−1𝑟−2^3/2. (2.12)

Theorem 2.1. Let 𝜏 and |ℎ|=ℎ²₁+⋯+ℎ²_𝑛 be sufficiently small numbers. Then, the solutions of difference scheme (2.7) and (2.11) satisfy the following stability estimate: max_{1≤𝑘≤𝑁}𝑢^ℎ_𝑘𝐿2ℎ≤𝐶₁max_{1≤𝑘≤𝑁}𝑓^ℎ_𝑘𝐿2ℎ. (2.13) Here 𝐶₁ does not depend on 𝜏,ℎ, and 𝑓^ℎ_𝑘, 1≤𝑘≤𝑁.

Proof. For the solution of difference scheme (2.7), we have the following formulas: 𝑢^ℎ_𝑘(𝑥)=^𝑘_𝑠=1𝑅^{𝑘−𝑠+1}𝐹^ℎ_𝑠(𝑥)𝜏,1≤𝑘≤𝑁, (2.14) where 𝑅=𝐼+𝜏𝐴^𝑥_ℎ⁻¹,𝐹^ℎ_𝑘(𝑥)=𝑓^ℎ_𝑘(𝑥)−𝐷_𝜏^1/2𝑢^ℎ_𝑘(𝑥),𝐷_𝜏^1/2𝑢^ℎ_𝑘(𝑥)=1√𝜋^𝑘_𝑚=1Γ(𝑘−𝑚+1/2)()!𝜏𝑘−𝑚^−1/2−𝐷_𝜏^1/2𝑢^ℎ_𝑚(𝑥)+𝑓^ℎ_𝑚(𝑥). (2.15) The proof of (2.13) for (2.7) is based on (2.14) and estimate 𝐴^𝑥_ℎ𝑅^𝑘_{𝐿2ℎ⟶𝐿2ℎ}≤1,𝑅𝑘𝜏^{𝑘𝐿ℎ2⟶𝐿2ℎ}≤1,1≤𝑘≤𝑁, (2.16) and the triangle inequality.
In the same manner, we can obtain (2.13) for (2.11) using the inequality 𝐴^𝑥_ℎ𝐵^𝑘𝐶²_{𝐿2ℎ⟶𝐿2ℎ}≤1,𝐵𝑘𝜏^{𝑘𝐿2ℎ⟶𝐿2ℎ}≤1,1≤𝑘≤𝑁. (2.17)

Theorem 2.2. Let 𝜏 and |ℎ|=ℎ²₁+⋯+ℎ²_𝑛 be sufficiently small numbers. Then, the solutions of difference scheme (2.7) satisfy the following almost coercive stability estimate: max_{1≤𝑘≤𝑁}𝑢^ℎ_𝑘−𝑢^ℎ_𝑘−1𝜏_𝐿2ℎ+max^𝑚_{1≤𝑘≤𝑁𝑝=1}𝑢^ℎ_{𝑘𝑥𝑝𝑥𝑝,𝑗𝑝𝐿2ℎ}≤𝐶₂1ln𝜏+|ℎ|max_{1≤𝑘≤𝑁}𝑓^ℎ_𝑘𝐿2ℎ. (2.18) Here 𝐶₂ is independent of 𝜏, ℎ, and 𝑓^ℎ_𝑘, 1≤𝑘≤𝑁.

Proof. The proof of (2.18) for (2.7) is based on (2.14) and estimate (2.16) and the triangle inequality.

Theorem 2.3. Let 𝜏 and |ℎ|=ℎ²₁+⋯+ℎ²_𝑚 be sufficiently small numbers. Then, the solutions of difference scheme (2.11) satisfy the following almost coercive stability estimate: max_{1≤𝑘≤𝑁}𝑢^ℎ_𝑘−𝑢^ℎ_𝑘−1𝜏_𝐿2ℎ+max_{1≤𝑘≤𝑁}12^𝑚_𝑝=1𝑢^ℎ_𝑘+𝑢^ℎ_{𝑘−1𝑥𝑝𝑥𝑝,𝑗𝑝𝐿2ℎ}≤𝐶₃1ln𝜏+|ℎ|max_{1≤𝑘≤𝑁}𝑓^ℎ_𝑘𝐿2ℎ. (2.19) Here 𝐶₃ does not depend on 𝜏,ℎ, and 𝑓^ℎ_𝑘, 1≤𝑘≤𝑁.

Proof. The proof of (2.19) for (2.11) is based on (2.14) and estimate (2.17) and the triangle inequality.

Remark 2.4. The method of proofs of Theorems 2.1–2.3 enables us to obtain the estimate of convergence of difference schemes of the first and second orders of accuracy for approximate solutions of the initial-boundary-value problem 𝜕𝑢(𝜕𝑡,𝑥)−𝜕𝑡^𝑛_𝑝=1𝑎_𝑝(𝑥)𝑢_{𝑥𝑝𝑥𝑝}+^𝑛_𝑝=1𝑏_𝑝(𝑥)𝑢_𝑥𝑝+𝐷^𝛼_𝑡𝑢(𝑡,𝑥)=𝑓𝑡,𝑥;𝑢(𝑡,𝑥),𝑢_𝑥1(𝑡,𝑥),…,𝑢_𝑥𝑛(,𝑥𝑡,𝑥)𝑥=₁,…,𝑥_𝑛𝑢()∈Ω,0<𝑡<𝑇,0,𝑥=0,𝑥∈𝑢()∣Ω,𝑡,𝑥_𝑆1()=0,𝜕𝑢𝑡,𝑥𝜕𝑛⃗_𝑆2=0 (2.20) for semilinear fractional parabolic partial differential equations.

Note that one has not been able to obtain a sharp estimate for the constants figuring in the stability estimates of Theorems 2.1, 2.2, and 2.3. Therefore, our interest in the present paper is studying the difference schemes (2.7) and (2.11) 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 Applications

For numerical results we consider two examples.

Example 3.1. We consider the initial-boundary-value problem 𝜕𝑢(𝑡,𝑥)𝜕𝑡+𝐷_𝑡^1/2𝑢()−𝜕𝑡,𝑥()()𝜕𝑥1+𝑥𝜕𝑢𝑡,𝑥()(),𝑓()=√𝜕𝑥+𝑢𝑡,𝑥=𝑓𝑡,𝑥𝑡,𝑥3+𝑡+16𝑡5√𝜋+(1+𝑥)2𝜋²𝑡𝑡²sin²𝜋𝑥−(1+𝑥)2𝜋²𝑡³cos²𝜋𝑥−2𝜋𝑡³𝑢()sin𝜋𝑥cos𝜋𝑥,0<𝑡<1,0<𝑥<1,𝑡,0=𝑢_𝑥()𝑢()𝑡,1=0,0≤𝑡≤1,0,𝑥=0,0≤𝑥≤1, (3.1) for the one-dimensional fractional parabolic partial differential equation. The exact solution of problem (3.1) is 𝑢(𝑡,𝑥)=𝑡³sin²𝜋𝑥. (3.2)
First, applying difference scheme (2.7), we obtain 𝑢^𝑘_𝑛−𝑢_𝑛^𝑘−1𝜏+1√𝜋^𝑘_𝑟=1Γ(𝑘−𝑟+1/2)()!𝑢𝑘−𝑟^𝑟𝑛−𝑢_𝑛^𝑟−1𝜏^1/2−1ℎ1+𝑥_𝑛+1𝑢^𝑘_𝑛+1−𝑢^𝑘_𝑛ℎ−1+𝑥_𝑛𝑢^𝑘_𝑛−𝑢^𝑘_𝑛−1ℎ+𝑢^𝑘_𝑛=𝜑^𝑘_𝑛,𝜑^𝑘_𝑛=𝑓𝑡_𝑘,𝑥_𝑛,𝑡_𝑘=𝑘𝜏,1≤𝑘≤𝑁,𝑥_𝑛𝑢=𝑛ℎ,1≤𝑛≤𝑀−1,_𝑘0=0,𝑢^𝑘_𝑀−1=𝑢^𝑘_𝑀𝑢,0≤𝑘≤𝑁,_0𝑛=0,0≤𝑛≤𝑀. (3.3) We get the system of equations in the matrix form 𝐴𝑈_𝑛+1+𝐵𝑈_𝑛+𝐶𝑈_𝑛−1=𝐷𝜑_𝑛𝑈,1≤𝑛≤𝑀−1,₀=0,𝑈_𝑀−1=𝑈_𝑀, (3.4) where 𝐴=000⋯000𝑎_𝑛0⋯0000𝑎_𝑛⋯00⋯⋯⋯⋯⋯⋯000⋯𝑎_𝑛0000⋯0𝑎_{𝑛(𝑁+1)×(𝑁+1)}, (3.5)𝑏𝐵=¹¹𝑏00⋯00₂₁𝑏₂₂𝑏0⋯00₃₁𝑏₃₂𝑏₃₃𝑏⋯00⋯⋯⋯⋯⋯⋯_𝑁1𝑏_𝑁2𝑏_𝑁3⋯𝑏_𝑁𝑁0𝑏_𝑁+1,1𝑏_𝑁+1,2𝑏_𝑁+1,3⋯𝑏_𝑁+1,𝑁𝑏_{𝑁+1,𝑁+1(𝑁+1)×(𝑁+1)}, (3.6)𝐶=000⋯000𝑐_𝑛0⋯0000𝑐_𝑛⋯00⋯⋯⋯⋯⋯⋯000⋯𝑐_𝑛0000⋯0𝑐_{𝑛(𝑁+1)×(𝑁+1)}, (3.7)𝐷=000⋯00010⋯00001⋯00⋯⋯⋯⋯⋯⋯000⋯10000⋯01_{(𝑁+1)×(𝑁+1)}, (3.8)𝜑_𝑛=𝜑⁰_𝑛𝜑¹_𝑛𝜑²_𝑛⋮𝜑_𝑛^𝑁−1𝜑^𝑁_{𝑛(𝑁+1)×1},𝑈_𝑞=𝑈⁰_𝑞𝑈¹_𝑞𝑈²_𝑞⋮𝑈_𝑞^𝑁−1𝑈^𝑁_{𝑞(𝑁+1)×1},𝑞=𝑛±1,𝑛, (3.9)𝑎_𝑛=−1+𝑥_𝑛+1ℎ²,𝑐_𝑛=−1+𝑥_𝑛ℎ²,𝑏₁₁=1,𝑏₂₁1=−√𝜏−1𝜏,𝑏₂₂=1√𝜏+1𝜏+2+𝑥_𝑛+1+𝑥_𝑛ℎ²𝑏+1,₃₁Γ()=−1+1/2√𝜋𝜏,𝑏₃₂=Γ(1+1/2)−Γ(1/2)√−1𝜋𝜏𝜏,𝑏₃₃=1√𝜏+1𝜏+2+𝑥_𝑛+1+𝑥_𝑛ℎ²𝑏+1,_𝑖𝑗=−Γ(𝑖−2+1/2)√()!√𝜋𝜏𝑖−2,𝑗=1,Γ𝜋𝜏𝑖−𝑗+1/2!−Γ𝑖−𝑗𝑖−𝑗−1+1/2!√𝑖−𝑗−1,2≤𝑗≤𝑖−2,Γ1𝜋𝜏1+21−Γ2−1𝜏1,𝑗=𝑖−1,√𝜏+1𝜏+2+𝑥_𝑛+1+𝑥_𝑛ℎ²+1,𝑗=𝑖,0,𝑖<𝑗≤𝑁+1, (3.10) for 𝑖=4,5,…,𝑁+1 and 𝜑^𝑘_𝑛=√3+𝑘𝜏+16𝑘𝜏5√𝜋+(1+𝑛ℎ)2𝜋²𝑘𝜏(𝑘𝜏)²sin²()−()𝜋𝑛ℎ1+𝑛ℎ2𝜋²(𝑘𝜏)³cos²(𝜋𝑛ℎ)−2𝜋(𝑘𝜏)³sin(𝜋𝑛ℎ)cos(𝜋𝑛ℎ). (3.11)
So, we have the second-order difference equation with respect to 𝑛 matrix coefficients. This type system was developed by Samarskii and Nikolaev [20]. To solve this difference equation we have applied a procedure for difference equation with respect to 𝑘 matrix coefficients. Hence, we seek a solution of the matrix equation in the following form: 𝑈_𝑗=𝛼_𝑗+1𝑈_𝑗+1+𝛽_𝑗+1,𝑈_𝑀=𝐼−𝛼_𝑀⁻¹𝛽_𝑀,𝑗=𝑀−1,…,2,1, (3.12) where 𝛼_𝑗(𝑗=1,2,…,𝑀) are (𝑁+1)×(𝑁+1) square matrices and 𝛽_𝑗(𝑗=1,2,…,𝑀) are (𝑁+1)×1 column matrices defined by 𝛼_𝑗+1=−𝐵+𝐶𝛼_𝑗⁻¹𝛽𝐴,_𝑗+1=𝐵+𝐶𝛼_𝑗⁻¹𝐷𝜑_𝑗−𝐶𝛽_𝑗,𝑗=1,2,…,𝑀−1, (3.13) where 𝑗=1,2,…,𝑀−1,𝛼₁ is the (𝑁+1)×(𝑁+1) zero matrix and 𝛽₁ is the (𝑁+1)×1 zero matrix.
Second, applying difference scheme (2.11), we obtain the second order of accuracy difference scheme in 𝑡 and in 𝑥𝑢^𝑘_𝑛−𝑢_𝑛^𝑘−1𝜏+𝐷_𝜏^1/2𝑢^𝑘_𝑛−121+𝑥_𝑛𝑢^𝑘_𝑛+1−2𝑢^𝑘_𝑛+𝑢^𝑘_𝑛−1ℎ²+𝑢^𝑘_𝑛+1−𝑢^𝑘_𝑛−12ℎ−𝑢^𝑘_𝑛+1+𝑥_𝑛⋅𝑢^𝑘−1_𝑛+1−2𝑢_𝑛^𝑘−1+𝑢^𝑘−1_𝑛−1ℎ²+𝑢^𝑘−1_𝑛+1−𝑢^𝑘−1_𝑛−12ℎ−𝑢_𝑛^𝑘−1=𝜑^𝑘_𝑛,𝜑^𝑘_𝑛=𝑓𝑡_𝑘−𝜏2,𝑥_𝑛,𝑡_𝑘=𝑘𝜏,𝑥_𝑛𝑢=𝑛ℎ,1≤𝑘≤𝑁,1≤𝑛≤𝑀−1,_𝑘0=0,3𝑢^𝑘_𝑀−4𝑢^𝑘_𝑀−1+𝑢^𝑘_𝑀−2𝑢=0,0≤𝑘≤𝑁,_0𝑛=0,0≤𝑛≤𝑀. (3.14) Here 𝐷_𝜏^1/2𝑢^𝑘_𝑛 is defined same as in (2.9). We get the system of equations in the matrix form 𝐴𝑈_𝑛+1+𝐵𝑈_𝑛+𝐶𝑈_𝑛−1=𝐷𝜑_𝑛𝑈,1≤𝑛≤𝑀−1,₀=0,3𝑈_𝑀−4𝑈_𝑀−1+𝑈_𝑀−2=0, (3.15) where 𝑎𝐴=000⋯00_𝑛𝑎_𝑛0⋯000𝑎_𝑛𝑎_𝑛⋯00⋯⋯⋯⋯⋯⋯000⋯𝑎_𝑛0000⋯𝑎_𝑛𝑎_{𝑛(𝑁+1)×(𝑁+1)}, (3.16)𝑏𝐵=¹¹𝑏00⋯00₂₁𝑏₂₂𝑏0⋯00₃₁𝑏₃₂𝑏₃₃𝑏⋯00⋯⋯⋯⋯⋯⋯_𝑁1𝑏_𝑁2𝑏_𝑁3⋯𝑏_𝑁𝑁0𝑏_𝑁+1,1𝑏_𝑁+1,2𝑏_𝑁+1,3⋯𝑏_𝑁+1,𝑁𝑏_{𝑁+1,𝑁+1(𝑁+1)×(𝑁+1)}, (3.17)𝑐𝐶=000⋯00_𝑛𝑐_𝑛0⋯000𝑐_𝑛𝑐_𝑛⋯00⋯⋯⋯⋯⋯⋯000⋯𝑐_𝑛0000⋯𝑐_𝑛𝑐_{𝑛(𝑁+1)×(𝑁+1)}, (3.18)𝐷=000⋯00010⋯00001⋯00⋯⋯⋯⋯⋯⋯000⋯10000⋯01_{(𝑁+1)×(𝑁+1)}, (3.19)𝜑_𝑛=𝜑⁰_𝑛𝜑¹_𝑛𝜑²_𝑛⋮𝜑_𝑛^𝑁−1𝜑^𝑁_{𝑛(𝑁+1)×1},𝑈_𝑞=𝑈⁰_𝑞𝑈¹_𝑞𝑈²_𝑞⋮𝑈_𝑞^𝑁−1𝑈^𝑁_{𝑞(𝑁+1)×1},𝑞=𝑛±1,𝑛, (3.20)𝑎_𝑛1=−21+𝑥_𝑛ℎ²+12ℎ,𝑐_𝑛1=−21+𝑥_𝑛ℎ²−122ℎ,𝑑=√,𝑏𝜋𝜏₁(𝑟)=1𝑟+2−1𝑟−2,𝑏₂(𝑟)1=−31𝑟+2^3/2−1𝑟−2^3/2,𝑏₁₁=1,𝑏₂₁𝑑√=−23−1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₂₂=𝑑√23+1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₃₁=𝑑2√65,𝑏₃₂=𝑑√65−1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₃₃=𝑑√65+1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₄₁=𝑑1𝑏₁(1)+𝑏₂(1),𝑏₄₂=𝑑−3𝑏₁(1)−2𝑏₂(1−𝑑)6√2,𝑏₄₃=𝑑2𝑏₁(1)+𝑏₂(1𝑑)−46√2−1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₄₄𝑑=56√2+1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₅₁=𝑑2𝑏₁(2)+𝑏₂(2),𝑏₅₂=𝑑−5𝑏₁(2)−2𝑏₂(2)+1𝑏₁(1)+𝑏₂(1,𝑏₅₃)=𝑑3𝑏₁(2)+𝑏₂(2)−3𝑏₁(1)−2𝑏₂(1−𝑑)6√2,𝑏₅₄=𝑑2𝑏₁(1)+𝑏₂(1𝑑)−46√2−1𝜏+1+𝑥_𝑛ℎ²+12,𝑏₅₅𝑑=56√2+1𝜏+1+𝑥_𝑛ℎ²+12,𝑏_𝑖𝑗=𝑑(𝑖−3)𝑏₁(𝑖−3)+𝑏₂(𝑑)𝑏𝑖−3),𝑗=1,(5−2𝑖₁(𝑖−3)−2𝑏₂(𝑖−3)+(𝑖−4)𝑏₁(𝑖−4)+𝑏₂(𝑑𝑏𝑖−4),𝑗=2,𝑖−𝑗+1₁𝑖−𝑗+𝑏₂𝑖−𝑗+2𝑗−2𝑖+1𝑏₁𝑖−𝑗−1−2𝑏₂𝑖−𝑗−1+𝑖−𝑗−2𝑏₁𝑖−𝑗−2+𝑏₂𝑑𝑖−𝑗−2,3≤𝑗≤𝑖−3,3𝑏₁(2)+𝑏₂(2)−3𝑏₁(1)−2𝑏₂(1−𝑑)6√2𝑑,𝑗=𝑖−2,2𝑏₁(1)+𝑏₂(1𝑑)−46√2−1𝜏+1+𝑥_𝑛ℎ²+125𝑑,𝑗=𝑖−1,6√2+1𝜏+1+𝑥_𝑛ℎ²+12,𝑗=𝑖,0,𝑖<𝑗≤𝑁+1, (3.21) for 𝑖=6,7,…,𝑁+1 and 𝜑^𝑘_𝑛=√3+𝑘𝜏+16𝑘𝜏5√𝜋+(1+𝑛ℎ)2𝜋²𝑘𝜏(𝑘𝜏)²sin²()−()𝜋𝑛ℎ1+𝑛ℎ2𝜋²(𝑘𝜏)³cos²(𝜋𝑛ℎ)−2𝜋(𝑘𝜏)³sin(𝜋𝑛ℎ)cos(𝜋𝑛ℎ). (3.22) For the solution of the matrix equation (3.15), we use the same algorithm as in the first order of accuracy difference scheme, where 𝑢_𝑀=[3𝐼−4𝛼_𝑀+𝛼_𝑀−1𝛼_𝑀]⁻¹∗[(4𝐼−𝛼_𝑀−1)𝛽_𝑀−𝛽_𝑀−1].

Example 3.2. We consider the initial-boundary-value problem 𝜕𝑢(𝑡,𝑥)𝜕𝑡+𝐷_𝑡^1/2𝑢()−𝜕𝑡,𝑥()()𝜕𝑥1+𝑥𝜕𝑢𝑡,𝑥()(),𝑓()=√𝜕𝑥+𝑢𝑡,𝑥=𝑓𝑡,𝑥𝑡,𝑥3+𝑡+16𝑡5√𝜋+(1+𝑥)2𝜋²𝑡𝑡²sin²𝜋𝑥−(1+𝑥)2𝜋²𝑡³cos²𝜋𝑥−2𝜋𝑡³𝑢sin𝜋𝑥cos𝜋𝑥,0<𝑡<1,0<𝑥<1,_𝑥()()𝑢()𝑡,0=0,𝑢𝑡,1=0,0≤𝑡≤1,0,𝑥=0,0≤𝑥≤1, (3.23) for the one-dimensional fractional parabolic partial differential equation. The exact solution of problem (3.23) is 𝑢(𝑡,𝑥)=𝑡³sin²𝜋𝑥. (3.24)
First, applying difference scheme (2.7), we obtain 𝑢^𝑘_𝑛−𝑢_𝑛^𝑘−1𝜏+1√𝜋^𝑘_𝑟=1Γ(𝑘−𝑟+1/2)()!𝑢𝑘−𝑟^𝑟𝑛−𝑢_𝑛^𝑟−1𝜏^1/2−1ℎ1+𝑥_𝑛+1𝑢^𝑘_𝑛+1−𝑢^𝑘_𝑛ℎ−1+𝑥_𝑛𝑢^𝑘_𝑛−𝑢^𝑘_𝑛−1ℎ+𝑢^𝑘_𝑛=𝜑^𝑘_𝑛,𝜑^𝑘_𝑛=𝑓𝑡_𝑘,𝑥_𝑛,𝑡_𝑘=𝑘𝜏,1≤𝑘≤𝑁,𝑥_𝑛𝑢=𝑛ℎ,1≤𝑛≤𝑀−1,_𝑘0=𝑢^𝑘₁,𝑢^𝑘_𝑀𝑢=0,0≤𝑘≤𝑁,_0𝑛=0,0≤𝑛≤𝑀. (3.25) We get the system of equations in the matrix form 𝐴𝑈_𝑛+1+𝐵𝑈_𝑛+𝐶𝑈_𝑛−1=𝐷𝜑_𝑛𝑈,1≤𝑛≤𝑀−1,₀=𝑈₁,𝑈_𝑀=0, (3.26) where matrices 𝐴,𝐵,𝐶,𝐷,𝜑_𝑛,𝑈_𝑞(𝑞=𝑛∓1,𝑛) are defined same as in (3.5), (3.6), (3.7), (3.8), (3.9), respectively.
So, we have the second-order difference equation with respect to 𝑛 matrix coefficients. To solve this difference equation we have applied a procedure for difference equation with respect to 𝑘 matrix coefficients. Hence, we seek a solution of the matrix equation in the following form: 𝑈_𝑗=𝛼_𝑗+1𝑈_𝑗+1+𝛽_𝑗+1,𝑈_𝑀=0,𝑗=𝑀−1,…,2,1, (3.27) where 𝛼_𝑗(𝑗=1,2,…,𝑀) are (𝑁+1)×(𝑁+1) square matrices and 𝛽_𝑗(𝑗=1,2,…,𝑀) are (𝑁+1)×1 column matrices defined by 𝛼_𝑗+1=−𝐵+𝐶𝛼_𝑗⁻¹𝛽𝐴,_𝑗+1=𝐵+𝐶𝛼_𝑗⁻¹𝐷𝜑_𝑗−𝐶𝛽_𝑗,𝑗=1,2,…,𝑀−1, (3.28) where 𝑗=1,2,…,𝑀−1,𝛼₁ is the (𝑁+1)×(𝑁+1) identity matrix and 𝛽₁ is the (𝑁+1)×1 zero matrix.
Second, applying the formulas 𝑢_𝑥𝑡_𝑘=𝑢,0^𝑘1−𝑢^𝑘₀ℎ−ℎ2𝑢_𝑥𝑥𝑡_𝑘,0+𝑜ℎ²𝑢,0≤𝑘≤𝑁,_𝑥𝑡_𝑘,𝑀=3𝑢^𝑘_𝑀−4𝑢^𝑘_𝑀−1+𝑢^𝑘_𝑀−2ℎ2ℎ+𝑜²𝑢,0≤𝑘≤𝑁,_𝑡𝑡_𝑘=𝑢,0^0𝑘+1−𝑢₀^𝑘−1𝜏2𝜏+𝑜²𝑢,1≤𝑘≤𝑁−1,_𝑡𝑡_𝑁,0=3𝑢^𝑁₀−4𝑢₀^𝑁−1+𝑢₀^𝑁−2𝜏2𝜏+𝑜²,𝑘=𝑁, (3.29) and applying difference scheme (2.11), we obtain the second order of accuracy difference scheme in 𝑡 and in 𝑥𝑢^𝑘_𝑛−𝑢_𝑛^𝑘−1𝜏+𝐷_𝜏^1/2𝑢^𝑘_𝑛−121+𝑥_𝑛𝑢^𝑘_𝑛+1−2𝑢^𝑘_𝑛+𝑢^𝑘_𝑛−1ℎ²+𝑢^𝑘_𝑛+1−𝑢^𝑘_𝑛−12ℎ−𝑢^𝑘_𝑛+1+𝑥_𝑛𝑢^𝑘−1_𝑛+1−2𝑢_𝑛^𝑘−1+𝑢^𝑘−1_𝑛−1ℎ²+𝑢^𝑘−1_𝑛+1−𝑢^𝑘−1_𝑛−12ℎ−𝑢_𝑛^𝑘−1=𝜑^𝑘_𝑛,𝜑^𝑘_𝑛=𝑓𝑡_𝑘−𝜏2,𝑥_𝑛,𝑡_𝑘=𝑘𝜏,𝑥_𝑛𝑢=𝑛ℎ,1≤𝑘≤𝑁,1≤𝑛≤𝑀−1,₀₀−ℎ=0,𝑘=0,𝑢4𝜏^0𝑘−1+1ℎ+ℎ2𝐷_𝑡^1/2+ℎ2𝑢^𝑘₀+ℎ𝑢4𝜏^0𝑘+1=1ℎ𝑢^𝑘₁+ℎ2𝜑^𝑘₀ℎ,1≤𝑘≤𝑁−1,𝑢4𝜏^0𝑁−2−ℎ𝜏𝑢₀^𝑁−1+1ℎ+3ℎ+ℎ4𝜏2𝐷_𝑡^1/2+ℎ2𝑢^𝑁₀=1ℎ𝑢^𝑁₁+ℎ2𝜑^𝑁₀𝑢,𝑘=𝑁,_𝑘𝑀𝑢=0,0≤𝑘≤𝑁,_0𝑛=0,0≤𝑛≤𝑀. (3.30) Here 𝐷_𝜏^1/2𝑢^𝑘_𝑛 is defined similar to (2.9). We get the system of equations in the matrix form 𝐴𝑈_𝑛+1+𝐵𝑈_𝑛+𝐶𝑈_𝑛−1=𝐷𝜑_𝑛,1≤𝑛≤𝑀−1,𝐸𝑈₀=𝐹𝑈₁+𝑅𝜑₀,𝑈_𝑀=0, (3.31) where matrices 𝐴,𝐵,𝐶,𝐷,𝜑_𝑛,𝑈_𝑞(𝑞=𝑛∓1,𝑛) are defined same as in (3.16), (3.17), (3.18), (3.19), (3.20), respectively.
For the solution of the matrix equation (3.31), we use the same algorithm as in the first order of accuracy difference scheme, where 𝑢_𝑀=0,𝛼₁=𝐸⁻¹𝐹,𝛽₁=𝐸⁻¹𝑅𝜑₀,01𝐹=000⋯00ℎ10⋯0000ℎ1⋯00⋯⋯⋯⋯⋯⋯000⋯ℎ01000⋯0ℎ_{(𝑁+1)×(𝑁+1)},𝑅=00⋯001⋯0⋯⋯⋯⋯00⋯1_{(𝑁+1)×(𝑁+1)},𝑒𝐸=¹¹𝑒00⋯00₂₁𝑒₂₂𝑒0⋯00₃₁𝑒₃₂𝑒₃₃𝑒⋯00⋯⋯⋯⋯⋯⋯_𝑁1𝑒_𝑁2𝑒_𝑁3⋯𝑒_𝑁𝑁0𝑒_𝑁+1,1𝑒_𝑁+1,2𝑒_𝑁+1,3⋯𝑒_𝑁+1,𝑁𝑒_{𝑁+1,𝑁+1(𝑁+1)×(𝑁+1)},𝑒₁₁=1,𝑒₂₁ℎ=−−4𝜏4ℎ3√𝜋𝜏,𝑒₂₂=1ℎ+ℎ2+4ℎ3√𝜋𝜏,𝑒₂₃=ℎ,𝑒4𝜏₃₁=2√2ℎ15√𝜋𝜏,𝑒₃₂=−16√2ℎ15√−ℎ𝜋𝜏4𝜏,𝑒₃₃=1ℎ+ℎ2+14√2ℎ15√𝜋𝜏,𝑒₃₄=ℎ,𝑒4𝜏₄₁=𝑑ℎ211+2𝑏₁(1)+𝑏₂(1,𝑒₄₂)=𝑑ℎ2−4𝑏₁(1)−2𝑏₂(1)+12𝑏₁(0)+𝑏₂(0,𝑒₄₃)=𝑑ℎ212+2𝑏₁(1)+𝑏₂(1)−1−2−23−ℎ,𝑒4𝜏₄₄=1ℎ+ℎ2+𝑑ℎ211+2𝑏₁(0)+𝑏₂(0),𝑒₄₅=ℎ,𝑒4𝜏₅₁=𝑑ℎ212+2𝑏₁(2)+𝑏₂(2,𝑒₅₂)=𝑑ℎ2−23𝑏₁(2)−2𝑏₂(2)+11+2𝑏₁(1)+𝑏₂(1,𝑒₅₃)=𝑑ℎ212+1+2𝑏₁(2)+𝑏₂(2)−22𝑏₁(1)−2𝑏₂(1)+12𝑏₁(0)+𝑏₂(0,𝑒₅₄)=𝑑ℎ211+1+2𝑏₁(1)+𝑏₂(1)−2𝑏₁(0)−2𝑏₂(0−ℎ),𝑒4𝜏₅₅=1ℎ+ℎ2+𝑑ℎ211+2𝑏₁(0)+𝑏₂(0),𝑒₅₆=ℎ,𝑒4𝜏_𝑖𝑗=𝑑ℎ21𝑖−3+2𝑏₁(𝑖−3)+𝑏₂(𝑖−3),𝑗=1,𝑑ℎ2−2(𝑖−2)𝑏₁(𝑖−3)−2𝑏₂()+1𝑖−3𝑖−4+2𝑏₁(𝑖−4)+𝑏₂(𝑖−4),𝑗=2,𝑑ℎ21𝑖−𝑗+1+2𝑏₁𝑖−𝑗+𝑏₂𝑖−𝑗−2𝑖−𝑗𝑏₁𝑖−𝑗−1−2𝑏₂+1𝑖−𝑗−1𝑖−𝑗−2+2𝑏₁𝑖−𝑗−2+𝑏₂−ℎ𝑖−𝑗−2,3≤𝑗≤𝑖−2,+4𝜏𝑑ℎ212+2𝑏₁(1)+𝑏₂(1)−2𝑏₁(0)−2𝑏₂(01),𝑗=𝑖−1,ℎ+ℎ2+𝑑ℎ211+2𝑏₁(0)+𝑏₂(0ℎ),𝑗=𝑖,ℎ4𝜏,𝑗=𝑖+1,+4𝜏𝑑ℎ21𝑖−𝑁+2+2𝑏₁(𝑖−𝑁+1)+𝑏₂()()𝑏𝑖−𝑁+1−2𝑖−𝑁+1₁(𝑖−𝑁)−2𝑏₂()+1𝑖−𝑁𝑖−𝑁−1+2𝑏₁(𝑖−𝑁−1)+𝑏₂(−ℎ𝑖−𝑁−1),𝑗=𝑁−1,𝜏+𝑑ℎ212+2𝑏₁(1)+𝑏₂(1)−2𝑏₁(0)−2𝑏₂(01),𝑗=𝑁,ℎ+ℎ2+3ℎ+4𝜏𝑑ℎ211+2𝑏₁(0)+𝑏₂(0),𝑗=𝑁+1,0,𝑗>𝑖+1, (3.32) for 𝑖=6,7,…,𝑁+1 and 𝜑^𝑘₀=−2𝜋²(𝑘𝜏)³. (3.33)

3.1. Error Analysis

Finally, we give the results of the numerical analysis. The error is computed by the following formula:𝐸^𝑁_𝑀=max_{1≤𝑘≤𝑁,1≤𝑛≤𝑀−1}𝑢𝑡_𝑘,𝑥_𝑛−𝑢^𝑘_𝑛, (3.34) where 𝑢(𝑡_𝑘,𝑥_𝑛) represents the exact solution and 𝑢^𝑘_𝑛 represents the numerical solutions of these difference schemes at (𝑡_𝑘,𝑥_𝑛). The numerical solutions are recorded for different values of 𝑁 and 𝑀. Tables 1 and 2 are constructed for 𝑁=𝑀=25,50, and 100, respectively.

Table 1 

Error analysis of first and second order of accuracy difference schemes for Example 3.1.

Table 2 

Error analysis of first and second order of accuracy difference schemes for Example 3.2.

The results in Tables 1 and 2 show that, by using the Crank-Nicholson difference scheme, more accurate approximate results can be obtained rather than the first order of accuracy difference scheme.

Acknowledgment

The author is grateful to Professor Allaberen Ashyralyev (Fatih University, Turkey) for his comments and suggestions to improve the quality of the paper.