images
Computer Modeling in Engineering & Sciences
images

DOI: 10.32604/cmes.2023.044878

ARTICLE

A Novel Method for Linear Systems of Fractional Ordinary Differential Equations with Applications to Time-Fractional PDEs

Sergiy Reutskiy1, Yuhui Zhang2,*, Jun Lu3,* and Ciren Pubu4

1A. Pidhornyi Institute of Mechanical Engineering Problems of NAS of Ukraine, Kharkiv, 61046, Ukraine
2College of Mechanics and Materials, Hohai University, Nanjing, 211100, China
3Nanjing Hydraulic Research Institute, Nanjing, 210029, China
4Shigatse Everest Urban Investment and Development Group Co, Ltd., 857000, China
*Corresponding Authors: Yuhui Zhang. Email: zyh940511@outlook.com; Jun Lu. Email: lujunhhu1981@126.com
Received: 10 August 2023; Accepted: 17 November 2023; Published: 29 January 2024

Abstract: This paper presents an efficient numerical technique for solving multi-term linear systems of fractional ordinary differential equations (FODEs) which have been widely used in modeling various phenomena in engineering and science. An approximate solution of the system is sought in the form of the finite series over the Müntz polynomials. By using the collocation procedure in the time interval, one gets the linear algebraic system for the coefficient of the expansion which can be easily solved numerically by a standard procedure. This technique also serves as the basis for solving the time-fractional partial differential equations (PDEs). The modified radial basis functions are used for spatial approximation of the solution. The collocation in the solution domain transforms the equation into a system of fractional ordinary differential equations similar to the one mentioned above. Several examples have verified the performance of the proposed novel technique with high accuracy and efficiency.

Keywords: System of FODEs; numerical solution; Müntz polynomial basis; time fractional PDE; BSM collocation method.

1  Introduction

A novel numerical method for solving a linear system of fractional ordinary differential equations (FODEs).

A^Dt(α)[Υ(t)]+∑k=1KA^k(t)Dt(αk)[Υ(t)]=B^(t)Υ(t)+F(t),0≤t≤T,(1)

is proposed in this paper. Here: n−1<α≤n; 0≤αk<α; positive integer n defines the maximal order of the time derivative in the equation and so, defines the number of the initial conditions of the problem; Υ(t)=[Υ1(t),Υ2(t),….,ΥN(t)]T is the N−vector of unknowns; A^ is a constant non-singular N×N matrix and A^k(t), B^(t) are time-dependent N×N matrices; F(t)=[f1(t),f2(t),….,fN(t)]T. The operator Dt(ν) denotes the Caputo fractional derivative defined by [1,2]:

Dt(α)[f(t)]={1Γ(m−α)∫0tdt(m)f(τ)dτ(t−τ)α−m+1,m−1<α<m,dt(m)f(t),α=m,(2)

where m∈𝒩={1,2,….}, and Γ(z) denotes the gamma function. In particular, for the power functions we get:

Dt(α)[tz]={0,if z∈𝒩0 and z<n,Γ(z+1)Γ(z+1−α)tz−α,if z∈𝒩0 and z≥n or z∉𝒩0 and z>n−1,(3)

where 𝒩0={0,1,2,….}.

The equations similar to (1) often arise in the modeling of various physical phenomena such as the models of pollution in systems of lakes [3–5], of processing the Magnetic Resonance Imaging (MRI) data [6], of the spread of infections [7,8], and also in modeling the nuclear magnetic resonance [9,10]. Recently such problems have become very relevant due to the widespread use of the fractional-order mathematical model of the COVID-19 disease [11–14].

Besides, as is shown below, based on this technique an effective method for solving multi-term time fractional partial differential equations (TFPDEs) of the type

Dt(α)[v(x,t)]+∑k=1Kak(t)Dt(αk)[v(x,t)]=b(t)L(x)[v(x,t)]+f(x,t),(4)

can be developed. Here ai(t) and b(t) are smooth enough functions and

L(x)[v(x,t)]=∂∂x(c(x)∂v(x,t)∂x)=c(x)∂2v(x,t)∂x2+∂c(x)∂x∂v(x,t)∂x(5)

is a spatial differential operator of the second order defined for 0≤x≤1. The function c(x) has the physical sense of the diffusivity as the transport problem is considered. Let us note that Eq. (4) includes many different known equations as particular cases. For example:

– the time-fractional sub-diffusion equation [15]

Dt(α)u(x,t)=a∂2u(x,t)∂x2+f(x,t),0<α<1,a>0,(6)

– the time-fractional telegraph equation [16–18]

Dt(α)u(x,t)+α1Dt(α−1)u(x,t)+α2u(x,t)=∂2u(x,t)∂x2+f(x,t),(7)

– the multi-term time-fractional diffusion and diffusion-wave equations [19–21]

Dt(μ)u(x,t)+∑i=1nαiDt(νi)u(x,t)=Ke∂2u(x,t)∂x2+f(x,t),0<νi<μ<2,αi,Ke>0,(8)

– the time-fractional modified anomalous sub-diffusion equation [22–26]

∂u(x,t)∂t=[α1Dt(1−v1)+α2Dt(1−v2)]∂2u(x,t)∂x2+f(x,t),0<v1,v2<1,α1,α2>0.(9)

It has been shown by many researchers that the fractional equations are more suitable for modeling some real-world applications compared with the equations of the integral order. The reviews of some real-world applications of the fractional equations were provided by Almeida et al. [27] and Sun et al. [28], in physics [29], solid mechanics [30], and fluid mechanics [31]. The application of fractional equations can also be noted in the recently published books [32,33] and we refer readers to them and to the references therein.

The exact solutions of the fractional equations are critically important for revealing complex physical phenomena. Some well-known analytical methods have been proposed for this goal: the Laplace transform method [34,35], the Green function method [36], the Fourier transform method [37], the variational iteration method [38], the Adomian decomposition method [39], the method of separating variables [40], etc.

However, because analytical solutions are available only for a narrow class of fractional problems, a great number of numerical techniques have been developed. Currently, the finite difference (FD) and finite element (FE) techniques are still the most useful tools in this field.

A survey of the FD methods for solving FODEs and fractional PDEs was presented by Li et al. in [41]. Some non-standard FD techniques were proposed to solve complex fractional systems in [42–44]. The fast FD methods for the fractional equations were proposed for solving 2D/3D space-fractional diffusion equations in [45,46]. Similar fast FD techniques were proposed for distributed-order space-fractional problems in [47], for parameters identifying problems governed by fractional equations in [48], for time-dependent space-fractional diffusion equations with fractional boundary conditions in [49], for the nonlinear fractional wave equation in [50], for fractional equations with singularity in [51], etc. The FE techniques also are the most commonly used for solving fractional equations. The FE approach was used to solve 1D fractional equations in [52,53] and for 2D fractional equations in [54–56]. Many works focus on the error analysis of the FE methods such as [57–59].

Recently meshless methods have become the focused issues of the researchers in science and engineering. The meshless methods can be divided into two groups: the pure collocation techniques [60–63] and the methods based on the integration [64–67]. To improve the accuracy of the meshless methods combinations with semi-analytical techniques have been proposed. The Laplace transform method has been coupled with the Adomian decomposition method in [68]. The analytical and semi-analytical solutions of the time-fractional Cahn–Allen equation have been studied by Khater et al. in [69]. A semi-analytical solution for the time-fractional diffusion equation has been developed by Kazem et al. in [70]. The homotopy analysis transform method [71] and the fundamental solution method [72] belong to the same group of techniques. Five semi-analytical techniques for solving the fractional nonlinear telegraph equation have been studied in [73].

In this paper, a new semi-analytical meshless technique-the backward substitution method (BSM) [74,75] is proposed to solve multi-term linear systems of FODEs. Based on the method provided a flexible and efficient numerical technique is constructed to solve the TFPDE (4). Applying the collocation approach, the original TFPDE is transformed into the system of FODEs which can be handled by the proposed new technique. The performance of this approach has been thoroughly examined by typical numerical examples. The test results are compared with the exact solutions and with the data obtained by other numerical techniques.

The rest of the paper is organized as follows. The detailed scheme of the BSM for solving the system of FODEs is formulated in Section 2. The scheme for solving the TFPDEs is presented in Section 3. The numerical examples are given in Section 4. Finally, some conclusions are briefly discussed in Section 5.

2  Backward Substitution Method for FODEs

In this section, we propose a novel numerical scheme for solving the system (1) subjected to the initial conditions (ICs):

Υ(0)=Υ0, ∂tΥ(0)=Υ1… . ,∂t(n−1)Υ(0)=Υn−1.(10)

Let us define a new vector-function P(t) which satisfies the relation

P(t)=Υ(t)−Θn−1(t),(11)

where

Θn−1(t)=Υ0+Υ1t+….+Υn−1tn−1(n−1)!,(12)

is a known vector function of time. Substituting the relation (11) into the governing Eq. (1), one gets the equation for the new variable P(t)

A^Dt(α)[P(t)]+∑k=1KA^k(t)Dt(αk)[P(t)]=F1(t)+B^(t)P(t),(13)

where

F1(t)=F(t)+B^(t)Θn−1(t)−∑k=1KA^k(t)Dt(αk)[Θn−1(t)].(14)

In should be noted that Dt(α)[Θn−1(t)]=0 because n−1<α≤n. The new the system is subjected to the zero ICs:

P(0)=0, ∂tP(0)=0,…, ∂t(n−1)P(0)=0.(15)

Let us rewrite the system in the form:

A^Dt(α)[P(t)]=B^(t)P(t)+F1(t)−∑i=1IA^i(t)Dt(αi)[P(t)].(16)

Let φm(t) be the system of basis functions on [0,T] which are chosen in such a way that the right-hand side of Eq. (16) can be represented in the form of the series

B^(t)P(t)+F1(t)−∑k=1KA^k(t)Dt(αk)[P(t)]=A^∑m=1∞qmφm(t),(17)

where qm = [qm,1,….,qm,N]T are N−vectors to be determined.

Throughout the paper, we use the generalized power functions or the Müntz polynomials basis (MPB) [76,77]. A fractional derivative of a Müntz polynomial is again a Müntz polynomial. This is a crucial feature of this base for using it in the collocation methods for Fractional Differential Equations (FDEs). So, we take

φm(t)=tδm,δm=σ(m−1),0<σ≤1,m=1,2,3,⋯(18)

as the basis functions and the solution is sought in the class of functions which can be approximated by the MPB and for which there exist fractional derivatives of the original Eq. (1).

Here σ is the parameter of the MPB. The BSM which uses the Müntz polynomials has been developed for solving single FODE in [78–80]. The results show that the Müntz polynomials provide quite an accurate approximation of the solution in the range of 0.10≤σ≤0.3. Throughout the paper, we use σ from this interval.

Under the condition (17) the system (16) can be written in the form:

A^Dt(α)[P(t)]=A^∑m=1∞qmφm(t).(19)

Suppose that the matrix A^ is invertible. So, we obtain the reduced matrix equation

Dt(α)[P(t)]=∑m=1∞qmφm(t).(20)

As it follows from Eq. (3) the analytical expression

Φm(t)=Γ(δm+1)Γ(δm+α+1)tδm+α,(21)

satisfies the FODE

Dt(α)[Φm(t)]=tδm=φm(t).(22)

Because n−1<α≤n the function Φm(t) satisfies zero ICs

Φm(0)=0,∂tΦm(0)=0,⋯,∂t(n−1)Φm(0)=0.(23)

Therefore, the linear combination

P∞(t)=∑m=1∞qmΦm(t),(24)

is the semi-analytical solution of Eq. (20) for any qm. It satisfies zero ICs Eq. (15). Let us emphasize: in general case, Eqs. (13) and (20) are different ones. However, if the relation (17) is fulfilled, they are identical. In this case P∞(t) is also the solution of (13) for any sequence qm. So, to get the vectors qm we substitute P∞(t) into the relation (17) and get the infinite system:

B^(t)∑m=1∞qmΦm(t)−A^1(t)∑m=1∞qmΦm(α1)(t)−⋯−A^K(t)∑m=1∞qmΦm(αK)(t)+F1(t)=A^∑m=1∞qmφm(t),(25)

or

∑m=1∞[A^φm(t)+A^1(t)Φm(α1)(t)+⋯+A^K(t)Φm(αK)(t)−B^(t)Φm(t)]qm=F1(t).(26)

If the Eq. (26) is fulfilled at any time moment t∈[0,T], then the vector function P∞(t) given in (24) is the exact solution of the problem (13), (15) if it exists. Then the sum (11) is the exact solution of the original problem (1), (10).

In practical calculations, we consider the truncated series

PM(t)=∑m=1MqmΦm(t),(27)

as an approximate solution of the problem (13), (15). It satisfies the truncated analog of the system (20):

Dt(α)[PM(t)]=∑m=1Mqmφm(t).(28)

and the unknown vectors qm are obtained by applying the collocation procedure to the truncated analog of (26):

∑m=1M[A^φm(tj)+A^1(tj)Φm(α1)(tj)+⋯+A^K(tj)Φm(αK)(tj)−B^(tj)Φm(tj)]qm=F1(tj)=deffj,j=1,⋯,Nc,(29)

where Nc is the number of the collocation nodes tj∈[0,T]. We use the Gauss-Chebyshev collocation points:

tj=0.5T[1+cos⁡(π(2j−1)/2Nc)]∈[0,T],j=1,2,⋯Nc.(30)

The collocation system (29) can be written in the compact form:

C^Q=ℱ,(31)

where

Q=[q1,q2,…,qM]T=[q1,1,…,q1,N,q2,1,…,q2,N,…,qM,1,…,qM,N]T,(32)

ℱ=[f1,f2,…,fN𝒞]T=[f1,1,…,f1,N,f2,1,…,f2,N,…,fN𝒞,1,…,fN𝒞,N]T.(33)

The collocation matrix C^ contains NcM blocks

C^=[C^1,1C^1,2…C^1,m…C^1,MC^2,1C^2,2…C^2,m…C^2,M………………C^Nc,1C^Nc,2…C^Nc,m…C^Nc,M](34)

Here M is the number of the Müntz polynomials φm(t) and Nc is the number of the collocation points in the time interval [0,T].

The matrices A^, A^i(tj), B^(tj) are square matrices of the size N×N and so, their combinations with the scalar coefficients φm(tj), Φm(αi)(tj), Φm(tj),

Cj,m=A^φm(tj)+A^1(tj)Φm(α1)(tj)+….+A^K(tj)Φm(αK)(tj)−B^(tj)Φm(tj),(35)

are also the N×N square matrices. Let us note that to obtain a stable solution of the collocation system the number of the collocation points Nc is taken twice as large as the number of the Müntz polynomials M: Nc=2M. Finally, we get the over-determined linear system (31) with the collocation matrix C^ which contains 2M2 square N×N blocks Cj,m j=1,….,2M, m=1,….,M. The matrix C^ has 2NM rows and NM columns. The collocation system is solved by the standard least squares procedure.

So, the algorithm of the solution of the system (31) is as follows:

Step 1. Choose the parameter of the Müntz polynomials basis,σ.

Step 2. Choose the number of the Müntz polynomials M in the approximate solution.

Step 3. Define the functions φm(t) and Φm(t), m=1,…,M (see (21)).

Step 4. Calculate the collocation matrix C^ using (34) and (35).

Step 5. Calculate the vector of the right hand side ℱ using (29), (33).

Step 6. Solve the collocation system (31) for Q and find the vectors qm, m=1,…,M which form Q (see (32)).

Step 7. Getting the functions Φm(t) m=1,…,M and the vectors qm, m=1,…,M obtain the approximate solution PM(t) (see (27)).

Step 8. Obtain the approximate solution of the original problem (1), (10) as the sum ΥM(t)=PM(t)+Θn−1(t) (see (11)).

3  Numerical Scheme for TFPDEs

Let us consider the TFPDE of the form:

Dt(α)[v(x,t)]+∑k=1Kak(t)Dt(αk)[v(x,t)]=b(t)L(x)[v(x,t)]+f(x,t),x∈[0,1] , t∈[0,T],(36)

subjected to the Dirichlet boundary conditions (BCs)

v(0,t)=g0(t),v(1,t)=g1(t),(37)

where L(x) is a spatial differential operator given in (5). The ICs are determined as follows:

v(x,0)=v0(x),∂tv(x,0)=v1(x),…,∂t(n−1)v(x,0)=vn−1(x).(38)

Let us define the new function

u(x,t)=v(x,t)−Ψn−1(x,t),(39)

where

Ψn−1(x,t)=v0(x)+v1(x)t+…+tn−1(n−1)!vn−1(x).

This function satisfies the equation

Dt(α)[u]+∑k=1Kak(t)Dt(αk)[u]=f1+b(t)L(x)[u],(40)

under the BCs

u(0,t)=g0(t)−Ψn−1(0,t)≡h0(t),(41)

u(1,t)=g1(t)−Ψn−1(1,t)≡h1(t).(42)

and zero ICs

u(x,0)=0,∂tu(x,0)=0,…,∂t(n−1)u(x,0)=0.(43)

Here

f1(x,t)=f(x,t)+L(x)[Ψn−1(x,t)]−∑k=1Kak(t)Dt(αk)[Ψn−1(x,t)].(44)

Note: The last terms in (44) can be expressed in the analytical form. Indeed,

∑i=1Iai(t)Dt(αi)[Ψn−1]=v1(x)∑αk≤1Kak(t)Dt(αk)[t]+12v2(x)∑αk≤2Kak(t)Dt(αk)[t2],…,+1(n−1)vn−1(x)∑αk≤n−1Kak(t)Dt(αk)[tn−1],(45)

where the derivative Dt(αi)[tj] can be written using (3). The previous term in (44) can be written as follows:

L(x)[Ψn−1(x,t)]=L(x)[v0(x)]+L(x)[v1(x)]t+⋯+tn−1(n−1)!L(x)[vn−1(x)].(46)

So, (44) is the analytical expression.

Let us define the function ug(x,t)

ug(x,t)=h0(t)+x(h1(t)−h0(t)),(47)

which satisfies the BCs (41), (42) and introduce the new variable w(x,t):

w(x,t)=u(x,t)−ug(x,t).(48)

The function w(x,t) is the solution of the TFPDE

Dt(α)[w(x,t)]+∑k=1Kak(t)Dt(αk)[w(x,t)]=b(t)L(x)[w(x,t)]+f2(x,t),(49)

where

f2(x,t)=f1(x,t)+L(x)[ug(x,t)]−Dt(α)[ug(x,t)]−∑k=1Kak(t)Dt(αk)[ug(x,t)].(50)

It is easily to prove that the function w(x,t) satisfies zero ICs and BCs:

w(x,0)=0,∂tw(x,0)=0,…,∂t(n−1)w(x,0)=0,(51)

w(0,t)=w(1,t)=0.(52)

Let us choose a set of linearly independent functions ψi(x), i=1,…,N defined in  [0,1]. For this goal we mainly use the Multiquadric radial basis function (MQ-RBF) throughout the paper. The centers ζj of the RBF are distributed in the solution region [0,1]:

ψj(x)=(x−ζj)2+c2,(53)

where c is the shape parameter. We place the centers of the RBFs at the Gauss-Chebyshev points

ζj=0.5[1+cos⁡(π(2j−1)/2N)]∈[0,1],j=1,2,…N.(54)

Based on the numerical experiments carried out we fix the shape parameter c=0.5 in all the calculations.

We define the modified basis functions

ϕj(x)=ψj(x)+bj,0+bj,1x,(55)

where the coefficients bj,0,bj,1 are chosen to satisfy BCs:

ϕj(0)=ϕj(1)=0.(56)

As a result we get the linear system

bj,0=−ψj(0),(57)

bj,0+bj,1=−ψj(1),(58)

for each pair of the coefficients bj,0,bj,1. The system can be solved easily in the analytical form. So, the modified basis functions ϕj(x) and their linear combinations satisfy zero boundary condition (56).

We seek the solution of the Eq. (49), in the form of the linear series over the modified basis functions ϕj(x)

wN(x,t)=∑j=1Nϕj(x)Υj(t).(59)

Substituting Eq. (59) into Eq. (49) we get

∑j=1Nϕj(x){Dt(α)[Υj(t)]+a1(t)Dt(α1)[Υj(t)]+…+aK(t)Dt(αK)[Υj(t)]}==b(t)∑j=1NL(x)[ϕj(x)]Υj(t)+f2(x,t).(60)

Let xi∈[0,1], i=1,…,N be the collocation points distributed inside the solution domain. Applying the collocation procedure at these points, we get the system of the FODEs:

∑j=1Nϕj(xi){Dt(α)[Υj(t)]+a1(t)Dt(α1)[Υj(t)]+…+aK(t)Dt(αK)[Υj(t)]}==b(t)∑j=1NL(xi)[ϕj(xi)]Υj(t)+f2(xi,t),(61)

We take the centers of the RBFs as the collocation points: xi=ζi. Let us rewrite the system of the FODEs in the vector form:

A^Dt(α)[Υ(t)]+A^1(t)Dt(α1)[Υ(t)]+…++A^K(t)Dt(αK)[Υ(t)]=B^(t)Υ(t)+F(t),t∈[0.T],(62)

where A^, A^i(t),B^(t) are N×N matrices with the components

A^ =[ϕj(xi)]i,j=1N,A^k(t)=ak(t)A^,B^(t)=b(t)[L(xi)[ϕj(xi)]]i,j=1N,(63)

and Υ(t)=[Υ1(t),Υ2(t),….,ΥN(t)]T, F(t)=[f2(x1,t),f2(x2,t),…,f2(xN,t)]T are N-vectors. Note that the derivatives in the term L(xj)[ϕl(xj)] can be obtained in the analytical form

∂xϕj(x)=x−ζjψj(x)+cj,1,∂xxϕj(x)=1ψj(x)−(x−ζj)2(ψj(x))3.(64)

So, the system (62) takes the same form as the linear system of FODEs (1) and can be solved by the algorithm described in Section 2. It should be noted that taking into account (51), the vector Υ(t) satisfies zero ICs:

Υ(0)=0, ∂tΥ(0)=0,…,∂t(n−1)Υ(0)=0.(65)

It means that Θn−1(t)=0 (see (11)) and Υ(t)=P(t). Using M the Müntz polynomials, we get the approximate solution in the form:

wN,M(x,t)=∑j=1Nϕj(x)PM,j(t),(66)

where PM,j(t) is jth component of the vector PM(t) given in (27). Then, the approximate solution of the original problem (36)–(38) can be written in the form

vN,M(x,t)=uN,M(x,t)+Ψn−1(x,t)=wN,M(x,t)+ug(x,t)+Ψn−1(x,t)=∑j=1Nϕj(x)PM,j(t)+ug(x,t)+Ψn−1(x,t).(67)

3.1 Nonlinear Problem

Let us consider TFPDE (4) with a nonlinear term

Dt(α)[v(x,t)]+∑k=1Kak(t)Dt(αk)[v(x,t)]=b(t)L(x)[v(x,t)]+G(v)+f(x,t).(68)

Let v0(x,t) be a function considered as the initial approximate solution. Linearization of the function G(v) in the vicinity of v0(x,t)

G(v)≃G(v0)+∂G∂v(v0)[v−v0]=∂G∂v(v0)v++G(v)−∂G∂v(v0)v0

transforms (68) into a sequence of linear TFPDEs each of those can be solved by the technique described above. As a result, we get the iteration procedure. The iterations are stopped with the control of the error maxx,t|v(x,t)−v(x,t)| or after achieving the prescribed number of iterations.

4  Numerical Examples

In this section several numerical examples are provided to show the accuracy of the proposed scheme. To demonstrate the performance of this technique we consider the different types of errors for systems of FODEs and TFPDEs.

emax(Υk,Nt)=max1≤i≤Nt|Υk,ex(ti)−Υk,M(ti)|,(69)

eRMS(Υk,Nt)=1Nt∑i=1Nt(Υk,ex(ti)−Υk,M(ti))2.(70)

The errors (69), (70) are used in solving systems of the FODEs to estimate the approximate solution of each component of the vector Υ(t)=[Υ1(t), Υ2(t),…,ΥN(t)]T (see Eq. (1)).

Emax(t,Nt)=max1≤i≤Nt|vN,M(xi,t)−vex(xi,t)|,(71)

ERMS(t,Nt)=1Nt∑i=1Nt[(vex(xi,t)−vN,M(xi,t))2+(∂vex∂x(xi,t)−∂vN,M∂x(xi,t))2].(72)

The errors (71), (72) are used in solving TFPDE (36) with the BCs (37) and the ICs (38). The error ERMS(t,Nt) also includes the error in the approximation of the first derivative of the solution. The subscript ex of v indicates the comparison with the analytical solution or with the data taken from the literature. The subscript N, M of v defines the numerical solution. For (1+1) dimensional problems Nt =4N. The numerical experiments show that this relation guarantees the accuracy in the calculation of the errors. To carry out convergence research, we define the convergence order (CO).

CO=log⁡(Emax(N1)/Emax(N2))log⁡(N1/N2).(73)

4.1 Numerical Experiments for Systems of FODEs

Example 4.1 Let us consider the system given in [35]

Dt(α)[Υ1(t)Υ2(t)]=B^[Υ1(t)Υ2(t)],B^=[2,−14,−3],0<α≤1,(74)

with the general exact solution

[Υ1(t)Υ2(t)]=c1[11]Eα(tα)+c2[14]Eα(−2tα),(75)

where Eα(t) is the one parameter Mittag-Leffler function

Eα(t)=∑j=0∞tjΓ(αj+1),

and c1 and c2 are free parameters which define the ICs. We solve this problem by the suggested scheme for α=0.25, 0.5 and 0.75. The obtained results are reported in Table 1. We use the parameter of the MPB σ=0.3 and Nt=50000 test points distributed randomly inside [0,1] to calculate the errors (69), (70).

Example 4.2 Let us consider the system described in [35]

Dt(α)[Υ1(t)Υ2(t)Υ3(t)]=B^[Υ1(t)Υ2(t)Υ3(t)],B^=[1,1,12,1,−10,−1,1],0<α≤1,(76)

with the general exact solution

[Υ1(t)Υ2(t)Υ3(t)]=c1[−342]Eα(−tα)+c2[01−1]Eα(2tα)(77)

+c3{[01−1]tα∂tEα(2tα)+[101]Eα(2tα)},(78)

where ∂tEα(t) is the first derivative of the one parameter of the Mittag-Leffler function

∂tEα(t)=∑j=0∞(j+1)tjΓ(α(j+1)+1),

and c1, c2, c3 are free parameters which define the ICs. We solve this problem by the suggested scheme for α=0.5 and c1=c2=c3=1. The obtained results are reported in Table 2.

Example 4.3 Consider the following initial value problem for the inhomogeneous Bagley–Torvik equation [81]

d2V(t)dt2+Dt(3/2)[V(t)]+V(t)=t+1,V(0)=1,dV(0)dt=1,(79)

with the exact solution V(t)=t+1. As it is shown in [82] the original Bagley–Torvik equation can be rewritten as the system of FODEs of order 1/2

Dt(1/2)[Υ1(t)Υ2(t)Υ3(t)Υ4(t)]=B^[Υ1(t)Υ2(t)Υ3(t)Υ4(t)]+F(t),B^=[0,1,0,00,0,1,00,0,0,1−1,0,0,−1],F(t)=[000t+1],(80)

with the initial condition

Υ1(0)=1,Υ2(0)=0 , Υ3(0)=1,Υ4(0)=0,(81)

where the solution of the original Bagley–Torvik equation V(t)=Υ1(t). It can be easily proved that the exact solution of the system (80) is

Υex(t)=[1+t,2t1/2/π,1,0]T.(82)

According to the method described above, the approximate solution can be written in the form:

ΥM(t)=Υex(0)+∑m=1MqmΦm(t),(83)

where

Φm(t)=Γ(δm+1)Γ(δm+3/2)tδm+1/2,δm=σ(m−1).(84)

As it follows from (82), (83)

ΥM(t)−Υex(0)=[t,2t1/2/π,0,0]T=∑m=1MqmΦm(t).

Thus, the approximate solution (83) contains the exact solution Υex(t), if the sequence Φm(t), m=1,…,M contains the functions t1/2 and t. As it follows from (84), this condition is fulfilled when the sequence δm+1/2=σ(m−1)+1/2, m=1,….,M contains the values 1/2 and 1. For example, if σ=1/2, then δm+1/2=1/2,1,3/2,…. and the expression (83) contains the exact solution Υex(t) beginning from M=2. For σ=1/4, we get the sequence δm+1/2=1/2, 3/4, 1,… and for σ=0.1, we get the sequence δm+1/2=1/2,1/2+0.1,1/2+0.2,1/2+0.3,1/2+0.4,1/2+0.5=1,…. Thus, for σ=1/4 the exact solution Υex(t) is included in (83) beginning from M=3 and for σ=0.1 beginning from M=6.

The data placed in Table 3 demonstrate that if the approximate solution (83) contains the exact solution Υex(t), then the method calculates it up to the machine precision.

The data placed in the last rows of Table 3 correspond to the general case when the information of the solution is absent. We take σ=π/20 and there are no sequences Φm(t), m=1,…,M which contain the functions t1/2 and t. As a result, the error decreases slowly and gradually. The method also provides a high accuracy but with larger values of M.

The same problem has been studied in [81] on the time interval t∈[0,5] using fractional linear multistep methods based on the Adams-type predictor-corrector technique. As demonstrated in the paper, the errors depend on time step size Δt. For time step Δt=0.5 emax(Υ1)=−0.15131473519232, and for Δt=0.0625 emax(Υ1)=−0.00562770408881. These data also are placed in Table C.3 of [2]. Note that the component Υ1(t) corresponds to the solution of the original Bagley–Torvik equation (79).

Table 4 shows the results of the calculation by the proposed method on the time interval t∈[0,5] with the parameter of the MPB σ=π/20. It is obvious that the method presented provides a much more accurate solution. With a special choice of parameter σ the accuracy is even higher. For example, if σ=0.25, M=3, then emax(Υ1) = 4.44E-15 and eRMS(Υ1) = 1.91E-15.

Example 4.4 Let us consider the multi-term system with time-dependent matrices

A^Dt(3)[Υ1(t)Υ2(t)]+A^1(t)Dt(π/2)[Υ1(t)Υ2(t)]+A^2(t)Dt(0.75)[Υ1(t)Υ2(t)]+A^3(t)Dt(0.25)[Υ1(t)Υ2(t)]=B^(t)[Υ1(t)Υ2(t)]+F(t).(85)

The ICs are

Υ1(0)=1,Υ2(0)=0.(86)

Her

A^=[1,0.10.1,1],A^1(t)=11+tA^,A^2(t)=t1+t2A^,A^3(t)=t21+t3A^,B^(t)=et[1,−0.1−0.1,1].(87)

The vector F(t) corresponds to the exact solution

Υex(t)=[exp⁡(t)sin⁡(t)].(88)

The data placed in Table 5 demonstrate the behavior of the error of the approximate solution with the growth of M for σ=0.1, 0.2 and 0.3. The method converges quite fast for all σ.

Similar to (83), (84), the approximate solution can be written in the form:

ΥM(t)=Υex(0)+∂tΥex(0)t+∑m=1MqmΦm(t).(89)

where

Φm(t)=Γ(δm+1)Γ(δm+3+1)tδm+3.(90)

In this case we use the additional information that the components of the solution Υ1(t), Υ2(t) are analytical functions and can be approximated by the polynomials 1, t, t2,…, tM+1. If we set δm=m+1−3, then Φm(t)∼tm+1 and so, the exact solution Υex(t) belongs to linear span of the polynomials 1, t, t2,…, tM+1.

Table 6 demonstrates a dramatic decrease in the errors with the growth of M for this special choice of δm when an additional information on the solution is added.

4.2 Numerical Experiments for TFPDEs

Example 4.5 Let us consider the multi-term TFPDE [83]

Dt(α)[v]+Dt(0.2)[v]=∂xxv+f(x,t),0≤x,t≤1.(91)

Here the source term f(x,t), IC and BCs conform to the exact solution v(x,t)=(1+t2)(x2−x).

Fig. 1 and Table 7 show the behavior of the errors of the approximate solution as the functions of N (see (59)) with the fixed M. For M=10 the errors decrease with the growth of N in the whole range 2≤N≤26. This means that the error of the spatial approximation is the dominant error. While for M=5 the accuracy does not improve for N>10. Therefore, in this case, the dominant error is the error in the solution of the system of FODEs in time. Fig. 1 shows that all the curves Emax(N), ERMS(N) originally lay on the same curve. They move away from this curve depending on the value of M when the error due to approximation in time becomes dominant.

Figure 1: The maximal absolute Emax (left) and ERMS (right) errors as functions of the number of RBFs used in the approximate solution. α=0.95, σ=0.3

Fig. 2 and Table 8 show the behavior of the errors as functions of M with the fixed N. It is evident that the proposed scheme converges fast with the increase of M. For the larger M more accurate results can be obtained. The same problem was studied by Jin et al. in [83] using the Galerkin FE method and FD discretization of the time-fractional derivatives. Using h=2−10 mesh size and the time step size τ=1/160, they obtained the data placed in the last row of Table 7. The comparison shows that the method presented provides a much more accurate solution.

Figure 2: The maximal absolute Emax (left) and ERMS (right) errors vs. M with fixed N. α=0.95, σ=0.3

Example 4.6 Let us consider the multi-term TFPDE

Dt(π/4)[v]+t1+tDt(0.5)[v]+t21+t2Dt(0.25)[v]=sinh⁡(t)Dt(0.2)[L(x)[v]]+cosh⁡(t)Dt(0.1)[L(x)[v]]+(1+t2)L(x)[v]+f(x,t),0≤x,t≤1,

with the spatial operator L(x)[v(x,t)]=∂x(cos⁡(x)∂xv(x,t)). The Dirichlet BCs and IC conform to the exact solution v(x,t)=sin⁡(x+t).

Table 9 shows the errors, convergence order, and CPU time as the functions of N with the fixed M. The data also are illustrated by the graphics in Fig. 3. With increasing of N the proposed method converges fast, and we can obtain the errors around 10−9 with M=10 and N=20. It should also be noted that with a small number of N=4, the computed errors are around 10−3, which should be sufficient for engineering applications.

Figure 3: The maximal absolute Emax (left) and ERMS (right) errors vs. N for different M

Table 10 and Fig. 4 show the errors vs. M with the fixed N to verify the performance of the proposed scheme. The order of convergence is larger than 3.

Figure 4: The maximal absolute Emax (left) and ERMS (right) errors vs. M for different fixed N

Example 4.7 In the following examples we consider three cases for α>1 to verify the performance of the proposed scheme.

Case 1: Consider the following equation:

Dt(π/2)[v]+11+tDt(e/2)[v]+t1+t2Dt(π/4)[v]+t21+t3Dt(e/4)[v]=∂xx2v+f(x,t),0≤x≤1,0≤t≤T.(92)

Here we have 1<α<2 and the equation needs two ICs

v(x,0)=cos⁡x,∂v(x,0)∂t=−sin⁡x.

The boundary conditions are

v(0,t)=cos⁡(t),v(1,t)=cos⁡(t+1).

The exact solution of the problem is v(x,t)=cos⁡(x+t).

Case 2: Consider the following equation:

Dt(5)[v]+cos⁡t1+tDt(e/2)[v]+sin⁡t1+t2Dt(π/4)[v]+t21+t3Dt(e/4)[v]=etL(x)[v]+f(x,t),0≤x≤1,0≤t≤T,

with the spatial operator L(x)[u(x,t)]=∂x(cosh⁡(x)∂xu(x,t)). Here we have 2<α<3 and the equation needs three ICs

v(x,0)=sin⁡x,∂v(x,0)∂t=cos⁡x,∂2v(x,0)∂t2=−sin⁡x.

The boundary conditions are

v(0,t)=sin⁡(t),v(1,t)=sin⁡(t+1).

The exact solution of the problem is v(x,t)=sin⁡(x+t).

Case 3: Consider the following equation:

Dt(π)[v]+cosh⁡t1+tDt(5)[v]+sinh⁡t1+t2Dt(e/2)[v]+t21+t3Dt(π/4)[v]=etL(x)[v]+f(x,t),0≤x≤1,0≤t≤T,(93)

with the spatial operator L(x)[v(x,t)]=∂x(sinh⁡(x)∂xv(x,t)). Here we have 3<α<4 and the equation needs four ICs

v(x,0)=sin⁡x,∂v(x,0)∂t=cos⁡x,∂2v(x,0)∂t2=−sin⁡x,∂3v(x,0)∂t3=−cos⁡x.(94)

The boundary conditions are

v(0,t)=sin⁡(t),v(1,t)=sin⁡(t+1).(95)

The exact solution of the problem is v(x,t)=sin⁡(x+t).

Tables 11, 12 show the errors with increasing of N with the fixed M. It is evident that the proposed scheme provides very accurate results. Furthermore, for a small number of N, we can also get moderately accurate results with errors around 10−6. For larger N the proposed scheme converges faster. The accuracy of the proposed method is also demonstrated in Fig. 5. Finally, Table 13 shows the error when the Gaussian RBF exp⁡(−(xc)2) is used in spatial aproximation. For a small number of N the Gaussian RBF provides a more accurate solution than the MQ-RBF. For large N both RBFs provide the results with errors of the same level of accuracy.

Figure 5: The domain absolute errors

Example 4.8 Consider the nonlinear time-fractional the Huxley-Burgers’ equation of the following form:

Dt(α)[v]=∂xxv+v(1−v)(v−γ)+f(x,t),0≤x,t≤1.

The Dirichlet BCs and IC conform to the exact solution v(x,t)=x3(t+1).

Tables 14, 15 and Fig. 6 show the behaviour of the errors with the growth of N and the fixed M. The data are obtained after 3 iterations of the quazilinearization procedure. The same problem was considered by Hadhoud et al. in [84] using a numerical technique based on the cubic B-spline collocation method and the mean value theorem for integrals. The maximal absolute errors obtained there for the mesh size Δx=0.01 and the time step Δt=0.01 are shown in the last rows of the tables. The last columns of the tables contain the data corresponding to α=1, i.e., to the solution of the equation of the integer order. These data demonstrate that the proposed method can be used for solving equations of fractional order as well as equations of integer order without any modification of the algorithm.

Figure 6: The maximal absolute Emax (left) and ERMS (right) errors as functions of the number of the RBFs N used in the approximate solution. The data correspond to α=0.5, γ=0.1

Example 4.9 Consider the TFPDE

Dt(α)[v]=∂xxv+f(x,t),0≤x,t≤1.

The source function f(x,t), Dirichlet BCs v(0,t)=v(1,t)=0 and IC v(x,0)=0 conform the exact solution v(x,t)=x4(1−x)tα with the strong singularity at t=0.

Table 16 shows the behavior of the errors with the growth of N with the fixed M=12 and with the parameter of the MPB σ=0.3. The same problem was considered by Ferrás et al. in [85] using a numerical technique based on the combination of the method of lines with the hybrid collocation method. The maximal absolute errors obtained there are shown in the last row of the table.

5  Conclusion

This paper presents a new meshless technique for solving multi-term linear systems of fractional equations. These systems have been used in modeling various phenomena in different branches of engineering and science. Using substitution (11), we transform the original system into the one for the vector variable P(t) which satisfies zero ICs. Then P(t) is sought in the form of the finite series over the Müntz polynomials basis. Applying the collocation procedure in the domain, we get the linear algebraic system solved by the standard numerical procedure. Then, on the base of this technique, the method of solving the TFPDE has been developed. The collocation at the centers of the RBF transforms the TFPDE into a system of FODEs similar to the one considered in Section 2.

In the authors’ opinion, the main results achieved in the paper are: (1) The effective method for solving systems of the FODEs with time-dependent coefficients has been developed and tested. (2) On the base of this technique the method of solving TFPDEs of the high fractional order has been proposed. The method has been tested on the problems with the highest derivative of the orders: 1<α<2, 2<α<3 and 3<α<4. (3) The technique has been extended to nonlinear the Huxley-Burgers’ TFPDE. It should be stressed that the proposed method can be used for solving both equations of fractional and integer order without any modification of the algorithm.

Some remarks: (1) In this paper the MQ-RBF is mainly used for spatial approximation. However, the last example demonstrates that the Gaussian RBF is suitable for this purpose. The other global RBFs, compactly supported RBFs and B-splines also can be used for spatial approximation in the framework of the proposed technique. (2) Only the Dirichlet BCs are considered in this study. However, the proposed technique can be extended to the problems with the boundary conditions of the general type by some modification of the Eqs. (56)–(58). (3) Only (1 + 1) dimensional problems have been considered. However, using multidimensional RBSs, this approach can be extended to the (2 + 1) and (3 + 1) dimensional problems.

It should be remarked that the limitation of the presented technique is caused by the fast growth of the size of the collocation matrix C^ of the linear system (31) with the growth of N and M. Because C^ is the dense matrix the problem of ill-conditioning also arises.

To overcome this problem we presuppose the use of a localized scheme of the spatial approximation based on the compactly supported radial basis functions (CSRBF) in the future to avoid dense and ill-conditioning matrices.

To overcome the problems of calculations on the large time interval [0,Tmax], we think of using the approach developed in [86].

Let [0,Tmax] represent a sum of the subintervals:

[0,Tmax]=[0,T]∪[T,2T]∪….∪[(L−1)T,LT].

Solving the Eq. (1) in the first subinterval [0,T], we use the ICs (10) of the original problem. As a result of the solution we get the vectors

Υ01=Υ(T),Υ11=∂tΥ(T)… . ,Υn−11=∂t(n−1)Υ(T),

which can be used as the initial data for solving the equation in the second subinterval [T,2T]. Applying the transform

t=τ+T,

We get the equation for the unknown vector Υ2(τ) defined on the interval τ∈[0,T]:

A^Dt(α)[Υ2(τ)]+∑k=1KA^k(τ+T)Dt(αk)[Υ2(τ)]=B^(τ+T)Υ2(τ)+F(τ+T),0≤τ≤T,

with the ICs

Υ2(0)=Υ01,∂tΥ2(0)=Υ11… . ,∂t(n−1)Υ2(0)=Υn−11.

Then, the vectors Υ2(T), ∂tΥ2(T),…, ∂t(n−1)Υ2(T) are used as the ICs for solving on the interval t∈[2T,3T], etc. So, we solve the original equation at the same time interval τ∈[0,T] with the time-dependent coefficients computed at the shifted time τ+lT. The calculations are continued till t=LT=Tmax. We can choose T small enough to reduce the value M and so the size of the collocation matrix. It is worth emphasizing that in the paper mentioned above the system of 3 FODEs has been solved on the time interval [0,Tmax]=[0,5000]. All these items mentioned above will be the subjects of the further study.

Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions whose efforts have helped to improve the quality of this paper.

Funding Statement: This research was funded by the National Key Research and Development Program of China (No. 2021YFB2600704), the National Natural Science Foundation of China (No. 52171272), and the Significant Science and Technology Project of the Ministry of Water Resources of China (No. SKS-2022112).

Availability of Data and Materials: Sergiy Reutskiy: Methodology and Writing–Original Draft; Yuhui Zhang: Writing–Review & Editing; Jun Lu: Writing–Review & Editing; Ciren Pubu: Validation.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Podlubny, I. (1999). Fractional differential equations. New York: Academic Press.

2. Diethelm, K. (2010). The analysis of fractional differential equations. In: Lecture notes in mathematics, vol. 2004. Berlin: Springer.

3. Biazar, J., Farrokhi, L., Islam, M. R. (2006). Modeling the pollution of a system of lakes. Applied Mathematics and Computation, 178(2), 423–430.

4. Khader, M. M., El Danaf, T. S., Hendy, A. S. (2013). A computational matrix method for solving systems of high order fractional differential equations. Applied Mathematical Modelling, 37(6), 4035–4050.

5. El-Dessoky Ahmed, M. M., Khan, M. A. (2020). Modeling and analysis of the polluted lakes system with various fractional approaches. Chaos, Solitons & Fractals, 134, 109720.

6. Qin, S., Liu, F., Turner, I., Vegh, V., Yu, Q. et al. (2017). Multi-term time-fractional Bloch equations and application in magnetic resonance imaging. Journal of Computational and Applied Mathematics, 319, 308–319.

7. Cardoso, L. C., Dos Santos, F. L. P., Camargo, R. F. (2018). Analysis of fractional-order models for hepatitis B. Computational and Applied Mathematics, 37, 4570–4586.

8. Ding, Y., Ye, H. (2009). A fractional-order differential equation model of HIV infection of CD4+ T-cells. Mathematical and Computer Modelling, 50(3–4), 386–392.

9. Magin, R., Feng, X., Baleanu, D. (2009). Solving the fractional order Bloch equation. Concepts in Magnetic Resonance Part A: An Educational Journal, 34(1), 16–23.

10. Yu, Q., Liu, F., Turner, I., Burrage, K. (2014). Numerical simulation of the fractional Bloch equations. Journal of Computational and Applied Mathematics, 255, 635–651.

11. Xu, C., Yu, Y., Ren, G., Sun, Y., Si, X. (2023). Stability analysis and optimal control of a fractional-order generalized SEIR model for the COVID-19 pandemic. Applied Mathematics and Computation, 457, 128210. https://doi.org/10.1016/j.amc.2023.128210

12. Batiha, I. M., Momani, S., Ouannas, A. (2022). Fractional-order COVID-19 pandemic outbreak: Modeling and stability analysis. International Journal of Biomathematics, 15(1), 2150090.

13. Ndaïrou, F., Area, I., Nieto, J. J., Silva, C. J., Torres, D. F. M. (2021). Fractional model of COVID-19 applied to Galicia, Spain and Portugal. Chaos, Solitons & Fractals, 144, 110652. https://doi.org/10.1016/j.chaos.2021.110652

14. Ogunrinde, R. B., Nwajeri, U. K., Fadugba, S. E., Ogunrinde, R. R., Oshinubi, K. I. (2021). Dynamic model of COVID-19 and citizens reaction using fractional derivative. Alexandria Engineering Journal, 60, 2001–2012. https://doi.org/10.1016/j.aej.2020.09.016

15. Zeng, F. H., Li, C. P., Liu, F. W., Turner, I. (2015). Numerical algorithms for time-fractional subdiffusion equation with second-order accuracy. SIAM Journal on Scientific Computing, 37, A55–A78.

16. Das, S., Vishal, K., Gupta, P. K., Yildirim, A. (2011). An approximate analytical solution of time-fractional telegraph equation. Applied Mathematics and Computation, 217, 7405–7411.

17. Hosseini, V. R., Chen, W., Avazzadeh, Z. (2014). Numerical solution of fractional telegraph equation by using radial basis functions. Engineering Analysis with Boundary Elements, 38, 31–39.

18. Mollahasani, N., Mohseni Moghadam, M., Afrooz, K. (2016). A new treatment based on hybrid functions to the solution of telegraph equations of fractional order. Applied Mathematical Modelling, 5–6, 2804–2814.

19. Liu, F., Meerschaert, M. M., McGough, R. J., Zhuang, P., Liu, Q. (2013). Numerical methods for solving the multi-term time-fractional wave-diffusion equations. Fractional Calculus and Applied Analysis, 16, 9–25.

20. Zheng, M., Liu, F., Anh, V., Turner, I. (2016). A high-order spectral method for the multi-term time-fractional diffusion equations. Applied Mathematical Modelling, 40, 4970–4985.

21. Dehghan, M., Safarpoor, M., Abbaszadeh, M. (2015). Two high-order numerical algorithms for solving the multi-term time fractional diffusion-wave equations. Journal of Computational and Applied Mathematics, 290, 174–195.

22. Liu, L., Yang, C., Burrage, K. (2009). Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term. Journal of Computational and Applied Mathematics, 231, 160–176.

23. Liu, Q., Liu, F., Turner, I., Anh, V. (2011). Finite element approximation for a modified anomalous subdiffusion equation. Applied Mathematical Modelling, 35, 4103–4116.

24. Chen, W., Sun, H., Zhang, X., Korošakb, D. (2010). Anomalous diffusion modeling by fractal and fractional derivatives. Journal of Computational and Applied Mathematics, 59, 1754–1758.

25. Mohebbi, M., Abbaszadeh, M., Dehghan, M. (2013). A high-order and unconditionally stable scheme for the modified anomalous fractional sub-diffusion equation with a nonlinear source term. Journal of Computational Physics, 240, 36–48.

26. Dehghan, M., Abbaszadeh, M., Mohebbi, A. (2016). Legendre spectral element method for solving time fractional modified anomalous sub-diffusion equation. Applied Mathematical Modelling, 40(5–6), 3635–3654.

27. Almeida, R., Bastos, N. R. O., Monteiro, M. T. T. (2016). Modeling some real phenomena by fractional differential equations. Mathematical Methods in the Applied Sciences, 239(16), 4846–4855.

28. Sun, H., Zhang, Y., Baleanu, D., Chen, W., Chen, Y. (2018). A new collection of real world applications of fractional calculus in science and engineering. Communications in Nonlinear Science and Numerical Simulation, 64, 213–231.

29. Hilfer, R. (2000). Applications of fractional calculus in physics. Singapore: World Scientific.

30. Rossikhin, Y. A., Shitikova, M. V. (2010). Application of fractional calculus for dynamic problems of solid mechanics: Novel trends and recent results. Applied Mechanics Reviews, 63(1), 010801.

31. Kulish, V. V., Lage, J. L. (2002). Application of fractional calculus to fluid mechanics. Journal of Fluids Engineering, 124(3), 803–806.

32. Yang, X. J. (2019). General fractional derivatives: Theory, methods and applications. Boca Raton, FL: CRC Press.

33. Mainardi, F. (2022). Fractional calculus and waves in linear viscoelasticity: An introduction to mathematical models. Singapore: World Scientific.

34. Kexue, L., Jigen, P. (2011). Laplace transform and fractional differential equations. Applied Mathematics Letters, 24(12), 2019–2023.

35. Odibat, Z. M. (2010). Analytic study on linear systems of fractional differential equations. Computers & Mathematics with Applications, 59(3), 1171–1183.

36. Mamchuev, M. O. (2017). Solutions of the main boundary value problems for the time-fractional telegraph equation by the Green function method. Fractional Calculus and Applied Analysis, 20, 190–211.

37. Bailey, D. H., Swarztrauber, P. N. (1991). The fractional Fourier transform and applications. SIAM Review, 33(3), 389–404.

38. Das, S. (2009). Analytical solution of a fractional diffusion equation by variational iteration method. Computers & Mathematics with Applications, 57(3), 483–487.

39. Momani, S., Odibat, Z. (2006). Analytical solution of a time-fractional Navier-Stokes equation by Adomian decomposition method. Applied Mathematics and Computation, 177(2), 488–494.

40. Wu, C., Rui, W. (2018). Method of separation variables combined with homogenous balanced principle for searching exact solutions of nonlinear time-fractional biological population model. Communications in Nonlinear Science and Numerical Simulation, 63, 88–100.

41. Li, C., Zeng, F. (2012). Finite difference methods for fractional differential equations. International Journal of Bifurcation and Chaos, 22(4), 1230014.

42. Moaddy, K., Momani, S., Hashim, I. (2011). The non-standard finite difference scheme for linear fractional PDEs in fluid mechanics. Computers & Mathematics with Applications, 61(4), 1209–1216.

43. Baleanu, D., Zibaei, S., Namjoo, M., Jajarmi, A. (2011). A nonstandard finite difference scheme for the modeling and nonidentical synchronization of a novel fractional chaotic system. Advances in Difference Equations, 2021(1), 308.

44. Hajipour, M., Jajarmi, A., Baleanu, D. (2018). An efficient nonstandard finite difference scheme for a class of fractional chaotic systems. Journal of Computational and Nonlinear Dynamics, 13, 021013.

45. Wang, H., Basu, T. S. (2012). A fast finite difference method for two-dimensional space-fractional diffusion equations. SIAM Journal on Scientific Computing, 34(5), A2444–A2458.

46. Wang, H., Du, N. (2013). A fast finite difference method for three-dimensional time-dependent space-fractional diffusion equations and its efficient implementation. Journal of Computational Physics, 253, 50–63.

47. Jia, J., Wang, H. (2018). A fast finite difference method for distributed-order space-fractional partial differential equations on convex domains. Computers & Mathematics with Applications, 75(6), 2031–2043.

48. Chen, S., Liu, F., Jiang, X., Turner, I., Burrage, K. (2016). Fast finite difference approximation for identifying parameters in a two-dimensional space-fractional nonlocal model with variable diffusivity coefficients. SIAM Journal on Numerical Analysis, 54(2), 606–624.

49. Zhao, M., Wang, H., Cheng, A. (2018). A fast finite difference method for three-dimensional time-dependent space-fractional diffusion equations with fractional derivative boundary conditions. Journal of Scientific Computing, 74, 1009–1033.

50. Lyu, P., Liang, Y., Wang, Z. (2020). A fast linearized finite difference method for the nonlinear multi-term time-fractional wave equation. Applied Numerical Mathematics, 151, 448–471.

51. Shen, J., Sun, Z., Du, R. (2018). Fast finite difference schemes for time-fractional diffusion equations with a weak singularity at initial time. East Asian Journal of Applied Mathematics, 8(4), 834–858.

52. Zhuang, P., Liu, F., Turner, I., Gu, Y. (2014). Finite volume and finite element methods for solving a one-dimensional space-fractional Boussinesq equation. Applied Mathematical Modelling, 38(15–16), 3860–3870.

53. Cao, W., Hao, Z., Zhang, Z. (2022). Optimal strong convergence of finite element methods for one-dimensional stochastic elliptic equations with fractional noise. Journal of Scientific Computing, 91(1), 1.

54. Zhao, Y., Bu, W., Huang, J., Liu, D. Y., Tang, Y. (2015). Finite element method for two-dimensional space-fractional advection-dispersion equations. Applied Mathematics and Computation, 257, 553–565.

55. Bu, W., Tang, Y., Wu, Y. Yang. J. (2015). Finite difference/finite element method for two-dimensional space and time fractional Bloch-Torrey equations. Journal of Computational Physics, 293, 264–279.

56. Feng, L., Liu, F., Turner, I., Yang, Q., Zhuang, P. (2018). Unstructured mesh finite difference/finite element method for the 2D time-space Riesz fractional diffusion equation on irregular convex domains. Applied Mathematical Modelling, 59, 441–463.

57. Jin, B., Lazarov, R., Pasciak, J., Zhou, Z. (2014). Error analysis of a finite element method for the space-fractional parabolic equation. SIAM Journal on Numerical Analysis, 52(5), 2272–2294.

58. Huang, C., Stynes, M. (2021). α-robust error analysis of a mixed finite element method for a time-fractional biharmonic equation. Numerical Algorithms, 87, 1749–1766.

59. Li, X., Yang, X., Zhang, Y. (2017). Error estimates of mixed finite element methods for time-fractional Navier-Stokes equations. Journal of Scientific Computing, 70, 500–515.

60. Mohebbi, A., Abbaszadeh, M., Dehghan, M. (2013). The use of a meshless technique based on collocation and radial basis functions for solving the time fractional nonlinear Schrodinger equation arising in quantum mechanics. Engineering Analysis with Boundary Elements, 37(2), 475–485.

61. Kumar, A., Bhardwaj, A., Kumar, B. V. R. (2019). A meshless local collocation method for time fractional diffusion wave equation. Computers & Mathematics with Applications, 78(6), 1851–1861.

62. Faghih, A., Mokhtary, P. (2021). A new fractional collocation method for a system of multi-order fractional differential equations with variable coefficients. Journal of Computational and Applied Mathematics, 383, 113139. https://doi.org/10.1016/j.cam.2020.113139

63. Hu, W., Fu, Z., Tang, Z., Gu, Y. (2022). A meshless collocation method for solving the inverse Cauchy problem associated with the variable-order fractional heat conduction model under functionally graded materials. Engineering Analysis with Boundary Elements, 140, 132–144.

64. Habibirad, A., Hesameddini, E., Azin, H., Heydari, M. H. (2023). The direct meshless local Petrov-Galerkin technique with its error estimate for distributed-order time fractional Cable equation. Engineering Analysis with Boundary Elements, 150, 342–352.

65. Abbaszadeh, M., Dehghan, M. (2020). Direct meshless local Petrov-Galerkin (DMLPG) method for time-fractional fourth-order reaction-diffusion problem on complex domains. Computers & Mathematics with Applications, 79(3), 876–888.

66. Li, X., Li, S. (2021). A fast element-free Galerkin method for the fractional diffusion-wave equation. Applied Mathematics Letters, 122, 107529.

67. Li, X. (2023). A stabilized element-free Galerkin method for the advection–diffusion–reaction problem. Applied Mathematical Letters, 146, 108831. https://doi.org/10.1016/j.aml.2023.108831

68. Shah, K., Alqudah, M. A., Jarad, F., Abdeljawad, T. (2020). Semi-analytical study of Pine Wilt Disease model with convex rate under Caputo-Febrizio fractional order derivative. Chaos, Solitons & Fractals, 135, 109754.

69. Khater, M. M., Bekir, A., Lu, D., Attia, R. A. (2021). Analytical and semi-analytical solutions for time-fractional Cahn-Allen equation. Mathematical Methods in the Applied Sciences, 44(3), 2682–2691.

70. Kazem, S., Dehghan, M. (2019). Semi-analytical solution for time-fractional diffusion equation based on finite difference method of lines (MOL). Engineering with Computers, 35, 229–241.

71. Arafa, A., Hagag, A. (2022). A new semi-analytic solution of fractional sixth order Drinfeld-Sokolov-Satsuma-Hirota equation. Numerical Methods for Partial Differential Equations, 38(3), 372–389.

72. Mainardi, F. (1996). The fundamental solutions for the fractional diffusion-wave equation. Applied Mathematics Letters, 9(6), 23–28.

73. Khater, M. M., Park, C., Lee, J. R., Mohamed, M. S., Attia, R. A. (2021). Five semi analytical and numerical simulations for the fractional nonlinear space-time telegraph equation. Advances in Difference Equations, 2021(1), 227.

74. Reutskiy, S. Y. (2016). The backward substitution method for multipoint problems with linear Volterra-Fredholm integro-differential equations of the neutral type. Journal of Computational and Applied Mathematics, 296, 724–738.

75. Zhang, Y., Rabczuk, T., Lu, J., Lin, S., Lin, J. (2022). Space-time backward substitution method for nonlinear transient heat conduction problems in functionally graded materials. Computers & Mathematics with Applications, 124, 98–110.

76. Esmaeili, S., Shamsi, M., Luchko, Y. (2011). Numerical solution of fractional differential equations with a collocation method based on Müntz polynomials. Computers & Mathematics with Applications, 62(3), 918–929.

77. Mokhtary, P., Ghoreishi, F., Srivastava, H. M. (2016). The M üntz-Legendre Tau method for fractional differential equations. Applied Mathematical Modelling, 40(2), 671–684.

78. Reutskiy, S. Y., Lin, J. (2018). A semi-analytic collocation method for space fractional parabolic PDE. International Journal of Computer Mathematics, 95(6–7), 1326–1339.

79. Lin, J., Zhang, Y., Reutskiy, S. (2021). A semi-analytical method for 1D, 2D and 3D time fractional second order dual-phase-lag model of the heat transfer. Alexandria Engineering Journal, 60(6), 5879–5896.

80. Lin, J., Bai, J., Reutskiy, S., Lu, J. (2023). A novel RBF-based meshless method for solving time-fractional transport equations in 2D and 3D arbitrary domains. Engineering with Computers, 39(3), 1905–1922.

81. Diethelm, K., Ford, J. (2002). Numerical solution of the Bagley-Torvik equation. BIT Numerical Mathematics, 42, 490–507.

82. Bhrawy, A. H., Zaky, M. A. (2016). Shifted fractional-order Jacobi orthogonal functions: Application to a system of fractional differential equations. Applied Mathematical Modelling, 40(2), 832–845.

83. Jin, B., Lazarov, R., Liu, Y., Zhou, Z. (2015). The Galerkin finite element method for a multi-term time-fractional diffusion equation. Journal of Computational Physics, 281, 825–843.

84. Hadhoud, A. R., Abd Alaal, F. E., Abdelaziz, A. A., Radwan, T. (2021). Numerical treatment of the generalized time-fractional Huxley-Burgers’ equation and its stability examination. Demonstratio Mathematica, 54(1), 436–451. https://doi.org/10.1515/dema-2021-0040

85. Ferrás, L. L., Ford, N., Morgado, M. L., Rebelo, M. (2021). High-order methods for systems of fractional ordinary differential equations and their application to time-fractional diffusion equations. Mathematics in Computer Science, 15, 535–551. https://doi.org/10.1007/s11786-019-00448-x

86. Reutskiy, S., Fu, Z. J. (2018). A semi-analytic method for fractional-order ordinary differential equations: Testing results. Fractional Calculus and Applied Analysis, 21, 1598–1618. https://doi.org/10.1515/fca-2018-0084

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.