key: cord-0634157-9dw7dtl2
authors: Ngondiep, Eric
title: Unconditional Stability Of A Two-Step Fourth-Order Modified Explicit Euler/Crank-Nicolson Approach For Solving Time-Variable Fractional Mobile-Immobile Advection-Dispersion Equation
date: 2022-05-06
journal: nan
DOI: nan
sha: 155ddbcc4b64d87f2c8d7c238bd4da656d09481e
doc_id: 634157
cord_uid: 9dw7dtl2

This paper considers a two-step fourth-order modified explicit Euler/Crank-Nicolson numerical method for solving the time-variable fractional mobile-immobile advection-dispersion model subjects to suitable initial and boundary conditions. Both stability and error estimates of the new approach are deeply analyzed in the $L^{infty}(0,T;L^{2})$-norm. The theoretical studies show that the proposed technique is unconditionally stable with convergence of order $O(k+h^{4})$, where $h$ and $k$ are space step and time step, respectively. This result indicate that the two-step fourth-order formulation is more efficient than a broad range of numerical schemes widely studied in the literature for the considered problem. Numerical experiments are performed to verify the unconditional stability and convergence rate of the developed algorithm.

In the last twenty years, fractional calculus and fractional differential equations have found a broad range of applications in fluid flow, sound, electrodynamics, elasticity, biology, finances, geology, electrostatics, heat, hydrology and medical problems [2, 17, 6, 16, 8, 46, 43, 42] . A large class of complex models are deeply described via variable-order derivatives. Most recently, the time variable fractional order telegraph equation has been shown to be a suitable model for various physical phenomena. Since the equations modeled by the time fractional partial differential equations (FPDEs) are highly complex, there is no method that can compute an exact solution. The big challenge with such a set of equations is the development of fast and efficient numerical approaches in an approximate solution. In the literature, abundant numerical schemes have been proposed for solving time variable (or constant) order FPDEs such as: finite difference schemes of first order accuracy in time and spatial second-order convergence, compact finite difference scheme with convergence order O(τ + h 4 ) and some numerical procedures of higher order O(τ 2−γ + h 4 ). All these techniques were one-step methods [4, 48, 50, 9, 41, 49, 3] . The author [21] developed a twostep numerical scheme with convergence order O(τ 2− γ 2 + h 4 ) for solving the time-fractional convectiondiffusion-reaction equation with constant order derivative. For classical integer order ordinary/partial differential equations such as: Navier-Stokes equations, systems of ODEs, mixed Stokes-Darcy model, Shallow water problem, convection-diffusion-reaction equation, advection-diffusion model and conduction equation [18, 5, 20, 24, 39, 26, 29, 37, 28, 27, 14, 32, 23, 50] , a wide set of numerical techniques have been developed and deeply analyzed. For more details, we refer the readers to [30, 38, 19, 45, 33, 34, 44, 22] and references therein. This class of partial differential equations (PDEs) also have a large set of applications. Furthermore, some methods developed for solving inter order PDEs can be used to efficiently compute approximate solutions of time FPDEs with low computational costs [49, 9, 4] . In this work, we develop a two-step modified explicit Euler/Crank-Nicolson approach in an approximate solution of time-variable fractional mobile-immobile advection-dispersion model with subjects to suitable initial and boundary conditions. The proposed technique is unconditionally stable, convergence with order O(τ + h 4 ), fast and more efficient than a large class of numerical schemes widely studied in the literature for the considered problem [48, 13, 15, 36, 11, 10] . A time variable fractional mobile-immobile advection-dispersion equation describes a broad range of problems in physical or mathematical systems which include ocean acoustic propagation and heat conduction through a solid [1] . In such an equation, the integer order time derivative term is added to describing the motion of particles conveniently [1] . Furthermore, this equation lies in a class of second order PDEs that govern continuous time random walks with heavy tailed random waiting times.

In this paper, we propose a two-step fourth-order modified explicit Euler/Crank-Nicolson formulation for the time-variable fractional mobile-immobile advection-dispersion equation [48] and describing by the following initial-boundary value problem

with initial condition u(x, 0) = u 0 (x), on [L 0 , L],

and boundary condition u(L 0 , t) = g 1 (t) and u(L, t) = g 2 (t), on [0, T ],

where u t , u x , and u 2x denote ∂u ∂t , ∂u ∂x , and ∂ 2 u ∂x 2 , respectively. f = f (x, t) represents the source term, u 0 is the initial condition whereas g 1 and g 2 designate the boundary conditions. cD β(x,t) 0t u, where (0 < β 1 ≤ β(x, t) ≤ β 2 < 1), is denotes the variable-order fractional derivative which has been introduced in various physical fields and defined (for example, in [47] ) as

For the sake of discretization and error estimates, we assume that the analytical solution of the initialboundary value problem (1)-(3) is regular enough. We remind that the goal of this work is to develop an efficient numerical method for solving the time-variable fractional equation (1) subjects to initial condition (2) and boundary condition (3) . Specifically, the attention is focused on the following three items:

(i1) full description of a two-step fourth-order modified explicit Euler/Crank-Nicolson approach for timevariable FPDE (1) with initial-boundary conditions given by (2) and (3), respectively, (i2) analysis of the unconditional stability and convergence rate of the new technique, (i3) a broad range of numerical evidences which confirm the theoretical results.

The outline of the paper is as follows. In Section 2, we construct the two-step fourth-order modified explicit Euler/Crank-Nicolson numerical scheme for computing an approximate solution to the initial-boundary value problem (1)- (3) . Section 3 analyzes both stability and error estimates of the new approach using the L ∞ (0, T ; L 2 )-norm. Some numerical examples that confirm the theoretical study are presented and discussed in Section 4. Finally, in Section 5 we draw the general conclusions and provide our future investigations.

2 Development of the two-step modified explicit Euler/crank-Nicolson numerical approach

In this section we develop a two-step fourth-order modified explicit Euler/crank-Nicolson numerical method for solving the initial-boundary value problem (1)- (3) . Starting with the explicit Euler scheme, the new approach is a two-step implicit method which approximates the time-variable fractional operator (4) using both forward and backward difference approximations in each step whereas the convection and diffusion terms are approximated using the central difference formulation. To compute a numerical solution of problem (1)-(3), we introduce a uniform grid of mesh points (x j , t n ), where x j = L 0 + jh, j = 0, 1, 2, ..., M , and t n = nk, for n = 0, 1, 2, ..., N . In this discretization, M and N are two positive integers, h = L−L0 M and k T N are the space step and the time step, respectively. The space of grid functions is defined as

For the convenience of writing, we set u(x j , t n ) = u n j . The exact solution and the computed one at the grid point (x j , t n ) are denoted by u n j and U n j , respectively. Furthermore, the values of the functions β and f at the mesh point (x j , t n ) are given by β n j and f n j , respectively. Finally we introduce the positive parameter α ∈ (0, 1 2 ).

Furthermore, we define the following operators

and we introduce the given norms

together with the inner products

where D = [L 0 , L] × [0, T ], | · | is the C-norm. The spaces L 2 (L 0 , L), L ∞ (0, T ; L 2 ) and C 6,3 D are equipped with the norms · 2 , | · | ∞,2 and | · | C 6,3 D , respectively, whereas the Hilbert space L 2 (L 0 , L) is endowed with the scalar product (·, ·).

The application of the Taylor series for u at the grid point (x j , t n+ 1 2 +α ) with time step k 2 using forward difference representation gives

Utilizing equation (1), this becomes

Expanding the Taylor series for u at the mesh point (x j , t n+ 1 2 +α ) with time step k 2 using both forward and backward difference formulations to get

Using equation (1), we obtain u n+1+α j = u

Subtracting the second equation from the first one and using the Mean value theorem, it is easy to see that

which is equivalent to

For the convenience of writing, we set γ

So the variable-order fractional time derivative given by (4) at the grid point (x j , t n+ 1 2 +α ) can be rewritten as

Let q 1,uj 2,j be the polynomial of degree two approximating u j (t) at the mesh points (t l , u l j ), (t l+ 1 2 , u 

where c∆ 

The summation index in equation (12) varies in the range: l = 0, 1 2 , 1, 3 2 , 2, ..., n, whereas the summation index in (13) satisfies: i = 0, 1, 2, 3, ..., n − 1. In relation (12) , the coefficients a α,β 

a α,β

for n ≥ 1 and i = 0, 1, 2, ..., n − 1. Furthermore,

where c∆

P 0,uj 1,s is the first-order polynomial interpolating u j at the mesh point (u 0 j , t 0 ) and (u α j , t α ).

In a similar manner, setting γ n+1+α j = Γ(1 − β n+1+α j ) −1 and replacing λ by β n+1+α j in the formulas obtained in [21] , pages: 8, 12, and 17, results in

where

with a α,β 1+α

and for n ≥ 1

Here the terms f α,β n+1+α j n+1,r and d α,β n+1+α j n+1,r , for n ≥ 1 and r = 0, 1, 2, ..., n + 1, are obtained by replacing λ by β n+1+α j in [21] , page 8-9, by

We recall that in equation (21) the summation index varies in the range l = 0, 1 2 , 1, 3 2 , 2, ..., n+ 1 2 . Furthermore, replacing in relation (52) given in [21] , λ with β n+1+α j , we obtain

where E 2,uj s,j (s) is the error associated with the quadratic polynomial interpolating u j (t) at the grid points (u l j , t l ), (u , for s = 0, 1 2 and 1. Using the Taylor series expansion for u about the grid point (x j , t n+s+α ) with step size h using both forward and backward representations, the author [21] has shown that u n+s+α

u n+s+α

for s ∈ {0, 1 2 , 1}, where

For s = 0, replacing n + 1 by n in (20) and plugging the obtained equation with (27) , (26) and (8) yields

Replacing n + 1 by n in (21), substituting the new equation into (30) and rearranging terms results in

Since γ n+α

Setting u 0x (x, t) = u(x, t) and applying the Taylor series expansion for the functions u mx (m = 0, 1, 2) at the mesh points (x j , t n ), (x j , t n+ 1 2 +α ) and (x j , t n+1+α ) with time step k 2 gives

Combining (33) and (26)- (27) to get

where we set δ 4 0x u n+α

and δ 4 0x u n j = u n j . For m = 0, 1, 2, substituting the last equation of (34) into (32) and performing simple computations to obtain

For n = 0, utilizing relation (19) and substituting (17) into (8) to get

Since

Replacing n and s by 0 into (26)- (27) and (33)- (34) and m by 0, 1, 2, into the first equations of (33) and (34) , combining the obtained equations together with (37) and utilizing equality δ α t u 0

where

Suppose U n = (U n 2 , U n 3 , ..., U n M−2 ) be the approximate solution vector at time level n and u n = (u n 2 , u n 3 , ..., u n M−2 ) be the analytical one at time t n . Truncating the error terms in both equations (35) and (38), we obtain the first-step of the new approach

In addition, setting 

where θ n+s+α

s ∈ {0, 2 −1 , 1}. Thus, equations (40) and (43) can be expressed in the matrix form as

where "*" denotes the componentwise usual multiplication between two vectors, A 0 and A 1 are two (M − 2) × (M − 2) "pentadiagonal" matrices defined by

where

To complete the full description of the desired algorithm, we should develop the second-step of the method. Combining equations (9), (11) and (20), direct calculations give

Since for s = 1 2 , 1, γ n+s+α

and m = 0, 1, 2, plugging equations (12), (21) , (33) , (34) and (48), straightforward computations result in

which is equivalent to

where θ n+s+α

is defined by (43) . Omitting the error terms k

equation (49) can be approximated as

We introduce the following "pentadiagonal" matrices A and

where

A combination of equations (42) , (43) and (50) provides the following matrix form

where "*" represents the componentwise usual multiplication between two vectors. We recall that the summation index "l" varies in the range l = 0, 1 2 , 1, 3 2 , 2, ..., n, n + 1 2 . Furthermore, equation (50) denotes the second-step of the proposed two-step fourth-order modified explicit Euler/Crank-Nicolson numerical scheme in a computed solution of the initial-boundary value problem (1)-(3).

An assembly of equations (44) , (45) and (53) provides the new algorithm for solving the problem (1)-(3), that is, for n = 1, 2, ..., N − 1,

with initial and boundary conditions U 0 j = u 0 j , j = 0, 1, ..., M ; U n 0 = g n 1 and U n M = g n 2 , for n = 0, 1, ..., N,

where the matrices A 0 , A 1 , A and A 2 are given by relations (46)- (47) and (51)-(52). To start the algorithm we should set U n 1 = U n 0 and U n M−1 = U n M , for n = 0, 1, ..., N . However, the terms U n 1 and U n M−1 can be obtained by using any one-step fractional approach such as the method analyzed in [12] .

It's worth noticing that the coefficients of the pentadiagonal matrices A and A i , for i = 0, 1, 2, come from the entries of a g-Toeplitz matrix (T n,g (f )) or g-circulant matrix (C n,g (f )), generated by a Lebesgue integrable function f defined over the domain (−π, π), where g is a nonnegative integer. Specifically, these matrices are called band Toeplitz matrices which represent a subclass of g-Toeplitz structures. For more details about g-Toeplitz, g-circulant and band Toeplitz matrices, we refer the readers to [25, 7, 31, 35] and references therein. Furthermore, since the matrices A and A 0 are not symmetric, at time level n or n+ 1 2 , each system of linear equations (54)-(57) can be efficiently solved using the Preconditioned Generalized Minimal Residual Algorithm [40] . where s = 1 2 , 1, β(x, t) = β is constant, cD λ 0t u n+s+α and c∆ λ 0t u i+ 1 2 +α , for 0 ≤ n ≤ N − 1, are defined by (11), (12) and (17)-(21), and C s for s ∈ { 1 2 , 1} are positive constants independent of the mesh size h and time step k.

Proof. Since 0 < β < 1 is constant and 0 < α < 1 2 , the proof of this Lemma can be found in [21] . 

Furthermore,

.., n, (resp., n + 1 2 ).

Lemma 3.3. Let (a α,β n+s,l ) l≤n+s be the generalized sequences defined by relations (14)- (15) and (22)- (23) . For every mesh function u(·, ·) defined on the grid space U hk = {u n j , 0 ≤ j ≤ M and n = 0, 1, 2, ..., N }, setting W β,l j = l r=0 a α,β n+s,r+ 1 2 δ t u r j , the following estimates hold for s = 1 2 , 1,

Furthermore,

Proof. The proof of (62) is obtained by replacing λ with β in the proof of Lemma 3.3 established in [21] . The proof of (63) is obvious since a α,β n+s,l < a α,β n+s,l+ 1 2 , for s = 1 2 , 1, and l = 1 2 , 1, 3 2 , 2, ..., n (resp., n + 1 2 ). In addition, it is easy to see that a α,β

Lemma 3.4. Given (a α,β n+s,l ) l be the generalized sequences defined by equations (14)- (15) and (22)-(23), for any grid function v(·, ·) defined on the grid space U hk , it holds m l=l0 a α,β

for m ∈ {n, n + 1 2 } and l = l 0 , l 0 + 1 2 , l 0 + 1, l 0 + 3 2 , ..., m, where l 0 is a nonnegative integer satisfying l 0 ≤ m. Proof. Expanding the left side of this equality and rearranging terms to obtain the result. 

for j = 2, 3, ..., M − 2, where q is any nonnegative rational number. If

where C p is a positive constant independent of the time step k and space step h.

Proof. In equations (28) and (29), replacing n + s + α by q to obtain

Since δ

Performing straightforward computations, it is not difficult to show that

Utilizing assumption v q 1 = 0 and v q M−1 = 0, this becomes

The last equality follows from the assumption v q 1 = 0 and v q M−1 = 0. Analogously, one easily shows that

Plugging equations (67) and (69)-(71) yields

Multiplying both sides of this equation by h, applying the Hölder and Cauchy-Schwarz inequalities, using the definition of L 2 -norm and the scalar product (·, ·), this results in

In a similar manner, one easily proves that

Combining equations (68) and (73), it is not hard to observe that

Multiplying this equation by h, utilizing the Hölder and Cauchy-Schwarz inequalities together with the definitions of L 2 -norm and scalar product (·, ·) to get

A combination of (65) and estimates (73)-(74) gives

This completes the proof of the first estimate in (66). The proof of the second estimate in (66) is obtained thanks to the Poincaré-Friedrich inequality.

Lemma 3.6. Suppose u ∈ C 2,0 D , be a function defined on D = [L 0 , L]×[0, T ], satisfying u(L 0 , t) = u(L, t) = 0, for any t ∈ [0, T ]. Let U (t) = (U 0 (t), U 1 (t), ..., U M (t)) be a grid function such that, U j (t) = u(x j , t), for j = 0, 1, ..., M . So, it holds

for every t ∈ [0, T ].

Proof. We should show that the operator −Lu = −u 2x + u x satisfies: (−Lu(t), u(t)) ≥ δ x u(t) 2 2 and then, use the discrete L 2 -norm defined in relation (6) to conclude. The application of the Taylor series expansion using backward difference formulation gives:

. Using the conditions u 1 (t) = u M−1 (t) = 0, for every t ∈ [0, T ] together with the summation by parts and the equality a(a − b)

, for any real numbers a and b, direct computations yield

The last equality follows from hδ x u M− 3 2 (t) = −u M−2 (t) and δ x u 3

Substituting approximation (77) into relation (76) to obtain

Furthermore, since u j (t) = U j (t), for j = 0, 1, ..., M , for h sufficiently small, neglecting the infinitesimal terms O(h 2 ) and O(h), equation (76) implies

This ends the proof of the first estimate in Lemma 3.6. Since

, the proof of Lemma 3.6 is completed thanks to the first estimate in (75) and the definition of the scalar product given by (7).

Armed with Lemmas 3.1-3.6, we should state and prove the main result of this work (Theorem 3.1).

Theorem 3.1. (Unconditional stability and Error estimates). Suppose U be the approximate solution provided by the proposed approach (54)-(58) and let u be the analytical solution of the initial-boundary value problem (1)- (3) . let β(x, t) = β ∈ (0, 2 3 ), for any (x, t) ∈ D, be a positive constant function, 0 < α < 1 2 be a parameter and let (a α,β ·,l ) l be the generalized sequences defined by relations (14)- (15) and (22)- (23) . Thus, the following estimates are satisfied

Furthermore, denote e = u − U be the error term, it holds

where C is a positive constants independent on the space size h and time step k.

We recall that estimate (78) suggests that the proposed technique (54)-(58) is unconditionally stable whereas inequality (79) shows that the developed numerical scheme is fourth-order spatial convergent and temporal accurate of order O(k).

Proof. Let e n+ 1 2 = u n+ 1 2 − U n+ 1 2 be the temporary error term and e n+1 = u n+1 − U n+1 be the exact one at time level n + 1. Subtracting equation (41) from (35) and utilizing (43) and (65) 

which is equivalent to

Multiplying both sides of this equation by e 

.., O(k 2 + k 3 + kh 4 )), multiplying both sides of this estimate by 2h, summing the obtained estimate up from j = 2, 3, ..., M − 2, and using the definition of the L 2 -norm and scalar product given by relations (6) and (7), respectively, this provides 2 , where C p denotes a positive constant independent of k and h. Using this, estimate (66) and the Hölder inequality, straightforward calculations give θ n+α l+1 − θ n+α

Estimate (88) is satisfied for any 0 < α < 2 −1 . For n ≥ 1, it comes from (15)- (16) and (23)- (24) , that α 1−β = f α,β n,n < a α,β n,n = f α,β n,n−1 + f α,β n,n = f α,β n+ 1 2 ,n−1 + f α,β n+ 1 2 ,n = a α,β

Since θ n+α θ n+α l+1 − θ n+α

where the summation index "l" varies in the range: 1 2 , 1, 3 2 , 2, ..., n − 1 2 , n. (90) is equivalent to

for every α ∈ (0, 2 −1 ). Setting C 1 = lim 

where C 2 = C 1 (1 + C 2 1 2 ).

In a similar way, a combination of equations (49) , (50) and (65) gives θ n+1+α

Applying the summation by parts and rearranging terms, this becomes 

Performing direct computations, using the Hölder and Poincaré-Friedrich inequalities, equations (11) , (20) and the second estimate in (66), it is not difficult to show that

where C 2 > 0, is a constant independent of the time step k and mesh size h. 

For small values of k, 1 − 1 2 θ n+ 1 2 +α n+1 > 0. Utilizing this and (59), relation (101) becomes

Letting C 3,α = max

It follows from (89) that a α,β n+ 1 2 ,n+ 1 2 = a α,β n+1,n+1 . This fact, together with (43) and (81) give θ n+1+α n+1 = θ n+ 1 2 +α n+1 . Since, 1 + 1 2 θ n+1+α n+1 −1 < 1, multiplying both sides of (102) by 1 + 1 2 θ n+1+α n+1 −1 to get e n+1 2 2 ≤ max

Taking the limit when α approaches 1 4 , (103) provides e n+1 2 2 ≤ max

It is worth noticing to remind that the summation index "l" varies in the range: l = 0, 1 2 , 1, ..., n + 1 2 , n + 1. Now, setting Z q = max 0≤l≤q e l 2 2 and C 4 = max{ C 2 , C 3 },

estimates (92) and (104) can be rewritten as

Substituting the first estimate into the second one gives

Summing this up from n = 0, 1, 2, ..., N − 1, to obtain

But, it comes from the initial condition (58) that e 0 j = u 0 j − U 0 j = 0, for j = 0, 1, 2, ..., M . Furthermore, since k = T N so, N k = T . These facts together with (105) and (106) yield max 0≤l≤N e l 2 2 ≤ e 0 2 2 + 2

This is equivalent to

Taking the square root to obtain

These estimates imply

for n = 0, 1, 2, ..., N − 1 (resp., N ). But | u q 2 − U q 2 | ≤ u q − U q 2 = e q 2 , thus

for n = 0, 1, 2, ..., N − 1 (resp., N ), which imply

This ends the proof of estimate (78) in Theorem 3.1. Now, since k < 1 and 2 − β > 1, so k 2−β , k 2 ≤ k, thus k + k 2−β + k 2 + h 4 ≤ 3(k + h 4 ). Using this, (107) implies

This completes the proof of Theorem 3.1.

This section considers some computational results to show the unconditional stability and the convergence order of the new approach (54)-(58) applied to time-variable fractional mobile-immobile advection-dispersion equation (1) subjects to initial and boundary value conditions (2) and (3), respectively. To demonstrate the efficiency and accuracy of the proposed algorithm, two examples are taken in [48] . We set k = h 4 , where h ∈ {2 −i , i = 1, 2, 3, 4} so, k = 2 −4r , r = 1, 2, 3, 4, and we compute the L ∞ -norm of the numerical solution: U n , the exact one: u n , and the corresponding error: e n , at time level n, using the following formulas

In each example, the numerical evidences are performed with two different order functions: β 1 (x, t) = 1 − 2 −1 e −xt and β 2 (x, t) = 4 3 − 5.10 −3 cos(xt) sin(xt). Furthermore, the convergence rate R(k, h) of the new algorithm is estimate using the formula

where we set k = h 4 . Finally, the numerical computations are carried out by the use of MATLAB R2013b. Figures 1-4 suggest that the proposed two-step technique (54)-(58) is unconditionally stable whereas Tables 1-4 indicate that the developed numerical method is temporal first-order accurate and spatial fourthorder convergent. These numerical studies confirm the theoretical results provided in Section 3, Theorem 3.1.

• Example 1. Let D = [0, 1] × [0, 1] be the domain. The parameter α = 0.25, 0.49 and the function β is defined as β(x, t) = 1 − 2 −1 e −xt . Consider the following time-variable fractional mobile-immobile defined in [48] by

. The analytical solution u is given by 1] . We assume that the parameters α ∈ {0.25, 0.49} and the function β is given by β(x, t) = 4 3 − 5.10 −3 cos(xt) sin(xt). We consider the following time-variable fractional mobile-immobile advection-dispersion model defined in [48] as

u(x, 0) = u 0 (x) = 5 sin(πx) for 0 ≤ x ≤ 1,

where f (x, t) = 5(1 + π 2 (1 + t)) sin(πx) + 5 sin(πx)t 1−β(x,t) Γ(2−β(x,t)) − 5π(1 + t) cos(πx). The exact solution u is defined as u(x, t) = 5(1 + t)x 2 sin(πx). We observe from this table that the proposed method is temporal second order convergent and spatial fourth order accurate.

This paper has developed a two-step fourth-order modified explicit Euler/Crank-Nicolson formulation for solving the time-variable fractional mobile-immobile advection-dispersion model subjects to suitable initial and boundary value conditions. Both stability and error estimates of the new technique have been deeply analyzed in L ∞ (0, T ; L 2 )-norm. The theory has shown that the proposed approach is unconditionally stable, first-order convergence in time and fourth-order accurate in space (see Theorem 3.1). This theoretical analysis is confirmed by two numerical examples. Especially, the graphs (Figures 1-4) show that the new procedure is both unconditionally stable and convergent whereas Tables 1-4 indicate the convergence rate of the developed algorithm (convergence with order O(k) in time and fourth-order accurate in space). Furthermore, the theoretical and numerical studies suggested that the proposed numerical scheme (54)-(58) is faster and efficient than a wide set of numerical methods [48, 13, 15, 36, 11, 10] developed for the considered problem (1)-(3). Solving the two-dimensional time-variable fractional problems by the use of a two-step fourth-order explicit/implicit scheme will be the topic of our future work. In addition, we will also be interested in the analysis of the Preconditioned Generalized Minimal Residual Method in an efficient solution of linear systems of equations (54)-(58). Specifically, in view to construct efficient preconditioners for such systems of linear equations, our study will consider the spectral behavior of the sequences of coefficient matrices 

A simple and efficient random walk solution of multi-rate mobile/immobile mass transport equations

Fractional calculus in hydrologic modeling: a numerical perspective

A high-order compact exponential scheme for the fractional convection-diffusion equation

Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation

A meshless local Petrov-Galerkin method for the time-dependent Maxwell equations

Application of collocation method for solving a parabolic-hyperbolic free boundary problem which models the growth of tumor with drug application

A note on the (regularizing) preconditioning of g-Toeplitz sequences via g-circulants

Novel numerical analysis of multi-term time fractional viscoelastic non-Newtonian fluid models for simulating unsteady and MHD and Couette flow of a generalized Oldroyd-B fluid

A compact finite difference scheme for the fractional sub-diffusion equations

Efficient multistep methods for tempered fractional calculus: Algorithms and simulations

A RBF meshless approach for modeling a fractal mobile/immobile transport model

Finite difference/finite element method for a nonlinear time-fractional fourth-order reaction-diffusion problem

A RBF meshless approach for modeling a fractal mobile/immobile transport model

Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term

Numerical methods and analysis for a class of fractional advectiondispersion models

Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time

Fractional green function for linear time-fractional inhomogeneous partial differential equations in fluid mechanics

Mathematical models of complete shallow water equations with source terms, stability analysis of Lax-Wendroff scheme

Unconditional stability over long time intervals of a two-level coupled MacCormack/Crank-Nicolson method for evolutionary mixed Stokes-Darcy model

Stability analysis of MacCormack rapid solver method for evolutionary Stokes-Darcy problem

A two-level fourth-order approach for time-fractional convection-diffusion-reaction equation with variable coefficients

A fourth-order two-level factored implicit scheme for solving two-dimensional unsteady transport equation with time dependent dispersion coefficients

A novel three-level time-split MacCormack scheme for two-dimensional evolutionary linear convection-diffusion-reaction equation with source term

A novel three-level time-split approach for solving two-dimensional nonlinear unsteady convection-diffusion-reaction equation

Spectral features and asymptotic properties of gcirculant and g-Toeplitz sequences

Long time stability and convergence rate of MacCormack rapid solver method for nonstationary Stokes-Darcy problem

An efficient three-level explicit time-split approach for solving 2D heat conduction equations

An efficient three-level explicit time-split scheme for solving two-dimensional unsteady nonlinear coupled Burgers equations

A robust three-level time-split MacCormack scheme for solving two-dimensional unsteady convection-diffusion equation

Long time unconditional stability of a two-level hybrid method for nonstationary incompressible Navier-Stokes equations

How to determine the eigenvalues of g-circulant matrices

A two-level factored Crank-Nicolson method for two-dimensional nonstationary advection-diffusion equation with time dependent dispersion coefficients and source/sink term

A three-level time-split MacCormack method for two-dimensional nonlinear reaction-diffusion equations

A robust numerical two-level second-order explicit approach to predict the spread of covid-2019 pandemic with undetected infectious cases

Spectral distribution in the eigenvalues sequence of product of g-Toeplitz structures

A high-order predictor-corrector method for solving nonlinear differential equations of fractional order

Numerical solution of fractional advection-diffusion equation with a nonlinear source term

Fourth-order exponential finite difference methods for boundary value problems of convective diffusion type

A meshless method to the numerical solution of an inverse reactiondiffusion-convection problem

GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems

The Sinc-Legendre collocation method for a class of fractional convection-diffusion equations with variable coefficients

Fractal mobile/immobile solute transport

High-order numerical scheme for the fractional heat equation with Dirichlet and Neumann boundary conditions

A fourth-order compact ADI method for solving two-dimensional unsteady convection-diffusion problem

A high-order exponential scheme for solving 1D unsteady convection-diffusion equations

A fractional model to describe the Brownian motion of particles and its analytical solution

Second-order approximations for variable order fractional derivatives: algorithms and applications

A novel numerical method for the timevariable fractional order mobile-immobile advection-dispersion model

Compact alternating direction implicit scheme for the twodimensional fractional diffusion-wave equation

New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation

Acknowledgment. This work has been partially supported by the Deanship of Scientific Research of Imam Mohammad Ibn Saud Islamic University (IMSIU) under the Grant No. 331203.