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

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)