iconOpen Access



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

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

1 A. Pidhornyi Institute of Mechanical Engineering Problems of NAS of Ukraine, Kharkiv, 61046, Ukraine
2 College of Mechanics and Materials, Hohai University, Nanjing, 211100, China
3 Nanjing Hydraulic Research Institute, Nanjing, 210029, China
4 Shigatse Everest Urban Investment and Development Group Co, Ltd., 857000, China

* Corresponding Authors: Yuhui Zhang. Email: email; Jun Lu. Email: email

(This article belongs to this Special Issue: New Trends on Meshless Method and Numerical Analysis)

Computer Modeling in Engineering & Sciences 2024, 139(2), 1583-1612. https://doi.org/10.32604/cmes.2023.044878


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.


1  Introduction

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


is proposed in this paper. Here: n1<α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]:


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 zn or z𝒩0 and z>n1,(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 [35], 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 [1114].

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


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


is a spatial differential operator of the second order defined for 0x1. 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]


– the time-fractional telegraph equation [1618]


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


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


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 [4244]. 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 [5456]. Many works focus on the error analysis of the FE methods such as [5759].

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 [6063] and the methods based on the integration [6467]. 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(n1)Υ(0)=Υn1.(10)

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




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)




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

P(0)=0, tP(0)=0,, t(n1)P(0)=0.(15)

Let us rewrite the system in the form:


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


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


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 [7880]. 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:


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


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


satisfies the FODE


Because n1<αn the function Φm(t) satisfies zero ICs


Therefore, the linear combination


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:




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


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


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


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


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





The collocation matrix C^ contains NcM blocks


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),


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)+Θn1(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)


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


Let us define the new function




This function satisfies the equation


under the BCs



and zero ICs




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


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


So, (44) is the analytical expression.

Let us define the function ug(x,t)


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


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




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



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]:


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


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


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


As a result we get the linear system



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)


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


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:


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:


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


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(n1)Υ(0)=0.(65)

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


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


3.1 Nonlinear Problem

Let us consider TFPDE (4) with a nonlinear term


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)


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.



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)).



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).


4.1 Numerical Experiments for Systems of FODEs

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


with the general exact solution


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


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]


with the general exact solution



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


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]


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


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


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




As it follows from (82), (83)


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=σ(m1)+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


The ICs are




The vector F(t) corresponds to the exact solution


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:




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+13, 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]


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

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 2N26. 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=210 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


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 109 with M=10 and N=20. It should also be noted that with a small number of N=4, the computed errors are around 103, 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:


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


The boundary conditions are


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

Case 2: Consider the following equation:


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


The boundary conditions are


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

Case 3: Consider the following equation:


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


The boundary conditions are


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 106. 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.



images images

Figure 5: The domain absolute errors


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


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


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(1x)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:


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)… . ,Υn11=t(n1)Υ(T),

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


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


with the ICs

Υ2(0)=Υ01,tΥ2(0)=Υ11… . ,t(n1)Υ2(0)=Υn11.

Then, the vectors Υ2(T), tΥ2(T),, t(n1)Υ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.


  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  6. Qin, S., Liu, F., Turner, I., Vegh, V., & Yu, Q. (2017). Multi-term time-fractional Bloch equations and application in magnetic resonance imaging. Journal of Computational and Applied Mathematics, 319, 308-319. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar]
  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. [Google Scholar] [PubMed] [CrossRef]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar] [PubMed]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  31. Kulish, V. V., & Lage, J. L. (2002). Application of fractional calculus to fluid mechanics. Journal of Fluids Engineering, 124(3), 803-806. [Google Scholar]
  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. [Google Scholar]
  35. Odibat, Z. M. (2010). Analytic study on linear systems of fractional differential equations. Computers & Mathematics with Applications, 59(3), 1171-1183. [Google Scholar]
  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. [Google Scholar]
  37. Bailey, D. H., & Swarztrauber, P. N. (1991). The fractional Fourier transform and applications. SIAM Review, 33(3), 389-404. [Google Scholar]
  38. Das, S. (2009). Analytical solution of a fractional diffusion equation by variational iteration method. Computers & Mathematics with Applications, 57(3), 483-487. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  41. Li, C., & Zeng, F. (2012). Finite difference methods for fractional differential equations. International Journal of Bifurcation and Chaos, 22(4), 1230014. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  66. Li, X., & Li, S. (2021). A fast element-free Galerkin method for the fractional diffusion-wave equation. Applied Mathematics Letters, 122, 107529. [Google Scholar]
  67. Li, X. (2023). A stabilized element-free Galerkin method for the advection–diffusion–reaction problem. Applied Mathematical Letters, 146, 108831. [Google Scholar] [CrossRef]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  72. Mainardi, F. (1996). The fundamental solutions for the fractional diffusion-wave equation. Applied Mathematics Letters, 9(6), 23-28. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  81. Diethelm, K., & Ford, J. (2002). Numerical solution of the Bagley-Torvik equation. BIT Numerical Mathematics, 42, 490-507. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]

Cite This Article

Reutskiy, S., Zhang, Y., Lu, J., Pubu, C. (2024). A Novel Method for Linear Systems of Fractional Ordinary Differential Equations with Applications to Time-Fractional PDEs. CMES-Computer Modeling in Engineering & Sciences, 139(2), 1583–1612.

cc 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.
  • 435


  • 163


  • 0


Share Link