This site uses cookies for learning about our traffic, we store no personal details. ACCEPT COOKIES DECLINE COOKIES What are cookies?
univerge site banner
Review Article | Open Access | Aust. J. Eng. Innov. Technol., 2022; 4(5), 95-108 | doi: 10.34104/ijmms.022.0950108

Numerical Solution of Diffusion Equation with Caputo Time Fractional Derivatives Using Finite-Difference Method with Neumann and Robin Boundary Conditions

Shafiullah Niazai* Mail Img ,
Ariana Abdul Rahimzai Mail Img ,
Mustafa Danesh Mail Img ,
Burhanuddin Safi Mail Img

Abstract

Many problems in various branches of science, such as physics, chemistry, and engineering have been recently modeled as fractional ODEs and fractional PDEs. Thus, methods to solve such equations, especially in the nonlinear state, have drawn the attention of many researchers. The most important goal of researchers in solving such equations has been set to provide a solution with the possible minimum error. The fractional PDEs can be generally classified into two main types, spatial-fractional, and time-fractional differential equations. This study was designed to provide a numerical solution for the fractional-time diffusion equation using the finite-difference method with Neumann and Robin boundary conditions. The time fraction derivatives in the concept of Caputo were considered, also the stability and convergence of the proposed numerical scheme have been completely proven a numerical test was also designed and conducted to assess the efficiency and precision of the proposed method. Eventually, it can be said that based on the findings, the present technique can provide accurate results. 

INTRODUCTION

Fractional diffusion equations have significantly drawn the attention of many scholars due to their use in different branches of science, including their use in describing some phenomena in physics (Metzler and Klafter, 2000), chemistry and biochemistry (Yuste and Lindenberg, 2002), mechanical engine-ering (Magin et al., 2009), medicine (Chen et al., 2010), and electronics (Kirane et al., 2013). The non-local characteristic appears to be as property and the most important advantage of these equations, indicating that the state of a complex system depends not only on its current state but also on its previous states (Tamsir et al., 2021). Mathematical modeling has been recognized to be one of the strong and fundamental solutions for quantitative and qualitative analysis of such phenomena. In general, the quantitative and qualitative behavior characteristics of complex systems in various science and engineering problems can be much better un-derstood by mathematical modeling using fractional differential equations (Tamsir et al., 2021;  Demir et al., 2020;  Demir and Bayrak, 2019;  Demir et al., 2019). 

Different types of definitions have been suggested for fractional derivatives, including the Riemann-Liouville fractional derivatives and Caputo fractional derivatives as two important applications of them. Thus, the Caputo fractional derivative is a type of fractional derivative in mathematical modeling with experimental data analysis, which is broadly used in different branches of science. Since the analytical solution of the fractional differential equation, seems to be impossible in many phenomena and their numerical solving would highly matter (Demir et al., 2019; Yavuz et al., 2020; Usta and Sarikaya, 2019). 

Considerable numerical methods have been designed for diffusion fractional-time equations. Different authors have employed the finite element (Ford et al., 2011), compact finite element (Jacobs, 2016), Crank-Nicolson (Sayevand et al., 2016), the B-spline-based (Sweilam et al., 2016), and the implicit difference (Zhuang and Liu, 2006) methods to solve fractional-time equations. Murio et al. devised a stable unconditional implicit numerical method to solve the one-dimensional linear diffusion fractional-time equation on a finite medium (Murio, 2008). Using the generalized Euler method (GEM), Khader et al. proposed numerical methods for solving the fractional Riccati and Logistic differential equations based on Chebyshev approximations (Khader, 2011).  Celik et al. overcame to examine a numerical method for approximating a fractional diffusion equation, using the Riesz-fractional derivative on finite domains, which has second-order accuracy in terms of time and place. The “fractional central derivative” approach was also used to approximate the Riesz-fractional derivative and the Crank-Nicholson method, which has been applied to the fractional diffusion equation (Çelik and Duman, 2012). In a study, Lin et al. analyzed a stable and high-order scheme to effectively solve the fractional-time diffusion equation. Their proposed method relies on finite-difference time scheme and the Legendre spectral methods (Lin and Xu, 2007). A block-oriented finite-difference scheme has been suggested for solution of the fractional-time diffu-sion equation on non-uniform networks where un-conditional stability and convergence have been proven theoretically (Zhai and Feng, 2016). We considered a Caputo fractional-time diffusion equation in this article using the finite-difference method (FDM).

Consider the fractional-time diffusion equation as equation (1) with initial conditions (2) and boundary conditions (3).

(∂^α u(x,t))/〖∂(t)〗^α =c(x,t)  (∂^2 u(x,t))/〖∂x〗^2 +f(x,t),0

u(x,0)=g(x),                0

u(0,t)=0,   ωu(L,t)+(c(x,t)  ∂u(x,t)/∂x) │_(x=L)=y(t),       0

In these equations, the condition 0 <α <1 is al-ways established, C (x,t)is the continuous positive coefficient of the diffusion, f (x,t) is the source function, and g (x) is a sufficiently smooth function. The presented equations (1) to (3) are assumed to have a unique sufficiently smooth answer for the numerical analysis of the above fractional-time diffusion equation. In equation (3), the boundary conditions have been presented, which is governing this diffusion equation, now we express the Neu-mann boundary fraction conditions and Robin boun-dary fraction conditions as  ω=0, and ω>0, res-pectively. In these equations,   ∂^α/〖∂t〗^α   is the Caputos left fractional derivative (Roul and Goura, 2020; Pathi-ranage, 2021; Sun et al., 2013).

D_t^α u(x,t)=(∂^α u(x,t))/〖∂t〗^α =1/(Γ(1-α))+∫_0^t▒∂u(x,s)/∂s  ds/〖(t-s)〗^α   ,                                                  (4)

The Taylor expansion has been used in this study to discretize points u (x,t) with the third-order accu-racy followed by providing a sequence for stability and convergence of the proposed scheme to be proven. In general, a discretization process has been presented in the section-2 of this literature with an implicit finite-difference, which compatibility has been assessed as well and then in section-3, the stability and convergence of the approach has been proven, as well as a numerical example has been covered insection-4, to know the accuracy and efficiency of the scheme. Finally, the research general conclusion is provided in section-5.

Providing a model of implicit finite-difference ac-companied by evaluating its compatibility

Here, model of implicit finite-difference is pre-sented, accompanied by evaluating its compatibility. Considering Equations (1) to (3), the temporal and spatial partitions t_m and x_nwere defined as equation (5) for their discretization and numerical approxi-mation.

x_n=n×h(n=0,1,2,…,N)&t_m=m×τ(m=0,1,2,…,M)                                                          (5)

Where M and N are positive integers and according to the represented partitions, size of the temporal and the spatial networks in the examined discretization would be clearly equal to τ=T/M, and h=L/N, respectively. In the process of solving fractional-time diffusion equation in this literature, U_i^m  denotes the accurate answer and u_i^m represents the approximate answer resulting from the numerical method at point (x_i,t_m) of the network. Based on assumption, the short equivalent w_i^m   is used in this literature for simplifying the writing of equations for all para-meters such as w(x_i,t_m)=w_0. The lemmas 1 and 2 are used for its discretization by the finite-difference method (Sun & Wu, 2006; Tian et al., 2015).

Lemma 1

It was supposed that 0<α<1 and d(t)∈c^2 [0,t_m ]  are established, in such a case, the inequality in  (6)  would be established as follow.

|1/Γ(1-α)  ∫_0^(t_m)▒〖(d^ (s))/(t_m-s)^α  ds-τ^(-α)/Γ(2-α) 〗 (b_0 d(t_m )-∑_(j=1)^(m-1)▒〖(b_(m-j-1)-b_(m-1) 〗)d(t_j )-b_(m-1) d(t_0 )|≤1/Γ(1-α)  ((1-α)/12+2^(2-α)/(2-α)-(1+2^(-α) ))   max┬(0≤t≤t_m )⁡〖|d^ (t)| τ^(2-α),〗                             (6)

Where,  b_j=〖(j+1)〗^(1-α)-j^(1-α)  ,j=0,1,2,…

Lemma 2

It was supposed that d(x)∈L^1 (R)  and 〖-∞〗^(D_x^4 ) d(x) is established and its Fourier transform belongs to L^1 (R), and the weighted and transferred Grunwald-Letnikov operator is defined as equation (7) where p and q are integers.

L^(D_(h,p,q)^2 ) d(x)=λ_1/h^2  ∑_(j=1)^(m-1)▒〖g_j^((2) ) d(x-〗 (j-p)h)+λ_2/h^2  ∑_(j=0)^∞▒〖g_j^((2) ) d(x-(j-q)h),     〗                (7)

In this equation, we knew that λ_1=(2p-2)/(2(p-q))،λ_2=(2-2q)/(2(p-q))،p≠qand g_j^((2) ) 〖=(-1)〗^j (■(2@j))  are binominal fractional coefficients. Thus, we would have L^(D_(h,p,q)^2 ) d(x)=〖-∞〗^(D_x^2 ) d(x)+O(h^2). Therefore, we discredited the Caputo fractional-time derivative for each x∈R according to Lemma 1 as equation (8).

(∂^α u(x,t))/〖∂t〗^α  │_((x_i,t_(m+1)))=   τ^(-α)/(Γ(2-α)) ∑_(j=0)^m▒b_j  (U_i^(m+1-j)-U_i^(m-j) )+O(τ^(2-α) )  ,                         (8)

Then, according to Lemma 2, the spatial-derivative of equation (1), can be as approximation of equation (9), thus for p=1 and q=0, the spatial-derivative would be discretized as equation (10).

(∂^2 u(x,t))/〖∂x〗^2  │_((x_i,t_(m+1)))=   λ_1/h^2  ∑_(j=0)^(i+p)▒g_j^((2))  U_(i-j+p)^(m+1)+λ_2/h^2  ∑_(j=0)^(i+q)▒g_j^((2))  U_(i-j+q)^(m+1)+O(h^2 )  .                       (9)

(∂^2 u(x,t))/〖∂x〗^2  │_((x_i,t_(m+1)))=   1/h^2  ∑_(j=0)^(i+1)▒w_j^((2))  U_(i-j+1)^(m+1)+O(h^2 )  ,                                                               (10)

Where, w_0^((2))=2/2 g_0^((2)),w_j^((2))=g_j^((2))  ,j=1,2,…,M and the first order spatial-derivative on x=L are approxi-mated as equation (11).

∂u(x,t)/∂x │_((x_N,t_(m+1)))=   1/h ∑_(j=0)^(N+1)▒w_j  U_(N-j+1)^(m+1)+O(h^2 )  .                                                                   (11)

Then, we managed to discretize value of U_(N+1)^(m+1) as equation (12) using the third-order Taylor expansion to ultimately decompose equation (11) into equation (13).

U_(N+1)^(m+1)=  3U_N^(m+1)-3U_(N-1)^(m+1)+U_(N-2)^(m+1)+O(h^3 )  .                                                                          (12)

∂u(x,t)/∂x │_((x_N,t_(m+1) ) )=   1/h ∑_(j=1)^N▒w_j  U_(N-j+1)^(m+1)+w_0/h (3U_N^(m+1)-3U_(N-1)^(m+1)+U_(N-2)^(m+1) )+O(h^2 ).       (13)

It was then supposed that z=(τ^α Γ(2-α))/h^2   is established, therefore in the case, the implicit finite-difference scheme is written as equations (14) to (17).

U_i^1-〖zc〗_i^1 ∑_(j=0)^(N+1)▒w_j^((2) )  U_(i-j+1)^1=  U_i^0+^α (2-α) f_i^1  ,     1≤i≤N-1                                 (14)

U_i^(m+1)-〖zc〗_i^(m+1) ∑_(j=0)^1▒w_j^((2) )  U_(i-j+1)^(m+1)=(1-b_1 ) u_i^m+∑_(j=1)^(m-1)▒〖(b_j 〗+b_(j+1)) u_i^(m-j)+b_m u_i^0+τ^α Γ(2-α) f_i^(m+1),1≤m≤M-1    &1≤i≤N-1                                                                                   (15)

U_0^(m+1)=0,                        0≤m≤M-1    & i=N                                                                      (16)

〖ωu〗_N^(m+1)+(c_N^(m+1))/h ∑_(j=1)^N▒w_j  U_(N-j+1)^(m+1)+(c_N^(m+1) w_0)/h (3U_N^(m+1)-3U_(N-1)^(m+1)+U_(N-2)^(m+1) )=y^(m+1),0≤m≤M-1    & i=N                                                                                                    (17)

After analyzing local truncation-error for 1≤i≤N with R_i^m, the value of error can be seen with using equations (8), (10), and (13) as follows. In fact, this issue implies that, compatibility of implicit finite-difference schemes defined by (14) to (17) is:

R_i^1=U_i^1-zc_i^1 ∑_(j=0)^i▒w_j^((2) )  U_(i-j+1)^1-U_i^0- τ^α Γ(2-α) f_i^1  = O(τ^(2-α)+h^2 ),1≤i≤N-1,            (18)

R_i^(m+1)=U_i^(m+1)-zc_i^(m+1) ∑_(j=0)^i▒w_j^((2) )  U_(i-j+1)^(m+1)-(1-b_1 ) U_i^m-∑_(j=1)^(m-1)▒(b_j-b_(j+1) )  U_i^(m-j)-b_m U_i^0-τ^α Γ(2-α) f_i^(m+1)=O(τ^(2-α)+h^2 ),    1≤i≤N-1,      1≤m≤M-1,                                                                                                                                  (19)

R_N^(m+1)=〖ωU〗_N^(m+1)+(c_N^(m+1))/h ∑_(j=1)^N▒w_j  U_(N-j+1)^(m+1)+(c_N^(m+1) w_0)/( h)(3U_N^(m+1)-3U_(N-1)^(m+1)+U_(N-2)^(m+1))-y^(m+1)=O(h^2 ),0≤m≤M-1.                                                                                                                     (20)

The stability analysis and convergence

In this section of literature, stability and convergence of finite numerical difference has been discussed for solving the fractional-time diffusion equation. The columnar vectors〖 U〗^m, Q^(m-1), and F^m  are defined in the form of a matrix shown in (21)  for elucidation of the text. By using provided matrix definition, we expressed the implicit finite-difference scheme discretized in Equations (14) to (17) as a matrix shown in Equations (22) and (23).

{█(U^m= 〖〖(u〗_1^m,u_2^m,…,u_N^m)〗^T,                                                                                      @V^(m-1)= 〖〖(u〗_1^(m-1),u_2^(m-1),…,u_(N-1)^(m-1),0)〗^T,                                                                 @W^m= 〖〖(τ^α Γ(2-α)f〗_1^m,〖τ^α Γ(2-α)f〗_2^m,…,〖τ^α Γ(2-α)f〗_(N-1)^m,hy^m)〗^T,        )┤     for  1≤m≤M                                                                                                                                                 (21)

AU^1=V^0+W^1  ,1≤m≤M-1                                                                                                                           (22)

AU^(m+1)=〖(1-b_1)V〗^m+〖∑_(j=1)^(m-1)▒〖(b_j-b_(j+1))〗 V〗^(m-j)+ b_m V^0  +W^(m+1),1≤m≤M-1          (23)

The coefficients of elements in above matrix can be as form of matrix A, which is elucidated in equation (24).

A={■(■(zc_i^(m+1) w_(i-j+1)^((2) ),   1≤j≤i-1,1≤i≤N-1,@1-zc_i^(m+1) w_1^((2) ),                           1≤j=i≤N-1,)@■(-zc_i^(m+1) w_0^((2) ),          j=i+1,           1≤i≤N-1,@0,                       i+2≤j≤N,         1≤i≤N-2,@c_N^(m+1) w_(N-j+1),                1≤j≤N-3,           i=N,)@■(c_N^(m+1) w_3+c_N^(m+1) w_0,     j=N-2,    i=N,               @c_N^(m+1) w_2-〖3c〗_N^(m+1) w_0,    j=N-1,    i=N,               @hω+c_N^(m+1) w_1+〖3c〗_N^(m+1) w_0,   j=i=N.                      ))┤                                                                    (24)  

Some lemmas, as bellow are needed to evaluate the sustainability of the mentioned scheme (Samko et al., 1993; Liu et al., 2015; Lin and Xu, 2007).

Lemma 3

It was supposed that ρ includes positive real numbers andn≥1 is an integer. In such a case, the coefficients f_j^((ρ) ) (j=0,1,..)  would have properties as following:

f_0^((ρ) )=1,f_j^((ρ) )=(1-(ρ+1)/j) f_(j-1)^((ρ) )              for j≥1,                                                                                                                                                          (i)

f_1^((ρ) )0         for  0<ρ<1,(ii)

f_2^((ρ) )>f_3^((ρ) )>⋯>0,∑_(j=0)^n▒f_j^((ρ) ) >0        for  1<ρ<2,                                                                                      (iii)

∑_(j=0)^n▒f_j^((ρ) ) =(-1)^n (■(ρ-1@n)),                                                                                                                                           (iv)

∑_(j=0)^n▒f_j^((ρ_1 ) ) =f_(n-j)^((ρ_2 ) )=f_n^((〖ρ_1+ρ〗_2 ) ).                                                                                                                                      (v)

Lemma 4

It was supposed that ρ is positive real value, where the case,〖 w〗_j^((ρ) ) (j=0,1,…)  would have some of the properties, as expressed bellow.

w_0^((ρ) )=ρ/2,w_0^((ρ) )=(2-ρ-ρ^2)/2,w_2^((ρ) )=ρ(ρ^2+ρ-4)/2,                                                                                     (i)

w_0^((ρ) )=ρ/2 f_j^((ρ) )+(2-ρ)/2 f_(j-1)^((ρ) ),for j≥3,                                                                                                            (ii)

〖 w〗_2^((ρ) )≤〖 w〗_3^((ρ) )≤⋯≤0,∑_(j=0)^n▒〖w_j^((ρ) )>0,〗      for 0<ρ<1,n≥1,                                                       (iii)

〖1>w〗_0^((ρ) )>〖 w〗_3^((ρ) )>〖 w〗_4^((ρ) )…>0,∑_(j=0)^n▒〖w_j^((ρ) )<0,〗     for 1<ρ<2,n≥2.                                               (iv)

Lemma 5

If ρ_1 and ρ_2  are supposed as two constants, in such case, the coefficients b_j (j=1,2,…)  would be dealt with properties, written below.

b_(j>0),                                                                                                                                                                                    (i)

b_j>b_(j+1),                                                                                                                                                                            (ii)

ρ_1 j^α≤(b_j )^(-1)≤ρ_2 j^α.                                                                                                                                                   (iii)

We considered the error value ε_i^m=u_i^m-u ̃_i^m, where stability of the numerical method of implicit finite-difference provided scheme for solving the fractional-time diffusion equation is evaluated. Hence, u ̃_i^m  is obtained as the approximate answer of the differential scheme with the initial condition u ̃_i^0. In such a situation, the error value of the matrix form will be as defined in equation (25).

ε^m=〖(ε_1^m,ε_1^m,…,ε_N^m)〗^T,‖ε^m ‖_∞=max┬(1≤i≤N)⁡|ε_i^m |.                                                                                        (25)

Then, according to the definition of the finite-difference scheme provided in this literature, we obtain:

ε_i^1-〖zc〗_i^1 ∑_(j=0)^i▒w_j^((2) )  ε_(i-j+1)^(m+1)=  ε_i^0  ,1≤i≤N-1,1≤m≤M-1,                                                       (26)

ε_i^(m+1)-〖zc〗_i^(m+1) ∑_(j=0)^i▒w_j^((2) )  ε_(i-j+1)^(m+1)=  (1-b_1 )ε_i^m+∑_(j=1)^(m-1)▒〖(b_j-b_(j+1) )ε_i^(m-j)+〖b_m ε〗_i^0 〗  ,1≤i≤N-1    ,1≤m≤M-1,                                                                                               (27)

〖ωε〗_N^(m+1)+(c_N^(m+1))/h ∑_(j=1)^N▒w_j  ε_(N-j+1)^(m+1)+(c_N^(m+1) w_0)/h (3ε_N^(m+1)-│3ε_(N-1)^(m+1)+ε_(N-2)^(m+1) )=0,i=N,ε_0^m (m=0,1,…,M),0≤m≤M-1.                                                                                                                               (28)

Then, using equation (28), we obtained equation (29). Afterward, considering the i=N-1, we achieved equations (30) and (31) by substituting equation (29) in equations (26) and (27).

ε_N^(m+1)=(-c_N^(m+1) ∑_(j=2)^N▒〖w_j ε_(N-j+1)^(m+1)+3c_N^(m+1) w_0 ε_(N-1)^(m+1)-c_N^(m+1) w_0 ε_(N-2)^(m+1) 〗)/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )  ,                                                           (29)

ε_(N-1)^1-zc_(N-1)^1 ∑_(j=1)^(N-1)▒〖(w〗_(N-1)^((2) )  〖-s^1 c〗_N^1 w_0^((2) ) w_(N-j+1))ε_j^1-3z〖〖s^1 c〗_(N-1)^1 c_N^1 w〗_0^((2) ) w_0 ε_(N-1)^1+z〖〖s^1 c〗_(N-1)^1 c_N^1 w〗_0^((2) ) w_0 ε_(N-2)^1=ε_(N-1)^0,                                                                                                                                                (30)

ε_(N-1)^(m+1)-zc_(N-1)^(m+1) ∑_(j=1)^N▒〖(w〗_(N-1)^((2) )  〖-s^(m+1) c〗_N^(m+1) w_0^((2) ) w_(N-j+1))ε_j^(m+1)-3z〖〖s^(m+1) c〗_(N-1)^(m+1) c_N^(m+1) w〗_0^((2) ) w_0 ε_(N-1)^(m+1)+z〖〖s^(m+1) c〗_(N-1)^(m+1) c_N^(m+1) w〗_0^((2) ) w_0 ε_(N-2)^(m+1)=(1-b_1 )ε_(N-1)^m+∑_(j=1)^(m-1)▒〖(b_j-b_(j+1) ) ε_(N-1)^(m-j)+b_m ε_(N-1)^0 〗,                                                               (31)

Following the discussion and we prove theorem-1 for stability of the mentioned numerical method as following.

Theorem 1

 The implicit finite-difference presented for the fractional-time diffusion equation (14) - (17) is un-conditionally stable.

Proof

We first defined matrix of the function ξ in the form ofξ^m=〖(ε_1^m,ε_2^m,…,ε_(N-1)^m)〗^T,1≤m≤M. In such a case, equations (26),(27),(30) and (31) can be obtained as follows.

Bξ^1=ξ^0  ,1≤m≤M-1 ,                                                                                                                                (32)

Bξ^(m+1)=(1-b_1 ) ξ^m+∑_(j=1)^(m-1)▒〖(b_j-b_(j+1) ) ξ^(m-j)+b_m ξ^0 〗  ,1≤m≤M-1,                                              (33)

Then, the coefficients of elements in matrix B will be calculated as shown in equation (34).

B={■(-zc_i^(m+1) w_(i-j+1)^((2) ),1≤j≤i-1,1≤i≤N-1,                                             @■(1-zc_i^(m+1) w_1^((2) ),1≤j=i≤N-2,                                                                         @■(-zc_i^(m+1) w_0^((2) ),j=i+1,1≤i≤N-1,                                                          @■(0,i+2≤j≤N-1,1≤i≤N-3,                                                            @■(〖-zc〗_(N-1)^(m+1) 〖(w〗_(N-j)^((2) )-〖s^(m+1) c〗_N^(m+1) w_0^((2) ) w_(N-j+1)),1≤j≤N-3,i=N-1,@■(@■(〖-zc〗_(N-1)^(m+1) 〖(w〗_2^((2) )-〖s^(m+1) c〗_N^(m+1) w_0^((2) ) w_3-〖s^(m+1) c〗_N^(m+1) w_0^((2) ) w_0),j=N-2,i=N-1,@■(@■(〖1-zc〗_(N-1)^(m+1) 〖(w〗_1^((2) )-〖s^(m+1) c〗_N^(m+1) w_0^((2) ) w_2+3〖s^(m+1) c〗_N^(m+1) w_0^((2) ) w_0),i=j=N-1.       @)))))))))┤                             (34)

For proof of the theorem, firstly it is needed to show that matrix B, does not have unique eigenvalue. To prove this issue, we indicated that eigenvalues of matrix B fall within circles, having center b_(i,j) and radius ∑_(j=1,j≠i)^(N-1)▒|b_(i,j) | ,  according to Gerschgorin theorem (Xie and Fang, 2019). Thus, based on Lemma 4:

b_(i,j)≤0 ,    j≠i,1≤i≤N-1,                                                                                                                        (35)

b_(i,j)≥1,1≤i≤N-1,                                                                                                                                                (36)

The result will be as:

b_(i,j)-∑_(j=1,j≠i)^(N-1)▒|b_(i,j) | =∑_(j=1)^(N-1)▒b_(i,j) =1-zc_i^(m+1) ∑_(j=0)^i▒w_j^((2) ) >1.                                                                                     (37)

On the other hand, according to Lemma 3 and Lemma 4, it can be written that:

s^(m+1) c_N^(m+1)=(c_N^(m+1))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )  ≤  (c_N^(m+1))/(c_N^(m+1) w_1+3c_N^(m+1) w_0 )=1/(w_1+3w_0 )<1.                                  (38)

Then, direct calculation led to:

w_3+w_0>0,〖-w〗_2^((2) )+w_0^((2) ) w_3+w_0^((2) ) w_0<0,〖-w〗_1^((2) )+w_0^((2) ) w_2+〖3w〗_0^((2) ) w_0>0.                        (39)

Therefore, we reached the following equations:

b_(N-1,j)=zc_(N-1)^(m+1) (〖w_(N-j)^((2) )-s^(m+1) c〗_N^(m+1) w_0^((2) ) w_(N-j+1) )≤0 ,1≤j≤N-3,i=N-1,                          (40)

b_(N-1,N-2)=zc_(N-1)^(m+1) (-〖w_2^((2) )+s^(m+1) c〗_N^(m+1) w_0^((2) ) (w_3+w_0  ))

b_(N-1,N-2)=1+zc_(N-1)^(m+1) (-〖w_1^((2) )+s^(m+1) c〗_N^(m+1) w_0^((2) ) w_2-3s^(m+1) c_N^(m+1) w_0^((2) ) w_0 )>1+zc_(N-1)^(m+1) (〖-w_1^((2) )+w〗_0^((2) ) w_2-3w_0^((2) ) w_0 )>1 ,i=N-1.                                       (42)

The above three equations came to the conclusion below:

b_(N-1,N-1)-∑_(j=1)^(N-2)▒|b_(i,j) | =∑_(j=1)^(N-1)▒b_(i,j) =1+zc_(N-1)^(m+1) (-∑_(j=1)^(N-1)▒w_j^((2) )  〖+s^(m+1) c〗_N^(m+1) w_0^((2) ) ∑_(j=2)^N▒w_j -2s^(m+1) c_N^(m+1) w_0^((2) ) w_0 )>1+zc_(N-1)^(m+1) (w_0^((2) )-s^(m+1) c_N^(m+1) w_0^((2) ) w_0+w_1)-2s^(m+1) c_N^(m+1) w_0^((2) ) w_0 )>1+zc_(N-1)^(m+1) (w_0^((2) )-(w_0^((2) ) (c_N^(m+1) w_1+3c_N^(m+1) w_0 ))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 ))>1.                                                        (43)

The inequalities shown in equations (37) and (43) suggests that matrix B, does not have unique eigenvalue. Thus, spectral radius of matrix  B^(-1), was found to be less than one, so by following the process, it is advised to show that equation (44) is established.

‖ξ^m ‖_∞  ≤   ‖ξ^0 ‖_∞,1≤m≤M.                                                                                                                     (44)

In fact, according to (32), we are writing:

‖ξ^1 ‖_∞  ≤   ‖ξ^0 ‖_∞,                                                                                                                                                      (45)

Assuming the establishment of equation (46) and according to equation(33), equation(47), it is obtained as:

‖ξ^k ‖_∞  ≤   ‖ξ^0 ‖_∞,     k=2,3,…,m.                                                                                                                     (46)

‖ξ^(k+1) ‖_∞  ≤   〖‖(1-b_1 ) ξ^k+∑_(j=1)^(k-1)▒(b_j-b_(j+1) )  ξ^(k-1)+b_k ξ^0 ‖_∞≤〗_∞ (1-b_1 ) ‖ξ_1^k ‖_∞+∑_(j=1)^(k-1)▒〖(b_j-b_(j+1) ) ‖ξ_1^(k-j) ‖_∞+b_k ‖ξ_1^0 ‖_∞≤((1-b_1 )+∑_(j=1)^(k-1)▒(b_j-b_(j+1) ) +b_k ) ‖ξ^0 ‖_∞=‖ξ^0 ‖_∞ 〗.                                                                                                                                             (47)

By using mathematical induction, it is managed to obtain equation (44) and on the other hand, value of |ε_N^(m+1) |  is calculated as:

|ε_N^(m+1) |=  |c_N^(m+1) ∑_(j=2)^N▒〖w_j ε_(N-j+1)^(m+1)-3c_N^(m+1) w_0 ε_(N-1)^(m+1) 〗+c_N^(m+1) w_0 ε_(N-2)^(m+1) |/|hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 | ,≤|(c_N^(m+1) ∑_(j=2)^N▒〖|w_j ||ε_(N-j+1)^(m+1) |+3c_N^(m+1) w_0 |ε_(N-1)^(m+1) | 〗+c_N^(m+1) w_0 |ε_(N-2)^(m+1) |)|/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 ),≤(c_N^(m+1) (∑_(j=2)^N▒|w_j | +4w_0 ))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )   〖max〗┬(1≤i≤N-1)⁡|ε_i^(m+1) |≤(c_N^(m+1) (w_1+5w_0 ))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )   〖max〗┬(1≤i≤N-1)⁡|ε_i^(m+1) |,≤(5/3 c_N^(m+1) (w_1+3w_0 ))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )   〖max〗┬(1≤i≤N-1)⁡〖|ε_i^(m+1) |≤5/3〗  〖max〗┬(1≤i≤N-1) |ε_i^(m+1) |≤5/3 ‖ξ^(m+1) ‖_∞≤5/3 ‖ξ^0 ‖_∞,i=N .                                                                                                                                                   (48)

Using equations (44) and (48), then it is obtained that:

‖ε^(m+1) ‖_∞=max{‖ξ^(m+1) ‖_∞,|ε_N^(m+1) |}≤C‖ξ^0 ‖_∞≤C‖ε^0 ‖_∞.                                                                            (49)

Here C is positive constant coefficient and it is in-dependent of τ and h. Hence, the argument was completed according to (Smith et al., 1985; Yu and Tan, 2003). Afterward, we defined e_i^m and its matrix as form of equations (50) and (51), to continue the assessment of numerical convergence of the method presented in this literature.

e_i^m=U_i^m-u_i^m,1≤i≤N,0≤m≤M,                                                                                                         (50)

e^m=〖(e_1^m,e_2^m,…,e_N^m)〗^T,‖e^m ‖_∞=max┬(1≤i≤N)⁡〖|e_i^m |.〗                                                                                                 (51)

Alternatively, according to the definition of finite-difference scheme for 1≤i≤N-1 and 1≤m≤M-1, we achieved following equations:

e_i^1-zc_i^1 ∑_(j=0)^i▒〖w_j^((2) ) e_(i-j+1)^1 〗=e_i^0+τ^α Γ(2-α)R_i^1,                                                                                                     (52)

e_i^(m+1)-zc_i^(m+1) ∑_(j=0)^i▒〖w_j^((2) ) e_(i-j+1)^(m+1) 〗=(1-b_1 ) e_i^m+∑_(j=1)^(m-1)▒〖(b_j-b_(j+1) )e_i^(m-j)+b_m e_i^0 〗+τ^α Γ(2-α) R_i^(m+1),                                                                                                                          (53)

Where, i=N and 0≤m≤M-1, ε_0^m (m=0,1,…,M):

〖ωe〗_N^(m+1)+(c_N^(m+1))/h ∑_(j=1)^N▒〖w_j e_(N-j+1)^(m+1)+(c_N^(m+1) w_0)/h(3e_N^(m+1) 〗-3e_(N-1)^(m+1)+e_(N-2)^(m+1))=R_N^(m+1),                                               (54)

Based on (54), the equation below was obtained:

e_N^(m+1)=(-c_N^(m+1) ∑_(j=2)^N▒〖w_j e_(N-j+1)^(m+1)+3c_N^(m+1) w_0 e_(N-1)^(m+1)-c_N^(m+1) w_0 e_(N-2)^(m+1)+hR_N^(m+1) 〗)/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 ),                                                (55)

Supposing, that i=N-1 is established and equations (56) and (57) are obtained by placing (55) in relations (52) and (53).

e_(N-1)^1-zc_(N-1)^1 ∑_(j=1)^(N-1)▒〖〖(w〗_(N-j)^((2) )-s^1 c_N^1 w_0^((2) ) w_(N-j+1))e_j^1 〗-3zs^1 c_(N-1)^1 c_N^1 w_0^2 w_0 e_(N-1)^1+zs^1 c_(N-1)^1 〖c_N^1 w〗_0^((2) ) w_0 e_(N-2)^1=e_(N-1)^0+τ^α Γ(2-α) R_(N-1)^(-m+1),                                                                                                           (56)

e_(N-1)^(m+1)-zc_(N-1)^(m+1) ∑_(j=1)^(N-1)▒〖〖(w〗_(N-j)^((2) )-s^(m+1) c_N^(m+1) w_0^((2) ) w_(N-j+1))e_j^(m+1) 〗-3zs^(m+1) c_(N-1)^(m+1) c_N^(m+1) w_0^((2) ) w_0 e_(N-1)^(m+1)+zs^(m+1) c_(N-1)^(m+1) 〖c_N^(m+1) w〗_0^((2) ) w_0 e_(N-2)^(m+1)=(1-b_1 )e_(N-1)^m+∑_(j=1)^(m-1)▒〖(b_j-b_(j+1) ) e_(N-1)^(m-j) 〗+b_m e_(N-1)^0+τ^α Γ(2-α) R_(N-1)^(-m+1),                           (57)

Where, R_(N-1)^(-m+1)≤C(τ^(2-α)+h^2) is established. Now here we present and prove Theorem-2 for convergence of the method.

Theorem 2 

The implicit finite-difference scheme expressed in equations (14) to (17) is an unconditional convergence and the constant α constant is independent of τ and h so that –

‖U_i^m-u_i^m ‖_∞≤C(τ^(2-α)+h^2 ),     1≤m≤M                                                                                                          (58)

Proof

We first put ζ^m=〖(e_1^m,e_2^m 〖,…,e〗_(N-1)^m)〗^T,1≤m≤M. Therefore, we could write equations (52), (53), (56), and (57) as follow:

Bζ^(m+1)=ζ^0+R^1,                                                                                                                                                            (59)

Bζ^(m+1)=〖(1-b_1)ζ〗^m+∑_(j=1)^(m-1)▒(b_j-b_(j+1) )  ζ^(m-j)+b_m ζ^0+R^m,                                                                            (60)

Where〖,R〗^m=τ^α Γ(2-α)〖(R_1^m,R_2^m,…,R_(N-2)^m,R_(N-1)^(-m+1))〗^T. In case of m = 0, it is followed as below:

max┬(1≤i≤N-1)⁡|e_i^1 |=‖ζ^1 ‖_∞≤‖ζ^0 ‖_∞+‖R^1 ‖_∞≤Cτ^α (τ^(2-α)+h^2 )=b_0^(-1) Cτ^α (τ^(2-α)+h^2 ),                                (61)

Where, C is positive constant, which is independent of τ and h. Now it is supposed that:

max┬(1≤i≤N-1)⁡|e_i^k |=‖ζ^k ‖_∞≤b_(k-1)^(-1) Cτ^α (τ^(2-α)+h^2 ),                                                                                                       (62)

Where, k=2,3,…,m and Cis constant independent of τ and h, due to b_m≤b_k≤1:

b_m^(-1)≥b_k^(-1).                                                                                                                                                                        (63)

Therefore, equation (62) can be expressed as form of equation (64).

max┬(1≤i≤N-1)⁡|e_i^k |=‖ζ^k ‖_∞≤b_m^(-1) Cτ^α (τ^(2-α)+h^2 ).                                                                                                         (64)

Therefore:

max┬(1≤i≤N-1)⁡|e_i^(k+1) |=‖ζ^(k+1) ‖_∞≤‖(1-b_1 ) ζ^k+∑_(j=1)^(k-1)▒(b_j-b_(j+1) )  ζ^(k-j)+b_k ζ^0+R^(k+1) ‖_∞,≤(1-b_1 ) ‖ζ^k ‖_∞+∑_(j=1)^(k-1)▒(b_j-b_(j+1) )  ‖ζ^(k-1) ‖_∞+‖R^(k+1) ‖_∞,≤((1-b_1 )+∑_(j=1)^(k-1)▒(b_j-b_(j+1) ) +b_m ) b_m^(-1) Cτ^α (τ^(2-α)+h^2 ),

The equation below can be obtained by the method of mathematical induction:

max┬(1≤i≤N-1)⁡|e_i^(m+1) |=‖ζ^(m+1) ‖_∞≤b_m^(-1) Cτ^α (τ^(2-α)+h^2 )≤C_1 (mτ)^α (τ^(2-α)+h^2 )≤C_2 (τ^(2-α)+h^2 ),                                                                                                                                (66)

Where, 0≤m≤M-1, and〖 C〗_1 and C_2  are positive constants independent of τ and h. also for i=N, we would have: 

|e_N^(m+1) |=|-c_N^(m+1) ∑_(j=2)^N▒〖w_j e_(N-j+1)^(m+1)+3c_N^(m+1) w_0 e_(N-1)^(m+1)-c_N^(m+1) w_0 e_(N-2)^(m+1)+hR_N^(m+1) 〗|/|hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 | ,≤  (c_N^(m+1) ∑_(j=2)^N▒〖|w_j ||e_(N-j+1)^(m+1) |+3c_N^(m+1) w_0 e_(N-1)^(m+1)+c_N^(m+1) w_0 |e_(N-2)^(m+1) |+h|R_N^(m+1) | 〗)/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )  ,≤(c_N^(m+1) (∑_(j=2)^N▒〖|w_j |+4w_0)〗)/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )   max┬(1≤i≤N-1)⁡〖|e_i^(m+1) |+h|R_N^(m+1) |/(c_N^(m+1) w_1+3c_N^(m+1) w_0 )  ,〗≤(c_N^(m+1) (w_1+5w_0 ))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )   max┬(1≤i≤N-1)⁡〖|e_i^(m+1) |+h|R_N^(m+1) |/(c_N^(m+1) w_1+3c_N^(m+1) w_0 )  ,≤(5/3 c_N^(m+1) (w_1+3w_0 ))/(hω+c_N^(m+1) w_1+3c_N^(m+1) w_0 )   max┬(1≤i≤N-1)⁡〖|e_i^(m+1) |+h|R_N^(m+1) |/(c_N^(m+1) w_1+3c_N^(m+1) w_0 )  ,≤〖5/3  max┬(1≤i≤N-1)〗⁡〖|e_i^(m+1) |+h|R_N^(m+1) |/(c_N^(m+1) w_1+3c_N^(m+1) w_0 )≤〖5/3 ‖ξ^(m+1) ‖_∞〗⁡〖+h|R_N^(m+1) |/(c_N^(m+1) w_1+3c_N^(m+1) w_0 )≤C_3 〗 〗 〗 〗 (τ^(2-α)+h^2 ).                                                                                                                                                      (67)

Where, 0≤m≤M-1 and C_3  is a positive coefficient, which is independent of τ and h. Based on (66) and (67), it was followed as below:

‖e^(m+1) ‖_∞=max{‖ζ^(m+1) ‖_∞,|e_N^(m+1) |}≤C (τ^(2-α)+h^2 ).                                                                                       (68)

Where, C is positive constant and it is independent of τ and  h. Thereby, the argument was completed.

Numerical examples

The efficiency and consistency of numerical method, presented for the Caputo fractional-time diffusion equation is determined in this section, which is followed by considering the maximum error indices〖 L〗_2 and L_max for this evaluation.

L_2=√(h∑_(j=0)^M▒|u_j^exact-u_j^numerical |^2 )   ,L_max=〖max〗┬(0≤j≤M)⁡|u_j^exact-u_j^numerical |,                                             (69)

The ROC of this problem is calculated as follow:

ROC=(log⁡(E^(h_1 )/E^(h_2 )))/(log⁡(h_1/h_2)),       ROC=(log⁡(E^(k_1 )/E^(k_2 )))/(log⁡(k_1/k_2))  ,                                                                                                (70)

Where, the errors E^(h_1 ) and E^(h_2 ) of size h_1 and h_2 are presented here respectively and as well as E^(k_1 ) and E^(k_2 ) indicate, errors in the mesh sizes k_1 and k_2 respectively.

Example

We initially considered the Caputo fractional-time diffusion equation:

(∂^α u(x,t))/〖∂(t)〗^α =c(x,t)  (∂^2 u(x,t))/〖∂x〗^2 +f(x,t),0

Considering the:

f(x,t)=(2/(Γ(3-α)) t^(2-α)+4π^2 t^2 ),     c(x,t)=1,     y(x,t)=2πt^2,    u(0,t)=0,     u(x,0)=0.

The accurate answer will be equal to  u(x,t)= t^2 sin⁡(2πx).

Result of the numerical and analytical solutions of this problem is illustrated in Figures 1-3 by using the method provided in this literature.

Fig. 1: Shows result of analytical and approximate solutions and as well as absolute error with Δx =2^(-10) and Δt = 2^(-7).

Fig. 2: Shows results of analytical and approximate solution at x≈1 for Δx =2^(-10) and  Δt = 2^(-7)at different t_S.

Fig. 3: Here result of analytical and approximate solution at  T=0.5 ,T=1 for Δx =2^(-10) and  Δt = 2^(-7) at different x_s  is shown.

In this literature, the Caputo temporal-fractional diffusion equation for different values of M and N is evaluated and Table 1 lists the error values〖 L〗_2  and L_max  for varying values of ∆t,∆x for α=0.5. Ac-cording to output of the problem, the error value decreases with the increased size of temporal and spatial meshes. However, variation of changes in the spatial meshes on the performance of the model seems to be higher, therefore the results reveal very good consistency between analytical and approxi-mate solutions.

Table 1: The error values of L_2 and L_max  are shown for various values of ∆t,∆x for α=0.5 and x=0.5.

CONCLUSION

In this work, a numerical method based on the finite-difference approach is provided for solution of fractional-time diffusion equation, using the Neu-mann and Robin Boundary Conditions. The time derivative, meaning Caputo, is described of the α order. It has been concluded that by stability analysis of the problem, applied method is unconditionally stable, also the numerical analysis of the problem led us to conclusion, that numerical answers have very good compliance with exact answer of the problem. Moreover, this method seems simple to implement and also requires relatively little memory meanwhile benefiting from the proper performance, which can be mentioned as one of its advantages. The method can be easily extended to solve fractional PDEs with higher dimensions.

ACKNOWLEDGMENT

We all authors are glad to bring our work online after long process and frequently review of the literature, therefore we would like to express grati-tude and appreciation to Prof. Ihsanullah Saqib (HOD of Mathematics and Physics Department of Nangarhar University), who revised our work for scientific and grammatical errors and as well we are appreciating all the editors and referees for their valuable comments and suggestions. 


CONFLICTS OF INTEREST

This research is contributed by all authors and no potential conflict of interest to publish it. 

Article References:

  1. Çelik, C. and Duman, M., (2012). Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional deri-vative. J. of comput. Phys., 231 (4), pp.1743-1750. https://doi.org/10.1016/j.jcp.2011.11.008 
  2. Chen, W., Ye, L. and Sun, H., (2010). Frac-tional diffusion equations by the Kansa met-hod. Computers & Mathematics with Appli-cations, 59(5), pp.1614-1620. https://doi.org/10.1016/j.camwa.2009.08.004 
  3. Demir, A. and Bayrak, M.A., (2019). A new approach for the solution of space-time frac-tional order heat-like partial differential equ-ations by residual power series method. Com-mun. in Math. and Appl, 10(3), pp.585-597.
  4. Demir, A., Bayrak, M.A. and Ebru, O., (2020). New approaches for the solution of space-time fractional Schrödinger equation. Advances in Difference Equations, 2020(1). https://doi.org/10.1186/s13662-020-02581-5 
  5. Demir, A., Bayrak, M.A. and Ozbilge, E., (2018). An approximate solution of the time-fractional Fisher equation with small delay by residual power series method. Mathematical Problems in Engineering, 2018. https://doi.org/10.1155/2018/9471910  
  6. Demir, A., Bayrak, M.A. and Ozbilge, E., (2019). A new approach for the approximate analytical solution of space-time fractional dif-ferential equations by the homotopy analysis method. Advances in Mathematical Physics, 2019. https://doi.org/10.1155/2019/5602565  
  7. Ford, N.J., Xiao, J. and Yan, Y., (2011). A finite element method for time fractional partial differential equations. Fractional Cal-culus and Applied Analysis, 14(3), pp.454-474.  https://doi.org/10.2478/s13540-011-0028-2 
  8. Jacobs, B.A., (2016). High order compact finite-difference and laplace transform method for the solution of time fractional heat equ-ations with dirchlet and neumann boundary conditions. Numerical Methods for Partial Differential Equations, 32(4), pp.1184-1199. https://doi.org/10.1002/num.22046 
  9. Khader, M.M., (2011). On the numerical solutions for the fractional diffusion equa-tion. Communications in Nonlinear Sci. and Numerical Simulation, 16(6), pp.2535-2542. https://doi.org/10.1016/j.cnsns.2010.09.007 
  10. Kirane, M., Malik, S.A. and Al Gwaiz, M.A., (2013). An inverse source problem for a two dimensional time fractional diffusion equation with nonlocal boundary conditions. Mathe-matical Methods in the Applied Sciences, 36 (9), pp.1056-1069. https://doi.org/10.1002/mma.2661 
  11. Lin, Y. and Xu, C., (2007). Finite difference/ spectral approximations for the time-fractional diffusion equation. Journal of computational physics, 225(2), pp.1533-1552.
  12. Lin, Y. and Xu, C., (2007). Finite difference/ spectral approximations for the time-fractional diffusion equation. Journal of computational physics, 225(2), pp.1533-1552. https://doi.org/10.1016/j.jcp.2007.02.001 
  13. Liu, F., Zhuang, P. and Liu, Q., (2015). Nu-merical methods of fractional partial dif-ferential equations and applications.
  14. Magin, R., Feng, X. and Baleanu, D., (2009). Solving the fractional order Bloch equation. Concepts in Magnetic Resonance Part A: An Educational Journal, 34(1), pp.16-23. 
  15. Metzler, R. and Klafter, J., (2000). Boundary value problems for fractional diffusion equ-ations. Physica A: Statistical Mechanics and its Applications, 278(1-2), pp.107-125. https://doi.org/10.1016/S0378-4371(99)00503-8  
  16. Murio, D.A., (2008). Implicit finite-difference approximation for time fractional diffusion equations. Computers & Mathematics with Applications, 56(4), pp.1138-1145.
  17. Pathiranage D. (2021). Numerical investiga-tion of dropwise condensation on smooth pla-tes with different wettability, Int. J. Mat. Math. Sci., 3(3), 60-73. https://doi.org/10.34104/ijmms.021.060073 
  18. Roul, P. and Goura, V.P., (2020). A high order numerical method and its convergence for time-fractional fourth order partial differential equations. Appl. Math. and Com, 366, p.124-727. https://doi.org/10.1016/j.amc.2019.124727 
  19. Samko, S.G., Kilbas, A.A. and Marichev, O.I., (1993). Fractional integrals and derivatives, 1, Yverdon-les-Bains, Switzerland: Gordon and breach science publishers, Yverdon. 
  20. Sayevand, K., Yazdani, A. and Arjang, F., (2016). Cubic B-spline collocation method and its application for anomalous fractional dif-fusion equations in transport dynamic systems. J. of Vibr. and Control, 22(9), pp. 2173-2186. https://doi.org/10.1177%2F1077546316636282 
  21. Smith, G.D., Smith, G.D. and Smith, G.D.S., (1985). Numerical solution of partial differen-tial equations: finite difference methods. Ox-ford university press. 
  22. Sun, H., Chen, W. and Sze, K., (2013). A semi-discrete finite element method for a class of time-fractional diffusion equations. Philo-sophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sci-ences, 371(1990), p.20120268. https://doi.org/10.1098/rsta.2012.0268 
  23. Sun, Z.Z. and Wu, X., (2006). A fully discrete difference scheme for a diffusion-wave sys-tem. Appl. Numer. Math., 56(2), pp.193-209. https://doi.org/10.1016/j.apnum.2005.03.003 
  24. Sweilam, N.H., Khader, M.M. and Mahdy, A.M.S., (2012). Crank-Nicolson finite-differ-ence method for solving time-fractional dif-fusion equation. J. of Fractional Calculus and Applications, 2(2), pp.1-9. 
  25. Tamsir, M., Nigam, D. and Chauhan, A., (2021). Approximation of Caputo time-fractional dif-fusion equation using redefined cubic expo-nential B-spline collocation technique. AIMS Mathematics, 6(4), pp.3805-3820.
  26. Tian, W., Zhou, H. and Deng, W., (2015). A class of second order difference approxima-tions for solving space fractional diffusion equations. Mathematics of Computation, 84 (294), pp.1703-1727. https://doi.org/10.1090/S0025-5718-2015-02917-2  
  27. Usta, F. and Sarıkaya, M.Z., (2019). The ana-lytical solution of Van der Pol and Lienard differential equations within conformable frac-tional operator by retarded integral inequa-lities. Demon. Mathem., 52(1), pp. 204-212. https://doi.org/10.1515/dema-2019-0017 
  28. Usta, F., (2021). Numerical analysis of frac-tional Volterra integral equations via Bernstein approximation method. Journal of Computa-tional and Applied Mathematics, 384, p.113-198. https://doi.org/10.1016/j.cam.2020.113198 
  29. Xie, C. and Fang, S., (2019). A second order finite difference method for fractional dif-fusion equation with Dirichlet and fractional boundary conditions. Numerical Methods for Partial Differential Equations, 35(4), pp.1383-1395. https://doi.org/10.1002/num.22355 
  30. Yavuz, M., Usta, F. and Bulut, H., (2020). Analysis and numerical computations of the fractional regularized long-wave equation with damping term, 1–1. https://doi.org/10.1002/mma.6343 
  31. Yu D and Tan H., (2003). Numerical methods of differential equations. Beijing. Sci. Pub.
  32. Yuste, S.B. and Lindenberg, K., (2002). Sub diffusion-limited reactions. Chemical physics, 284(1-2), pp.169-180. https://doi.org/10.1016/S0301-0104(02)00546-3 
  33. Zhai, S. and Feng, X., (2016). A block-cen-tered finite-difference method for the time-fractional diffusion equation on non uniform grids. Numerical Heat Transfer, Part B: Fun-damentals, 69(3), pp.217-233.https://doi.org/10.1080/10407790.2015.1097101  
  34. Zhuang, P. and Liu, F., (2006). Implicit dif-ference approximation for the time fractional diffusion equation. Journal of Applied Mathe-matics and Computing, 22(3), pp.87-99. https://doi.org/10.1007/BF02832039 

Article Info:

Academic Editor

Dr. Toansakul Tony Santiboon, Professor, Curtin University of Technology, Bentley, Australia.

Received

August 22, 2022

Accepted

September 28, 2022

Published

October 4, 2022

Article DOI: 10.34104/ijmms.022.0950108

Corresponding author

Shafiullah Niazai*

Assistant Professor, Dept. of Mathematics, Laghman University, Mehtarlam city, Laghman, Afghanistan.

Cite this article

Niazai S, Rahimzai AA, Danesh M, and Safi B. (2022). Numerical solution of diffusion equation with caputo time fractional derivatives using finite-difference method with Neumann and Robin boundary conditions, Aust. J. Eng. Innov. Technol., 4(5), 95-108. https://doi.org/10.34104/ijmms.022.0950108 

Views
407
Download
587
Citations
Badge Img
Share