A procedure of modified Gauss elimination method is applied for the solution of the first and second order of accuracy difference schemes of one-dimensional fractional parabolic differential equations

Electronic Journal of Differential Equations, Vol. 2014 (2014), No. 97, pp. 1–17.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu

WELL-POSEDNESS OF FRACTIONAL PARABOLIC DIFFERENTIAL AND DIFFERENCE EQUATIONS

WITH DIRICHLET-NEUMANN CONDITIONS

ALLABEREN ASHYRALYEV, NAZAR EMIROV, ZAFER CAKIR

Abstract. We study initial-boundary value problems for fractional parabolic equations with the Dirichlet-Neumann conditions. We obtain a stable dif- ference schemes for this problem, and obtain theorems on coercive stability estimates for the solution of the first order of accuracy difference scheme. A procedure of modified Gauss elimination method is applied for the solution of the first and second order of accuracy difference schemes of one-dimensional fractional parabolic differential equations.

1. Introduction

Theory, applications and methods of solutions of problems for fractional dif- ferential equations have been studied extensively by many researchers (see, e.g., [1]–[8], [10]–[16], [18], [19], [21], [23]–[30], [32]–[34], [39]–[46] and the references given therein). In this article, we study the initial-boundary value problem

D^α_tu(t, x)−a(x)u_xx(t, x) +σu(t, x) =f(t, x), 0< x < l, 0< t < T, u(t,0) = 0, ux(t, l) = 0, 0≤t≤T,

u(0, x) = 0, 0≤x≤l

(1.1) for the fractional parabolic equation with the Dirichlet-Neumann conditions. Here D_t^α=D^α₀₊ is the standard Riemann-Louville’s derivative of orderα∈[0,1). Here a(x)(x∈(0, l)) andf(t, x)(t∈(0, T), x∈(0, l)) are given smooth functions,a(x)≥ a > 0, σ > 0. Theorem on coercive stability estimates for the solution of the initial-boundary value problem (1.1) is established. Stable difference schemes for the approximate solution of problem (1.1) are considered. Theorem on coercive stability estimates for the solution of the first order of accuracy in t difference scheme is proved. A procedure of modified Gauss elimination method is applied for the solution of the first and second order of accuracy difference schemes for the fractional parabolic equations.

The organization of the present paper as follows. The first section is introduction where we provide the history and formulation of the problem. In Section 2, theorem on coercivity stability of problem (1.1) is established. In Section 3, stable difference

2000Mathematics Subject Classification. 35R11, 35B35, 47B39, 47B48.

Key words and phrases. Fractional parabolic equations; Dirichlet-Neumann conditions;

positive operator; difference schemes; stability.

c

2014 Texas State University - San Marcos.

Submitted December 26, 2013. Published April 10, 2014.

1

schemes for the approximate solution of problem (1.1) are considered. Theorem on coercivity stability for the first order of accuracy int difference scheme is proved.

In Section 4, the numerical application is given. Finally, Section 5 is conclusion.

2. Theorems on coercive stability We will give some statements which will be useful in the sequel.

Let E be a Banach space, and A : D(A) ⊂ E → E be a linear unbounded operator densely defined in E. We callA strongly positive in the Banach space E, if its spectrum σA lies in the interior of the sector of angle φ, 0 < 2φ < π, symmetric with respect to the real axis, and if on the edges of this sector,S1(φ) = {ρe^iφ : 0≤ρ≤ ∞} and S2(φ) ={ρe^−iφ : 0≤ρ≤ ∞}, and outside of the sector the resolvent (λ−A)⁻¹is subject to the bound

k(λ−A)⁻¹k_E→E≤ M

1 +|λ|. (2.1)

The infimum of such angles is called spectral angleϕ(A, E) ofA.

Throughout this article, positive constants have different values in time and they will be indicated withM On the other handM(α, β,· · ·) is used to focus on the fact that the constant depends only onα, β,· · ·.

For a positive operatorAin the Banach spaceE, let us introduce the fractional spacesEβ=Eβ(E, A)(0< β <1) consisting of thoseν ∈E for which the norm

kνkE_β = sup

λ>0

λ^βkA(λ+A)⁻¹νkE+kνkE

is finite.

Theorem 2.1([17, 31]). Let AandB be two commutative positive operators with ϕ(A, E) +ϕ(B, E) < π. Then it follows that there exists the bounded operator (A+B)⁻¹ defined on whole spaceE. Moreover, for every β ∈(0,1) and f, there exists a unique solution u=u(f)of the problem

Au+Bu=f and the following estimates hold

kAuk_Eβ(E,B)+kBuk_Eβ(E,B)+kBuk_Eβ(E,A)≤M(β)kfk_Eβ(E,B), kAuk_Eβ(E,A)+kBuk_Eβ(E,A)+kAuk_Eβ(E,B)≤M(β)kfk_Eβ(E,A).

Theorem 2.2 ([31]). Let A be the positive operator with ϕ(A, E)< π. Then for β≤ ¹₂, A^β is a positive operator withϕ(A^β, E)< ^π₂.

Theorem 2.3 ([3]). Let A be the operator acting in E = C[0, T] defined by the formula Av(t) = v⁰(t), with the domainD(A) = {v(t) :v⁰(t)∈C[0, T], v(0) = 0}.

ThenA is a positive operator in the Banach spaceE=C[0, T] and A^βf(t) =D^β_tf(t)

for allf(t)∈D(A).

From the above theorems it follows the following theorem.

Theorem 2.4. Let A and B be the positive operators with ϕ(A, E) < π and ϕ(B, E) ≤ ^π₂. Then for β ≤ ¹₂ it follows that there exists bounded (D^β+B)⁻¹

defined on whole space E. Moreover, for every f, there exists a unique solution u=u(f)of the problem

D^βu+Bu=f and the following estimate holds

kD^βukE_β(E,B)+kBukE_β(E,B)≤M(β)kfkE_β(E,B). Now, we consider the second order differential operator

B^xu(x) =−a(x)uxx(x) +σu(x) (2.2) with the domainD(B^x) ={u;u, u⁰, u⁰⁰∈C[0, l], u(0) = 0, u⁰(l) = 0}.

Let us introduce the Banach spaceC^γ[0, l], γ∈(0,1] of all continuous function ϕ(x) defined on [0, l] and satisfying a H¨older condition for which the following norm is finite

kϕk_Cγ[0,l]=kϕk_C[0,l]+ sup

x₁6=x₂

|ϕ(x1)−ϕ(x2)|

|x1−x₂|^γ ,

where C[0, l] is the Banach space of all continuous functionsϕ(x) defined on [0, l]

with the norm

kϕk_C[0,l] = max

x∈[0,l]|ϕ(x)|.

The positivity of the operatorB^xin the Banach spaceC[0, l] was established (see, [37, 38]). Moreover, we have that for any β ∈ (0,1/2) the norms in the spaces E_β(E, B) andC^2β[0, l] are equivalent.

Theorem 2.5. For β ∈ (0,1/2), the norms of the space Eβ(C[0, l], B^x) and the H¨older space C^2β[0, l] are equivalent.

The proof of Theorem 2.5 is based on the following estimates

|G^x(x, x₀;λ)| ≤ M(σ, a)

√σ+λ (e⁻¹²

√σ+λ

a (x−s), 0≤x0≤x, e⁻¹²

√σ+λ

a (x₀−x), x≤x0≤l,

|G^x_x(x, x₀;λ)| ≤M(σ, a) (e⁻¹²

√σ+λ

a (x−x0), 0≤x0≤x, e⁻¹²

√σ+λ

a (x₀−x), x≤x0≤l

for the Green’s function of the differential operatorB^xdefined by the formula (2.2) and it follows the scheme of the proof of the Theorem of paper [9].

Theorem 2.6. For the solution of problem (1.1)the coercive stability estimate

0≤t≤Tmax kuxx(t, .)k_Cβ[0,l]≤M(β)kf(t, .)k_Cβ[0,l]

holds, whereM(β)does not depend onf(t, x) (0≤t≤T,x∈[0, l])and0< β <1.

The proof of Theorem 2.6 is based on the positivity of differential operatorB^x defined by formula (2.2), on the Theorem 2.3 on connection of fractional derivatives with fractional powers of positive operators, on the Theorem 2.2 on spectral angle of fractional powers of positive operators, and on the Theorem 2.4 on fractional powers of coercively positive sums two operators.

3. 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, l]h= (xn=nh, 0≤n≤M, M h=l)

To the differential space operator B^x generated by formula (2.2), we assign the difference operatorB_h^x by the formula

B_h^xu^h=−a(x)u^h_x

n¯x_n+σu(x)^h (3.1)

acting in the space of grid functionsu^h(x), satisfying the conditionsu^h(x) = 0 for allx= 0 and D^hu^h(x) = 0 for x=l. Here D^hu^h(x) is the approximation of ux. With the help ofB_h^xwe arrive at the initial boundary value problem

D^α_tv^h(t, x) +B_h^xv^h(t, x) =f^h(t, x), 0< t < T, x∈[0, l]h, v^h(0, x) = 0, x∈[0, l]h

(3.2) for a finite system of ordinary fractional differential equations.

In the second step, applying the first order of approximation formula (see [3]) D^α_τuk = 1

Γ(1−α)

k

X

r=1

Γ(k−r−α+ 1) (k−r)!

ur−ur−1

τ^α ,1≤k≤N for

D^α_τu(t_k) = 1 Γ(1−α)

Z tk

0

(t_k−s)^−αu⁰(s)ds

and using the first order of accuracy stable difference scheme for parabolic equa- tions, one can present the first order of accuracy difference scheme with respect to t,

1 Γ(1−α)

k

X

r=1

Γ(k−r−α+ 1) (k−r)!

u^h_r(x)−u^h_r−1(x)

τ^α +B_h^xu^h_k(x) =f_k^h(x), f_k^h(x) =f^h(tk, x), tk=kτ, 1≤k≤N, N τ =T, x∈[0, l]h,

u^h₀(x) = 0, x∈[0, l]h

(3.3)

for the approximate solution of problem (1.1). Moreover, applying the second order of approximation formula: fork= 1,

D_τ^αuk=−d 2^α−1

(2−α)(1−α)u0+d 2^α−1

(2−α)(1−α)u1, fork= 2,

D^α_τuk =dh3^5−α 2^4−α

1

(1−α)(2−α)(3−α)−73^2−α 2^3−α

1 (1−α)(2−α)

i u0

+dh

−3^4−α 2^2−α

1

(1−α)(2−α)(3−α)+3^2−α 2^−α

1 (1−α)(2−α)

i u₁ +dh3^4−α

2^4−α

1

(1−α)(2−α)(3−α)−3^2−α 2^3−α

1 (1−α)(2−α)

iu₂,

for 3≤k≤N, D_τ^αuk=d

k−1

X

m=2

nh(k−m)

1−α ξ(k−m)−η(k−m) 2−α

i u_m−2 +h(2m−2k−1)

1−α ξ(k−m) +2η(k−m) 2−α

u_m−1 +(k−m+ 1)

1−α ξ(k−m)−η(k−m) 2−α

um

o

+dh

−2^α−2

2−αuk−2−2^α−1

1−α− 2^α−1 2−α

uk−1+2^α−1

1−α− 2^α−2 2−α

uk

i (3.4)

for

D_t^αu(t_k−τ /2) = 1 Γ(1−α)

Z tk−τ /2 0

(t_k−τ /2−s)^−αu⁰(s)ds,

and using a Crank-Nicholson difference scheme for parabolic equations, one can present the second order of accuracy difference scheme with respect totandx,

D_τ^αu^h_k(x) +1 2B_h^x

u^h_k(x) +u^h_k−1(x)

=f_k^h(x), x∈[0, l]h, f_k^h(x) =f(t_k−τ /2, x), t_k=kτ, 1≤k≤N, N τ =T,

u^h₀(x) = 0, x∈[0, l]h

(3.5)

for the approximate solution of problem (1.1). Here, d= τ^−α

Γ(1−α),ξ(r) = r+ 1/21−α

− r−1/21−α

, η(r) = r+ 1/22−α

− r−1/22−α

. Now, we consider the equation

B_h^xu^h+λu^h=f^h (3.6)

in the casea(x) = 1.

Lemma 3.1. Let λ >0. Then (3.6)is uniquely solvable, and the formula u^h= (B^x_h+λ)⁻¹f^h=n^MX⁻¹

i=1

G(k, i;λ+σ)fihoM 0

(3.7) is valid, where

G(k, i;λ+σ) = h(R^M−i−R^M+i)(R^M−k−R^M+k)

(1−R²)(1 +R^2M−1) +h(R^|k−i|+1−R^k+i+1) (1−R²) for1≤i≤M −1, and1≤k≤M,

R= (1 +δh)⁻¹, δ= 1

2 h(λ+σ) +p

(λ+σ)(4 +h²(λ+σ)) .

The grid function G(k, i;λ+σ) is called the Green function of equation (3.6) and by the formulas forR andδ, we get

M−1

X

i=1

G(k, i;λ+σ)h= 1

λ+σ− 1 λ+σ

R^k+R^2M−k−1

1 +R^2M−1 , 1≤k≤M. (3.8) To prove the positivity on B_h^x in the Banach space C_h, we need the following auxiliary lemmas [13].

Lemma 3.2. The following estimate holds

|δ| ≥maxn|λ+σ|h 2 ,p

|λ+σ|o

. (3.9)

Lemma 3.3. The following estimate

|R| ≤ 1

1 +p

|λ+σ|hcosθ <1 (3.10) is valid, where |θ|< π/2.

Theorem 3.4. For all λ in the sector Σ_θ = {λ : |argλ| ≤ θ,0 ≤ θ < π/2} the resolvent(λI+B_h^x)⁻¹ defined by (3.7)satisfies the following estimate

k(λI+B_h^x)⁻¹k_Ch→Ch≤ M(µ, θ, σ)

1 +|λ| . (3.11)

Proof. First, we consider the operator B^x_h defined by formula (3.1) in the case a(x) = 1. Let us set k=M. Since

uM = h²R(1−R^M−1)(1 +R^M−1) (1−R)(1 +R^2M−1) f_M−1

+ 1

(1−R)(1 +R^2M−1)

M−2

X

i=1

R^M−i−R^M+i h²f_i, we have that

uM

≤2

R 1−R

h²|fM−1|+ 1 1− |R|

M−2

X

i=1

|R|^M−i+|R|^M+i h²

fi

≤2h²kf^hkC_h

n| R

1−R|+ |R²| 1− |R|²

o .

Now, let us 1≤k≤M −1. Then by formula (3.7) and the triangle inequality, we obtain

|uk| ≤ |R|^M−k+|R|^M+k

|1−R²|

1 +R^2M−1

M−1

X

i=1

|R|^M−i+|R|^M+i h²

f_i

+ 1

|1−R²|

M−1

X

i=1

|R|^|k−i|+1+|R|^k+i+1 h²

fi

≤ 2

|1−R²|

M−1

X

i=1

|R|^M−i+1+|R|^M+i+1 h²

f_i

+ 1

|1−R²|

M−1

X

i=1

|R|^|k−i|+1+|R|^k+i+1 h²

fi

≤ 4h²

|1−R²|kf^hk_Ch

M−1

X

i=1

|R|^M−i+1

+ 2h²

|1−R²|kf^hkC_h

n^k−1X

i=1

|R|^k−i+1+|R|+

M−1

X

i=k+1

|R|^i−k+1o

≤ 2h²

|1−R²|kf^hkC_h

n 2|R|²

1− |R| + 2|R|²

1− |R|+|R|o

≤Mn |R|² 1− |R|²

h² 1 +R

+

R 1−R

1 1 +R

h²o

. From estimate (3.10) it follows that

|R|²

1− |R|² ≤

1 1+√

|λ+σ|hcosθ

1− ¹

1+√

|λ+σ|hcosθ

²

= 1

p|λ+σ|hcosθ ²

. (3.12) Clearly, we have that

|λ+σ|=|ρcosθ+iρsinθ+σ| =p

ρ²+ 2ρσcosθ+σ²

≥p

ρ²cos²θ+ 2ρσcosθ+σ²=|λ|cosθ+σ.

Thus 1

|λ+σ| ≤ 1

|λ|cosθ+σ ≤ 1

|λ|cosθ+σcosθ

=

1 cosθ

|λ|+σ =

1 σcosθ

1 + _σ¹|λ|

≤M(σ, θ) 1 +|λ|.

(3.13)

Combining estimates (3.12) and (3.13), we obtain that h²|R|²

1− |R|2 ≤

1 cos²θ

|λ+σ| ≤ M(σ, θ)

1 +|λ|. (3.14)

From the definition ofR and estimate (3.9), it follows that

R 1−R

h²= h

|δ| ≤ 2

1 +|λ|. (3.15)

Combining estimates (3.14) and (3.15), we obtain ku^hkC_h ≤ M(µ, σ, θ)

1 +|λ| kf^hkC_h.

This concludes the proof of Theorem 3.4 in the casea(x) = 1. Second, noted that the proof of this statement is based on estimates for the Green’s function. Under one more assumption thatσ >0 is sufficiently large number, applying a fixed point Theorem, same estimates for the Green’s function can be obtained. Therefore, this statement of theorem is true also for difference operator B^x_h defined by formula

(3.1). Theorem 3.4 is proved.

Theorem 3.5. Let 0 < β < ¹₂. Then, the norms of spaces Eβ(Ch, B_h^x) and C_h^2β are equivalent uniformly inh,0< h < h0.

Proof. From (3.7) and (3.8) it follows that λ^βB_h^x(B^x_h+λ)⁻¹f^h

k = σλ^β

λ+σfk+ λ^β+1 λ+σ

R^k+R^2M−k−1 1 +R^2M−1 fk

+λ^β+1

M−1

X

i=1

G(k, i;λ+σ)h(f_k−f_j).

Applying the triangle inequality, we obtain

λ^βB_h^x(B^x_h+λ)⁻¹f^h

k

≤ σλ^β λ+σ

fk

+ λ^β+1 λ+σ

fk

+λ^β+1

M−1

X

i=1

|G(k, i;λ+σ)|h|fk−fj|

≤h σλ^β

λ+σ+ λ^β+1

λ+σ+M(σ) λ^β+1

√λ+σ

M−1

X

i=1

R^|k−i||(k−i)h|^2βhi

kf^hk_C2β h

≤M₁(σ)kf^hk_C2β h

for anyλ >0 andx∈[0, l]. Therefore,f^h∈Eβ(Ch, B_h^x) and kf^hkE_β(C_h,B^x_h)≤M1(σ)kf^hk_C2β

h

.

Now, we prove the reverse inequality. For any positive operatorB_h^x, we can write v=

Z ∞ 0

M−1

X

i=1

G(k, i;λ+σ)B_h^x(B_h^x+λ)⁻¹f_ih₁dt.

Consequently, fk−fk+r=

Z ∞ 0

M−1

X

i=1

λ^−β[G(k+r, i;λ+σ)−G(k, i;λ+σ)]λ^βA^x_h(A^x_h+λ)⁻¹fih1dt, hence

|fk−fk+r| ≤ Z ∞

0

λ^−β

M−1

X

i=1

|G(k+r, i;λ+σ)−G(k, i;λ+σ)|h1dtkf^hk_Eβ(Ch,B^x

h). Let

Th=|rh1|^−2β Z ∞

0

λ^−β

M−1

X

i=1

|G(k+r, i;λ+σ)−G(k, i;λ+σ)|h1dt.

The proof of estimate

|f_k−f_k+r|

|rh1|^2β ≤T_hkf^hkE_β(C_h,B_h^x)

is based on the Lemmas 3.2 and 3.3. Thus, for any 1≤k < k+r≤N−1 we have established the inequality

|fk−fk+r|

|rh₁|^2β ≤ M

β(1−2β)kf^hk_Eβ(Ch,B^x

h). This means that

kf^hk_C2β

h ≤ M

β(1−2β)kf^hkE_β(C_h,B_h^x).

Theorem 3.5 in the casea(x) = 1 is proved. Now, leta(x) be continuous functions and letx, x₀∈[0,1] be arbitrary fixed points. Clearly, we have that

k(B^x_h−B^x_h⁰)(B^x_h⁰)kC_h→Ch ≤M.

From the formula

B^x_h(B_h^x+λ)⁻¹f^h=B_h^x0(B_h^x0+λ)⁻¹f^h

+λ(λ+B_h^x)⁻¹[B_h^x−B_h^x0](B^x_h⁰)⁻¹B_h^x0(B_h^x0+λ)⁻¹f^h

it follows that

λ^βB^x_h(B_h^x+λ)⁻¹f^h

≤ kf^hk_E

β(C_h,B^x_h⁰)+M λk(λ+B^x_h)⁻¹kC_h→Chkf^hk_E

β(C_h,B_h^x0)

≤M1kf^hk_E

β(C_h,B^x_h⁰). Then

kf^hk_Eβ(Ch,B^x

h)≤M1kf^hk_E

β(C_h,B^x_h⁰).

Theorem 3.5 is proved.

Theorem 3.6([3]). LetA_τ be the operator acting inE_τ =C[0, T]_τ defined by the formula Aτv^τ ={^vk−v_τ^k−1}^N₁, with v0 = 0. Then Aτ is a positive operator in the Banach spaceEτ =C[0, T]τ and

A^β_τf^τ=n 1 Γ(1−β)

k

X

r=1

Γ(k−m−β+ 1) (k−m)!

f_m−f_m−1 τ^β

oN 1

. By the definition of fractional difference derivative

D^β_τf^τ :=n 1 Γ(1−β)

k

X

r=1

Γ(k−m−β+ 1) (k−m)!

fm−f_m−1 τ^β

oN 1

.

Theorem 3.7. Let Aτ be the operator acting in Eτ = C[0, T]τ defined by the formula Aτv^τ(t) ={^vk−v_τ^k−1}^N₁ with the domain

D(A_τ) ={v^τ: v_k−v_k−1

τ ∈C[0, T]_τ, v₀= 0}.

ThenA is a positive operator in the Banach spaceE_τ =C[0, T]_τ, and A^β_τf^τ(t) =D^β_τf^τ(t)

for allf^τ(t)∈D(Aτ).

Thus, we have the following result on coercive stability of difference scheme (3.5).

Theorem 3.8. Let τ and hbe sufficiently small positive numbers and 0< β <1.

Then the solution of difference scheme (3.5)satisfies the following coercive stability estimate:

1≤k≤Nmax Bigknu^k_n+1−2u^k_n+u^k_n−1 h²

oM−1 n=1

k_Cβ[0,l]_h≤M(β) max

1≤k≤N

f_k^hk_Cβ[0,l]_h. Here,M(β)does not depend onτ, handf_k^h,1≤k≤N.

The proof of Theorem 3.8 is based on the Theorem 3.4 on positivity of difference space operatorB^x_h defined by formula ((3.1), on the Theorem 3.5 on the structure of fractional space Eβ(Ch, B_h^x0), on the Theorem 2.3 on connection of fractional derivatives with fractional powers of positive operators, on the Theorem 2.2 on spectral angle of fractional powers of positive operators, and on the Theorem 2.1 on fractional powers of coercively positive sums two operators.

4. A numerical application For numerical results, we consider the example

D^α_tu(t, x)−u_xx(t, x) +u(t, x) =f(t, x), f(t, x) = 6 sin²(πx)t^3−α

Γ(4−α) −2π²t³cos(2πx) +t³sin²(πx), 0< t <1, 0< x <1,

u(0, x) = 0,0≤x≤1, u(t,0) =u_x(t,1) = 0, 0≤t≤1

(4.1)

for the one-dimensional fractional parabolic partial differential equation with 0<

α < 1. The exact solution of problem (4.1) isu(t, x) =t³sin²πx. Note that this function is independent ofα, butf(t, x) depends onα.

Applying the difference scheme (3.3) for the numerical solution of (4.1), we obtain

1 Γ(1−α)

k

X

m=1

Γ(k−m−α+ 1) (k−m)!

u^m_n −u^m−1_n

τ^α −u^k_n+1−2u^k_n+u^k_n−1

h² +u^k_n=φ^k_n, φⁿ_k =f(t_k, x_n), t_k=kτ, 1≤k≤N, N τ =T,

xn=nh, 1≤n≤M −1, u⁰_n= 0, 0≤n≤M, u^k₀= 0, u^k_M−1=u^k_M, 0≤k≤N.

(4.2) We get the system of equations in the matrix form

AUn+1+BUn+CU_n−1=Dφn, 1≤n≤M−1,

U0=e0, UM−1=UM, (4.3)

where

A=























0 0 0 . . . 0 0

0 an 0 . . . 0 0

0 0 a_n . . . 0 0

. . . .

0 0 0 . . . an 0

0 0 0 . . . 0 an























(N+1)x(N+1)

,

B =























b₁₁ 0 0 . . . 0 0

b21 b22 0 . . . 0 0

b31 b32 b33 . . . 0 0

. . . . bN1 bN2 bN3 . . . bN N 0 b_N+1,1 b_N+1,2 b_N+1,3 . . . b_N+1,N b_N+1,N+1























(N+1)x(N+1)

,

C=























0 0 0 . . . 0 0

0 cn 0 . . . 0 0

0 0 c_n . . . 0 0

. . . .

0 0 0 . . . c_n 0

0 0 0 . . . 0 cn























(N+1)x(N+1)

,

D=























0 0 0 . . . 0 0

0 1 0 . . . 0 0

0 0 1 . . . 0 0

. . . .

0 0 0 . . . 1 0

0 0 0 . . . 0 1























(N+1)x(N+1)

,

φ_n=























 φ⁰_n φ¹_n φ²_n ... φ^N_n⁻¹

φ^N_n

























(N+1)x(1)

, U_n =

























 U_q⁰ U_q¹ U_q² ... U_q^N−1

U_q^N



























(N+1)x(1)

, q={n±1, n},

an=−1

h², cn=−1

h², b11= 1, b21=− 1

τ^α, b22= 1

τ^α + 1 + 2 h², b31=− Γ(2−α)

Γ(1−α)τ^α, b32= Γ(2−α) Γ(1−α)τ^α− 1

τ^α, b33= 1

τ^α + 1 + 2 h², and

bij =



















−Γ(1−α)(i−2)!τ^{Γ(i−1−α) α}, j= 1,

1 Γ(1−α)τ^α

h_{Γ(i−j+1−α)}

(i−j)! −^{Γ(i−j−α)}_(i−j−1)!i

, 2≤j≤i−2,

Γ(2−α)−Γ(1−α)

Γ(1−α)τ^α , j=i−1,

1

τ^α + 1 +_h²2, j=i,

0, i < j≤N+ 1

(4.4)

fori= 4,5, . . . , N+ 1 and φ^k_n =6 sin²(πnh)(kτ)^3−α

Γ(4−α) −2π²(kτ)³cos(2πnh) + (kτ)³sin²(πnh).

To solve the difference problem (4.3), a procedure of modified Gauss elimination method is applied. Hence, we seek a solution of the matrix equation in the following form:

U_j=α_j+1U_j+1+β_j+1, U_M = (I−α_M)⁻¹β_M, 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

αj+1=−(B+Cαj)⁻¹A,

βj+1= (B+Cαj)⁻¹(Dφ−Cβj), j= 1,2, . . . , M−1

wherej = 1,2, . . . , M−1,α1 is the (N + 1)×(N + 1) zero matrix, andβ1 is the (N+ 1)×1 zero matrix.

Second, applying the difference scheme (3.5), we obtain the second order of accuracy difference scheme intand inxand the Crank-Nicholson difference scheme for parabolic equations, one can represent the second order of accuracy difference scheme with respect intand inx

D^α_τu^k_n−1 2

hu^k_n+1−2u^k_n+u^k_n−1

h² +u^k−1_n+1−2u^k−1_n +u^k−1_n−1 h²

i +1

2 h

u^k_n+u^k−1_n i

=φ^k_n, φ^k_n=f(tk−τ

2, xn), tk=kτ, xn=nh, 1≤k≤N, 1≤n≤M −1,

u⁰_n= 0, 0≤n≤M,

u^k₀= 0, 3u^k_M −4u^k_M−1+u^k_M−2= 0, 0≤k≤N.

(4.5) HereD^α_τu^k_nis defined by (3.4) foru^k_n. We get the system of equations in the matrix form

AU_n+1+BU_n+CU_n−1=Dφ_n, 1≤n≤M−1,

U₀=e0, 3U_M−4U_M−1+U_M−2= 0, (4.6) where

A=























0 0 0 . . . 0 0

an an 0 . . . 0 0 0 a_n a_n . . . 0 0 . . . .

0 0 0 . . . an 0

0 0 0 . . . an an























(N+1)x(N+1)

,

B=























b₁₁ 0 0 . . . 0 0

b21 b22 0 . . . 0 0

b31 b32 b33 . . . 0 0

. . . . bN1 bN2 bN3 . . . bN N 0 b_N+1,1 b_N+1,2 b_N+1,3 . . . b_N+1,N b_N+1,N+1























(N+1)x(N+1)

,

C=























0 0 0 . . . 0 0

cn cn 0 . . . 0 0 0 c_n c_n . . . 0 0 . . . .

0 0 0 . . . c_n 0

0 0 0 . . . cn cn























(N+1)x(N+1)

,

D=























0 0 0 . . . 0 0

0 1 0 . . . 0 0

0 0 1 . . . 0 0

. . . .

0 0 0 . . . 1 0

0 0 0 . . . 0 1























(N+1)x(N+1)

,

φ_n=























 φ⁰_n φ¹_n φ²_n ... φ^N−1_n

φ^N_n

























(N+1)x(1)

, U_q=

























 U_q⁰ U_q¹ U_q² ... U_q^N−1

U_q^N



























(N+1)x(1)

, q={n±1, n},

an=− 1

2h², cn=− 1 2h², b11= 1, b21=−d 2^α−1

(2−α)(1−α)+ 1 h² +1

2, b22=d 2^α−1

(2−α)(1−α)+ 1 h²+1

2, b31=dh

(3/2)^5−α 1

1−α− 2

2−α+ 1 3−α

−73^2−α 2^3−α

1 (1−α)(2−α)

i , b32=dh

−3^4−α 2^3−α

1

1−α− 2

2−α+ 1 3−α

+3^2−α

2^−α

1 (1−α)(2−α)

i + 1

h² +1 2, b₃₃=dh3^4−α

2^5−α 1

1−α− 2

2−α+ 1 3−α

−3^2−α 2^3−α

1 (1−α)(2−α)

i+ 1 h² +1

2, b41=dh 1

1−αξ(1)− 1 2−αη(1)i

, b42=dh

− 5

1−αξ(1) + 2

2−αη(1)− 2^α−2 2−α

i , b43=dh 2

1−αξ(1)− 1

2−αη(1)− 2^α−1

1−α+ 2^α−1 2−α

i + 1

h² +1 2, b₄₄=dh2^α−1

1−α− 2^α−2 2−α

i + 1

h² +1

2, b₅₁=dh 2

1−αξ(2)− 1

2−αη(2)i ,

b₅₂=dh

− 5

1−αξ(2) + 2

2−αη(2) + 1

1−αξ(1)− 1 2−αη(1)i

, b₅₃=dh

− 3

1−αξ(1) + 2

2−αη(1) + 3

1−αξ(2)− 1

2−αη(2)− 2^α−2 2−α

i , b54=dh 2

1−αξ(1)− 1

2−αη(1)− 2^α−1

1−α+ 2^α−1 2−α

i + 1

h² +1 2, b55=dh2^α−1

1−α− 2^α−2 2−α

i + 1

h² +1 2, and

bij =

























































 dh

1

1−α(i−3)ξ(i−3)−_2−α¹ η(i−3)i

, j= 1,

dh

1

1−α(5−2i)ξ(i−3) +_2−α² η(i−3) +_1−α¹ (i−4)ξ(i−4)−_2−α¹ η(i−4)i

, j= 2,

dh

1

1−α(i−j+ 1)ξ(i−j)−_2−α¹ η(i−j)

+_1−α¹ (2j−2i+ 1)ξ(i−j−1) + _2−α² η(i−j−1) +_1−α¹ (i−j−2)ξ(i−j−2)−_2−α¹ η(i−j−2)i

, 3≤j≤i−3, dh

3

1−αξ(2)−_2−α¹ η(2)−_1−α³ ξ(1) +_2−α² η(1)−²_2−α^α−2i

, j=i−2,

dh_2ξ(1)

1−α −_2−α^η(1) −²_1−α^α−1 +²_2−α^α−1i

+_h¹2 +¹₂, j=i−1, dh

2^α−1

1−α −²_2−α^α−2i

+_h¹2 +¹₂, j=i,

0, i < j ≤N+ 1

fori= 6,7, . . . , N+ 1 and φ^k_n =6 sin²(πnh)(kτ)^3−α

Γ(4−α) −2π²(kτ)³cos(2πnh) + (kτ)³sin²(πnh).

For solving of the matrix equation (4.6), we use the same algorithm as in the (4.3) with

uM = [3I−4αM+α_M−1αM]⁻¹[(4I−αM−1)βM −βM−1].

Applying the difference schemes (3.3) and (3.5) for the numerical solution of (4.1), we constructed first and second order of accuracy difference schemes. The results of computer calculations show that the Crank-Nicholson difference scheme is more accurate than first order of accuracy difference scheme. Tables 11 and 2 are con- structed forN =M = 10,20,40,80, respectively.

Table 1. Error analysis of first and second order of accuracy dif- ference schemes forα= 1/2

Method N=M=10 N=M=20 N=M=40 N=M=80

1^st order of accuracy 1.1110 0.7049 0.3850 0.1998 2^nd order of accuracy 0.0953 0.0111 0.0017 3.332×10⁻⁴

Table 2. Error analysis of first and second order of accuracy dif- ference schemes forα= 1/3

Method N=M=10 N=M=20 N=M=40 N=M=80

1^st order of accuracy 1.1493 0.7333 0.4015 0.2086 2^ndorder of accuracy 0.1015 0.0121 0.0019 7.5456×10⁻⁴ Conclusion. In [12] the multidimensional fractional parabolic equation with the Dirichlet-Neumann conditions was studied. Stability estimates for the solution of the initial-boundary value problem for this fractional parabolic equation were given without proof. The stable difference schemes for this problem were presented.

Stability estimates for the solution of the first order of accuracy difference scheme were given without proof. The numerical result was given for the solution of first and second order of accuracy difference schemes of one-dimensional fractional parabolic differential equations without any discuss on the realization.

In the present study, coercive stability estimates for the solution of this initial- value problem for the fractional parabolic equation with the Dirichlet-Neumann conditions are established. Stable the first and second order of approximation in t and first order of approximation inx difference schemes for this problem are considered. Coercive stability estimates for the solution of the first order of accuracy difference scheme are obtained. A procedure of modified Gauss elimination method is applied for the solution of the first and second order of accuracy difference schemes of one-dimensional fractional parabolic differential equations. Moreover, applying this approach we can constructed the first and second of approximation intand a high order of approximation inxdifference schemes. Of course, coercive stability estimates for the solution of the first order of accuracy difference scheme can be obtained.

References

[1] R. P. Agarwal, B. de Andrade, C. Cuevas; Weighted pseudo-almost periodic solutions of a class of semilinear fractional differential equations. Nonlinear Analysis-Real World Applica- tions11(5) (2010), 3532–3554.

[2] D. Amanov, A. Ashyralyev;Initial-boundary value problem for fractional partial differential equations of higher order. Abstract and Applied Analysis (2012), Article ID 973102, doi:

10.1155/2012/973102, 1–14.

[3] A. Ashyralyev;A note on fractional derivatives and fractional powers of operators. Journal of Mathematical Analysis and Applications357(1) (2009), doi: 10.1016/j.jmaa.2009.04.012, 232–236.

[4] A. Ashyralyev;Well-posedness of fractional parabolic equations. Boundary Value Problems (2013), Article ID 31, doi: 10.1186/1687-2770-2013-31, 1–18.

[5] A. Ashyralyev;Well-Posedness of Parabolic Differential and Difference Equations with the Fractional Differential Operator. Malaysian Journal of Mathematical Sciences 6(S)(2012), 73–89.

[6] A. Ashyralyev;Well-posedness of the Basset problem in spaces of smooth functions. Applied Mathematics Letters24(7) (2011), doi: 10.1016/j.aml.2011.02.002, 1176–1180.

[7] A. Ashyralyev, B. Hicdurmaz;On the numerical solution of fractional Schrodinger differential equations with the Dirichlet condition. International Journal of Computer Mathematics89 (13-14) (2012), doi: 10.1080/00207160.2012.698841, 1927–1936.

[8] A. Ashyralyev, B. Hicdurmaz; A note on the fractional Schrodinger differential equations.

Kybernetes40(5-6) (2011), doi: 10.1108/03684921111142287, 736–750.

[9] A. Ashyralyev, D. Agirseven; Well-posedness of delay parabolic difference equations. Ad- vances in Difference Equations2014(18) (2014), doi: 10.1186/1687-1847-2014-18, 1–20.

[10] A. Ashyralyev, F. Dal;Finite difference and iteration methods for fractional hyperbolic partial differential equations with the neumann condition. Discrete Dynamics in Nature and Society (2012), Article ID 434976, doi: 10.1155/2012/434976.

[11] A. Ashyralyev, F. Dal, Z. Pinar; A note on the fractional hyperbolic differential and dif- ference equations. Applied Mathematics and Computation 217 (9) (2011), doi: 10.1016 /j.amc.2010.11.017, 4654–4664.

[12] A. Ashyralyev, N. Emirov, Z. Cakir; Fractional parabolic differential and difference equa- tions with the dirichlet-neumann conditions. First International Conference on Analysis and Applied Mathematics1470(2012), doi: 10.1063/1.4747641, 69–72.

[13] A. Ashyralyev, N. Nalbant, Y. Sozen;Structure of fractional spaces generated by second order difference operators. Journal of the Franklin Institute-Engineering and Applied Mathematics 351(2) (2014), doi: 10.1016/j.jfranklin.2013.07.009, 713–731.

[14] A. Ashyralyev, Y. A. Sharifov; Existence and uniqueness of solutions for the system of nonlinear fractional differential equations with nonlocal and integral boundary conditions.

Abstract and Applied Analysis (2012), doi: 10.1155/2012/594802, Article ID 594802.

[15] A. Ashyralyev, Z. Cakir;On the numerical solution of fractional parabolic partial differential equations with the dirichlet condition. Discrete Dynamics in Nature and Society (2012), doi:

10.1155/2012/696179, Article ID 696179.

[16] A. S. Berdyshev, A. Cabada, E. T. Karimov;On a non-local boundary problem for a parabolic- hyperbolic equation involving a Riemann-Liouville fractional differential operator. Nonlinear Analysis–Theory Methods and Applications 75(6) (2012), doi: 10.1016/j.na.2011.12.033, 3268–3273.

[17] S. Q. Bu, P. Clement, S. Guerre-Delabriere;Regularity of pairs of positive operators. Illinois Journal of Mathematics42(3) (1998), 357–370.

[18] Z. Cakir; Stability of difference schemes for fractional parabolic pde with the dirichlet- neumann conditions. Abstract and Applied Analysis (2012), doi: 10.1155/2012/463746, Ar- ticle ID 463746, 1–17.

[19] J. Y. Cao, C. J. Xu; A high order schema for the numerical solution of the frac- tional ordinary differential equations. Journal of Computational Physics 238(2013), doi:

10.1016/j.jcp.2012.12.013, 154–168.

[20] Ph. Clement; On (Lp-Lq) Coerciveness for a Class of Integrodifferential Equation on the Line. Two lectures at the 23rd annual Voronezh Winter School of Mathematics (USSR), Tech. Report 5-4-90, Voronezh, (1990).

[21] E. Cuesta, C. Lubich, C. Palencia;Convolution quadrature time discretization of fractional diffusion-wave equations. Mathematics of Computation75(254) (2006), doi: 10.1090/S0025- 5718-06-01788-1, 673–696.

[22] G. Da Prato, P. Grisvard; Sums of linear-operators and operational differential equations.

Journal De Mathematiques Pures Et Appliquees54(1975), 305–387.

[23] F. Dal;Application of variational iteration method to fractional hyperbolic partial differential equations. Mathematical Problems in Engineering (2009), doi: 10.1155/2009/824385, Article ID 824385.

[24] M. De la Sen; About robust stability of caputo linear fractional dynamic systems with time delays through fixed point theory. Fixed Point Theory and Applications (2011), doi:

10.1155/2011/867932, Article ID 867932.

[25] M. De la Sen; Positivity and stability of the solutions of caputo fractional linear time- invariant systems of any order with internal point delays. Abstract and Applied Analysis (2011), doi: 10.1155/2011/161246, Article ID 161246.

[26] M. De la Sen, R. P. Agarwal, A. Ibeas, et al;On the existence of equilibrium points, bound- edness, oscillating behavior and positivity of a sveirs epidemic model under constant and impulsive vaccination. Advances in Difference Equations (2011), doi: 10.1155/2011/748608, Article ID 748608.

[27] A. S. Erdogan, H. Uygun;A note on the inverse problem for a fractional parabolic equation.

Abstract and Applied Analysis (2012), doi: 10.1155/2012/276080, Article ID 276080.

[28] J. Fan, J. H. He; Fractal derivative model for air permeability in hierarchic porous media.

Abstract and Applied Analysis (2012), doi: 10.1155/2012/354701, Article ID 354701.

[29] A. A. Kilbas, H. M. Srivastava, J. J. Trujillo;Theory and Applications of Fractional Differ- ential Equations. Elsevier, Amsterdam (2006).