Loading [MathJax]/jax/output/SVG/jax.js
Research article

A novel numerical approach for solving fractional order differential equations using hybrid functions

  • Received: 13 January 2021 Accepted: 15 March 2021 Published: 22 March 2021
  • MSC : 26A33, 44A45

  • This article presents a novel numerical method for seeking the numerical solutions of fractional order differential equations using hybrid functions consisting of block-pulse functions and Taylor polynomials. The fractional integrals operational matrix of the hybrid function is conducted through projecting the hybrid functions onto block-pulse functions. Then, the fractional order differential equations are converted to a set of algebraic equations via the derived operational matrix. Then, the numerical solutions are obtained via solving the algebraic equations. Moreover, we perform error analysis of the algorithm and gives the upper bound of absolute error. Finally, numerical examples are given to show the effectiveness of the proposed method.

    Citation: Hailun Wang, Fei Wu, Dongge Lei. A novel numerical approach for solving fractional order differential equations using hybrid functions[J]. AIMS Mathematics, 2021, 6(6): 5596-5611. doi: 10.3934/math.2021331

    Related Papers:

    [1] Mohamed Obeid, Mohamed A. Abd El Salam, Mohamed S. Mohamed . A novel generalized symmetric spectral Galerkin numerical approach for solving fractional differential equations with singular kernel. AIMS Mathematics, 2023, 8(7): 16724-16747. doi: 10.3934/math.2023855
    [2] Imran Talib, Md. Nur Alam, Dumitru Baleanu, Danish Zaidi, Ammarah Marriyam . A new integral operational matrix with applications to multi-order fractional differential equations. AIMS Mathematics, 2021, 6(8): 8742-8771. doi: 10.3934/math.2021508
    [3] Chang Phang, Abdulnasir Isah, Yoke Teng Toh . Poly-Genocchi polynomials and its applications. AIMS Mathematics, 2021, 6(8): 8221-8238. doi: 10.3934/math.2021476
    [4] Zahra Pirouzeh, Mohammad Hadi Noori Skandari, Kamele Nassiri Pirbazari, Stanford Shateyi . A pseudo-spectral approach for optimal control problems of variable-order fractional integro-differential equations. AIMS Mathematics, 2024, 9(9): 23692-23710. doi: 10.3934/math.20241151
    [5] P. Pirmohabbati, A. H. Refahi Sheikhani, H. Saberi Najafi, A. Abdolahzadeh Ziabari . Numerical solution of full fractional Duffing equations with Cubic-Quintic-Heptic nonlinearities. AIMS Mathematics, 2020, 5(2): 1621-1641. doi: 10.3934/math.2020110
    [6] Ismail Gad Ameen, Dumitru Baleanu, Hussien Shafei Hussien . Efficient method for solving nonlinear weakly singular kernel fractional integro-differential equations. AIMS Mathematics, 2024, 9(6): 15819-15836. doi: 10.3934/math.2024764
    [7] Shazia Sadiq, Mujeeb ur Rehman . Solution of fractional boundary value problems by ψ-shifted operational matrices. AIMS Mathematics, 2022, 7(4): 6669-6693. doi: 10.3934/math.2022372
    [8] Fawaz K. Alalhareth, Seham M. Al-Mekhlafi, Ahmed Boudaoui, Noura Laksaci, Mohammed H. Alharbi . Numerical treatment for a novel crossover mathematical model of the COVID-19 epidemic. AIMS Mathematics, 2024, 9(3): 5376-5393. doi: 10.3934/math.2024259
    [9] Najat Almutairi, Sayed Saber . Chaos control and numerical solution of time-varying fractional Newton-Leipnik system using fractional Atangana-Baleanu derivatives. AIMS Mathematics, 2023, 8(11): 25863-25887. doi: 10.3934/math.20231319
    [10] Samia Bushnaq, Kamal Shah, Sana Tahir, Khursheed J. Ansari, Muhammad Sarwar, Thabet Abdeljawad . Computation of numerical solutions to variable order fractional differential equations by using non-orthogonal basis. AIMS Mathematics, 2022, 7(6): 10917-10938. doi: 10.3934/math.2022610
  • This article presents a novel numerical method for seeking the numerical solutions of fractional order differential equations using hybrid functions consisting of block-pulse functions and Taylor polynomials. The fractional integrals operational matrix of the hybrid function is conducted through projecting the hybrid functions onto block-pulse functions. Then, the fractional order differential equations are converted to a set of algebraic equations via the derived operational matrix. Then, the numerical solutions are obtained via solving the algebraic equations. Moreover, we perform error analysis of the algorithm and gives the upper bound of absolute error. Finally, numerical examples are given to show the effectiveness of the proposed method.



    Fractional order differential equations (FODEs) are the extension of classical integer order differential equations (IODEs) using the concept of fractional calculus (FC) [1]. Compared to integer calculus (FC), FC is non-local and long-time memory, which enables FODEs be an alternative powerful mathematical tool for modeling the nature of many physical systems and real processes. For example, the relaxation modulus in viscoelastic material exhibits power-law behavior, which indicates that in the relaxation process the material has history memory. In this case, the relationship between stress and strain can be modeled using a FODE as σ(t)=Eταdαε(t)dtα [2]. Also, the Warburg impedance has a fractional-power-law dependence on angular frequency, which is modeled as Z(s)=As1/2 [3].

    In the last decades, linear or nonlinear FODEs have emerged in many scientific and engineering fields such as neural networks [4,5], chaos [6], diffusion process [7], fluid flows in porous media [8] and optimal control [9]. Now, more and more applications of FODEs can be found in practice, there is a great need of finding the solution of FODEs. However, it is not an easy thing to obtain the exact or analytic solution for most FODEs, pursuing numerical solution of FODEs is an important thing. In the past, the finite difference method [10], predictor-corrector method [11], Adomian decomposition method [12], variational iteration method [13] and Homotopy analysis method [14] had been proposed to solve FODES.

    Though those methods are efficient and can provide good approximation to the real solution, however, they possess high computational complexities, which hampers its real application. So, developing efficient and effective methods for solving FODEs becomes an urgent task. Recently, the operational matrices had been widely adopted by researchers to solve different kinds of FODEs [15]. The core idea of operational matrices based methods is to transform the FODEs to a set of algebraic equations. As a result, the problem is dramatically simplified. In the literature, the operational matrices of fractional operators of block pulse functions (BPFs) [16], Bernstein polynomials [17], Bernoulli polynomials [18], Taylor polynomials [19] had been developed and adopted by several scholars to solve FODEs numerically. Among the above different polynomials or functions, the BPFs' fractional integral operational matrix is upper triangular matrix, in which the kth row can be obtained by shifting the (k1)th row to the right. This an appealing property which enable one to solve the set of algebraic equations more easily and drastically reduce the computation burden involved in the computation of matrix inverse. However, BPFs are not enough smooth, more BPFs are required if high approximation accuracy is desired, and in turn, the dimension of operational matrix increases.

    In recent years, the hybrid functions, which mainly combines BPFs or Haar function with other polynomials, had been widely adopted by researchers to solve different kind of FODEs. For example, in [20], the hybrid of BPFs and Bernoulli polynomials is adopted to find the numerical solution of nonlinear fractional integro-differential equations. In [21], the hybrid of Legendre polynomial and Haar function, called Legendre wavelet is used to solve distributed order FODEs. In [22], Chebyshev wavelets is adopted to solve nonlinear variable order FODEs. One of the advantage of hybrid functions is that they are piecewise polynomial instead of piecewise constant in each interval as BPFs. Thus, hybrid functions are smoother than BPFs, which leads to a more accurate approximation to a function than BPFs if an equal number of basis functions are used. Therefore, one can obtain more accurate solution for a FODEs using hybrid functions than BPFs.

    Observing the aforementioned facts, this paper presents a new numerical method for solving FODEs based on hybrid of BPFs and Taylor polynomials (HBT). The fractional integrals operational matrix of HBT is derived through projecting the HBT functions onto a set of BPFs. Then, the FODEs to be solved is transformed into a set of algebraic equations using the derived operational matrix. Through solving the algebraic equations, one can obtain the numerical solution of FODEs.

    The organization of this paper is arranged as follows. Some basics about fractional calculus are briefly reviewed in Section 2. In Sections 3, the basic formulation of HBT functions is given, and the calculation of fractional integral operational matrix of HBT is given in Section 4. The error analysis of the proposed method is given in Section 5. In Section 6, numerical experiments are presented to verify the effectiveness of the proposed method. In the last, the conclusive remarks are presented in Section 7.

    Unlike integer calculus, the definition of FC is not unique. There are several definitions for fractional integration and derivatives. Among these definitions, the widely used one are the Riemann-Liouville (R-L) definition and the Caputo definition.

    In R-L definition, the fractional integral of function f(t) is given as [23]

    (0Iαtf)(t)=1Γ(q)t0(ts)α1f(s)ds, (2.1)

    where α>0 is the order, 0 and t are the low and upper limits of fractional integrals, Γ() denotes the Gamma function. The fractional derivative of a function f(t) in the sense of Caputo definition is

    (0Dαtf)(t)=1Γ(nα)t0f(n)(s)(ts)α+1nds,n1<α<n,nN. (2.2)

    R-L fractional integral and Caputo fractional derivative has the following relationship

    (0Dαt0Itαf)(t)=f(t) (2.3)

    and

    (0Iαt0Dtαf)(t)=f(t)n1k=0f(k)(0)tkk!. (2.4)

    Since the derivation is the inverse operation of integration, the derivation of a function can easily be achieved from the result of integration.

    The Taylor polynomials defined on the interval [0,T] is [24]

    Tj(t)=tj,j=0,1,2, (3.1)

    A set of BPFs ψi(t),i=1,2,N is given as follows [25]

    ψi(t)={1,i1NTt<iNT,0,otherwise (3.2)

    The BPFs ψi(t) are disjoint and orthogonal, that is to say

    ψi(t)ψl(t)={0,il,ψi(t),i=l, (3.3)
    10ψi(t)ψl(t)dt={0,il,1N,i=l. (3.4)

    A function f(t) in Cn+1[0,T) can be expanded into an N-term BPF series as

    f(t)=[c1,c1,,cN]Ψ(t)=CTΨ(t), (3.5)

    where Ψ(t)=[ψ1(t),ψ2(t),,ψN(t)]T, the constant coefficients ci in Eq (3.2) are given by

    ci=NiNTi1NTf(t)dt. (3.6)

    On the interval [0,T), the hybrid of BPFs and Taylor polynomials are defined as [24]

    hij(t)={Tj(NTti+1),i1NTt<iNT,0,otherwise, (3.7)

    where i and j are the order of BPFs and the degree of Taylor polynomials respectively.

    Let H=Cn+1[0,T),{h10(t),h11(t),,hN(M1)(t)}H, be a set of HBT functions, X=span{h10(t),,h1(M1)(t),h20(t),,h2(M1)(t),,hN0(t),,hN(M1)(t)}, and f(t) be a function in H. Since X is a finite dimensional subspace of H, f(t) has the unique best approximation in X, such as f(t)X, that is

    x(t)X,f(t)f(t)∥≤∥f(t)x(t).

    Since f(t)X, there exist the unique coefficients c10,c11,,cN(M1) such that

    f(t)f(t)=Ni=1M1j=0cijhij(t)=CTH(t), (3.8)

    where T denotes transpose of vector or matrix, C and H(t) are N×M column vectors

    C=[c10,,c1(M1),c20,,c2(M1),,cN0,,cN(M1)]T, (3.9)

    and

    H(t)=[h10(t),,h1(M1)(t),h20(t),,h2(M1)(t),,hN0(t),,hN(M1)(t)]T. (3.10)

    To evaluate C, firstly let

    fnm=T0f(t)hnm(t)dt

    Using Eq (3.8) one gets

    fnm=Ni=1M1j=0cijT0hij(t)hnm(t)dt=Ni=1M1j=0cijdnmij,n=1,2,.N,m=0,1,M1,

    where

    dnmij=T0hij(t)hnm(t)dt.

    Therefore,

    fnm=CT[dnm10,,dnm1(M1),,dnm20,,dnm2(M1),,,dnmN0,,,dnmN(M1)]T

    or

    FT=CTD

    where

    F=[f10,,f1(M1),f20,,f2(M1),,fN0,,fN(M1)]T

    and

    D=[dnmij]

    is a matrix of order NM×NM and is given by

    D=T0H(t)HT(t)dt. (3.11)

    For HBT functions D has the following form

    D=(˜D1000˜D2000˜DN) (3.12)

    where

    ˜Di=1NT0T(t)TT(t)dt,i=1,2,,N (3.13)

    and T(t)=[T0(t),T1(t),,TM1(t)]T. Hence, CT in Eq (3.8) is given by

    CT=FTD1. (3.14)

    Here, the fractional integral operational matrix of HBT is derived. To simplify the derivation process, the HBT functions are first projected onto a set of BPFs, and the projection matrix is calculated. The HBT operational matrix is obtained by using matrix operations of the projection matrix and BPF operational matrix.

    The R-L fractional integration operator of hybrid function vector H(t) can be written as

    (0IαtH)(t)PαH(t), (4.1)

    where matrix Pα is the fractional integral operational matrix of HBT. To simplify the derivation process, the hybrid functions are projected onto a set of BPFs ψi(t),i=1,2,,m, and m=N×M. Specifically, the hybrid function vector defined in Eq (3.10) can be expressed as

    H(t)=ΦΨ(t) (4.2)

    where Ψ(t)=[ψ1(t),ψ2(t),,ψm(t)]T and Φ is the projection matrix which transform the hybrid functions onto BPFs.

    In general, for arbitrary M and N, one has

    Φ={Φ1=[aij]M×NM,0t<1NT,Φ2=[bij]M×NM,1NTt<2NT,ΦN=[dij]M×NM,N1NTt<T. (4.3)

    The expressions of aij, bij and dij are given in Appendix A. As an example, let N=2,M=3 and T=1, one has

    H(t)=[h10(t),h11(t),h12(t),h20(t),h21(t),h22(t)]T (4.4)
    Ψ(t)=[ψ1(t),ψ2(t),ψ3(t),ψ4(t),ψ5(t),ψ6(t)]T. (4.5)

    the function vector (4.4) can be expressed as

    h10(t)=ψ1(t)+ψ2(t)+ψ3(t)h11(t)=16ψ1(t)+12ψ2(t)+56ψ3(t)h12(t)=127ψ1(t)+727ψ2(t)+1927ψ3(t)},0t<12,
    h10(t)=ψ4(t)+ψ5(t)+ψ6(t)h11(t)=16ψ4(t)+12ψ5(t)+56ψ6(t)h12(t)=127ψ4(t)+727ψ5(t)+1927ψ6(t)},12t<1.

    In this case,

    Φ={Φ1=[aij]3×6,0t<12,Φ2=[bij]3×6,12t<1.

    where

    Φ1=(1110001/61/25/60001/277/2719/27000),
    Φ2=(0001110001/61/25/60001/277/2719/27).

    Hence, we have

    b14=a11,b15=a12b16=a13,b24=a21,b25=a22,b26=a23,b34=a31,b35=a32,b36=a33.

    The projection matrix is derived in detail in Appendix A and omitted here. It is easy to see from the above example that the obtained projection matrix is in block-diagonal form. This elegant property is useful for simplifying the computation of the operational matrices.

    From Eq (4.1) and Eq (4.2), one has

    (0IαtH)(t)=0Iαt(ΦΨ)(t)=Φ(0IαtΨ)(t). (4.6)

    According to Ref.[25],

    (0IαtΨ)(t)FαΨ(t), (4.7)

    where

    Fα=(Tm)α1Γ(α+2)=(1ξ1ξ2ξm101ξ1ξm2001ξm30000001) (4.8)

    is the BPFs operational matrix of fractional integration with ξk=(k+1)α+12kα+1+(k1)α+1. Substituting Eq (3.6) into Eq (4.6), one can get

    (0IαtH)(t)=ΦFαΨ(t)=ΦFαΦ1H(t). (4.9)

    Therefore,

    Pα=ΦFαΦ1. (4.10)

    Denote the HBT approximation to function f(t) at the level l as f(t)=Ni=1M1j=0cijhij(t), where l is determined by N and M(l=lN,M). We replace f(t) with f(t) in Eq (2.1) and call the resulting integral HBT approximation of the α order R-L fractional integral of f(t).

    0Iαtf(t)=1Γ(α)t0(ts)α1f(s)ds,

    then the absolute error between Iαf(t) and Iαf(t) is

    εl=|0Iαtf(t)0Iαtf(t)|.

    Theorem 1. Let f(t)CM[0,1), f(t)=Ni=1M1j=0cijhij(t) is the best approximation of f(t) in X, then

    1). |f(t)f(t)|TMM!NM|f(M)(t)|,t[i1NT,iNT),i=1,2,,N.

    2). εlKTα+MM!NMΓ(α+1),f(M)(t)∣≤K,K is a finite positive value.

    Proof. 1). Using Taylor formula, f(t) can be approximated in the ith interval [i1NT,iNT) with

    fM1i(t)=f(i1N)+f(i1N)(ti1N)+f(i1N)(ti1N)22++f(M1)(i1N)(ti1N)M1(M1)!,

    the truncation error is

    f(t)fM1i(t)=(ti1N)MM!f(M)(ξi),

    where ξi lies between i1N and t. Since the best approximation is unique [26], for all t[i1NT,iNT), the absolute error between f(t) and f(t) is given as

    |f(t)f(t)||f(t)fM1i(t)|TMM!NM|f(M)(t)|.

    2). The absolute error between Iαf(t) and Iαf(t) is

    εl=|0Iαtf(t)0Iαtf(t)|=1Γ(α)t0(ts)α1|f(s)f(s)|ds=1Γ(α)[ir=1rNTr1NT(ts)α1|f(s)f(s)|ds+tiNT(ts)α1|f(s)f(s)|ds]1Γ(α)[ir=1rNTr1NT(ts)α1TMM!NM|f(M)(s)|ds+tiNT(ts)α1TMM!NM|f(M)(s)|ds],

    according to the assumption f(M)(t)∣≤K, the absolute error between Iαf(t) and Iαf(t) can be estimated as

    εl1Γ(α)(TMKM!NM)[ir=1rNTr1NT(ts)α1ds+tiNT(ts)α1ds]KTα+MM!NMΓ(α+1),

    where K is a finite positive value.

    In order to verify the maximum absolute error arised by HBT approximation is smaller than the upper bound in theory given in Theorem 1, the function f(t)=t is selected as an example. The analytic expression of α order fractional integral of f(t)=t is given as

    (0Iαt)(t)=Γ(2)Γ(α+2)tα+1.

    Using Eq (4.9) with α=0.5,T=1 and N=2,M=3, the HBT estimation of (0Iαt)(t) is obtained as

    (0Iαt)(t)=[0.001067,0.112147,0.157513,0.267224,0.404722,0.082083]H(t),

    Using BPFs, (0Iαt)(t) is approximated as [27]

    (0Iαt)(t)=[0.033642,0.128795,0.269961,0.443942,0.645261,0.870557]Ψ(t).

    Table 1 presents the absolute errors obtained by the approximation of BPFs and HBT functions. The piecewise-polynomial property of HBT functions lead to a more accurate approximation of the fractional integral even with the equal number of basis functions as BPFs. Therefore, HBT approximation of R-L fractional integral is effective.

    Table 1.  Absolute errors using BPFs and HBT functions.
    t |Iαf(t)IαBPFt| |Iαf(t)IαHBTt|
    0 0.0336 0.0011
    0.2 0.0611 0.0038
    0.4 0.0797 0.0013
    0.6 0.0943 0.0018
    0.8 0.1070 0.0013

     | Show Table
    DownLoad: CSV

    To show the effectiveness of our presented method, the fractional integral operational matrix of HBT is used to solve several FODEs. The solutions obtained by our proposed method are compared with their exact solutions, those by block-pulse operational matrix (BPOM) method and fractional Taylor operational matrix (FTOM) method [19].

    A fractional Riccati equation as [28,29]

    Dαx(t)=2x(t)x2(t)+1,0<α1, (6.1)

    subject to the initial condition x(0)=0 is considered. Let

    Dαx(t)=CTH(t), (6.2)

    together with the initial condition, one has

    x(t)=CTPαH(t), (6.3)

    Since H(t)=ΦΨ(t), from Eq (6.3) one has

    x(t)=CTΦFαΨ(t), (6.4)

    Let

    CTPαΦ=[a1,a2,,am], (6.5)

    and using Eqs (3.3) and (3.4) one can obtain

    [x(t)]2=[a1ψ1(t)+a2ψ2(t)++amψm(t)]2=[a21,a22,,a2m]Ψ(t). (6.6)

    Substituting Eqs(4.2), (6.2), (6.4) and (6.6) into Eq (6.1), one has

    CTΦΨ(t)2CTP(α)ΦΨ(t)+[a21,a22,,a2m]Ψ(t)[1,1,,1]Ψ(t)=0. (6.7)

    Equation (6.7) is a set of algebraic equations. In this paper, the Matlab function fsolve is used to solve Eq (6.7). The numerical solution, for N=4,M=6 is shown in Figure 1. When α=1, the exact solution of this equation is

    x(t)=1+2tanh(2t+12ln(212+1))
    Figure 1.  Comparison of x(t) with N=4,M=6, α=0.7,0.8,0.9,1, with exact solution in Example 1.

    Table 2 shows comparisons of the approximate solutions at different values of t obtained by the present method with N=32,M=6, BPOM method with N=32, FTOM method with M=6 and the exact solution for α=1. In spite of the large amount of calculation, the solutions obtained by the present method are the same as exact solutions.

    Table 2.  Numerical results in Example 1.
    t BPOM method FTOM method Our method Exact solution
    0 0 0 0 0
    0.1 0.129822 0.110273 0.110295 0.110295
    0.2 0.259644 0.241973 0.241977 0.241977
    0.3 0.414392 0.395098 0.395105 0.395105
    0.4 0.594065 0.567813 0.567812 0.567812
    0.5 0.773738 0.755979 0.756014 0.756014
    0.6 0.973682 0.953510 0.953566 0.953566
    0.7 1.173626 1.152935 1.152949 1.152949
    0.8 1.361436 1.346364 1.346364 1.346364
    0.9 1.537112 1.1526844 1.526911 1.526911
    Time 0.99s 1.55s 8.25s -

     | Show Table
    DownLoad: CSV

    On the other hand, the numerical solutions obtained by the present method with N=4,M=6 for α=0.7,0.8,0.9,1 are shown in Figure 1. It is demonstrate that the approximate solution obtained by the present method is in good agreement with the exact solution when α=1. Moreover, as indicated in [28], the solution continuously depends on the time-fractional derivative.

    A fractional differential equation as

    aD2x(t)+bDα2x(t)+cDα1x(t)+e[x(t)]3=f(t),0<α1,1<α22 (6.8)

    and

    f(t)=2at+2bΓ(4α2)t3α2+2cΓ(4α1)t3α1+e[13t3]3

    subject to x(0)=x(0)=0 is considered [30]. The exact solution of this equation is x(t)=13t3.

    Let

    D2x(t)=CTH(t), (6.9)

    and take the initial states into consideration, one has

    Dα2x(t)=CTP2α2H(t), (6.10)
    Dα1x(t)=CTP2α1H(t), (6.11)

    Since H(t)=ΦΨ(t), from Eq (6.9) we have

    x(t)=CTP2ΦΨ(t). (6.12)

    Let

    CTP2Φ=[a1,a2,,am], (6.13)

    then

    [x(t)]3=[a31,a32,,a3m]Ψ(t). (6.14)

    Expanding f(t) onto HBT functions, one has

    f(t)=fTH(t), (6.15)

    where fT is a known constant vector. Substituting Eq (6.9)–(6.11) and Eq (6.14) into Eq (6.8), together with H(t)=ΦΨ(t), then

    CTΦΨ(t)+CTP2α2ΦΨ(t)+CTP2α1Φ(t)+[a31,a32,,a3m]Ψ(t)fTΦΨ(t)=0. (6.16)

    Eq (6.16) is a set of algebraic equations, which is solved via Matlab function fsolve. Figure 2 shows the numerical solution obtained by the present method with N=4,M=6 when a=1,b=1,c=1,e=1,α1=0.333,α2=1.234. It can be seen that the numerical solution is almost the same as the exact solution.

    Figure 2.  Numerical solution with N=4,M=6 and exact solution in Example 2.

    The comparison of absolute errors in x(t) obtained by the BPOM method with N=4, FTOM method with M=6 and present method with N=4,M=6 at given points for different values of N and M are shown in Table 3. The computational results illustrate that, compared with BPOM method and FTOM method, the present method can get more accurate approximate solutions although it has a large amount of calculation. Moreover, for the present method, the error is smaller and smaller when N and M increase. Therefore for higher accuracy of the approximation, lager N and M are recommended.

    Table 3.  Comparison of absolute errors for different methods in Example 2.
    t BPMO method FTOM method Our method
    (N=1,M=6)
    Our method
    (N=2,M=6)
    Our method
    (N=4,M=6)
    0 0 0 0 0 0
    0.1 8.07E-4 1.48E-9 6.09E-12 1.40E-15 8.56E-16
    0.2 0.0019 9.61E-9 1.04E-11 3.02E-15 1.55E-15
    0.3 0.0019 1.76E-8 1.27E-11 5.19E-16 3.61E-14
    0.4 0.0026 1.93E-8 1.56E-11 9.89E-16 4.89E-14
    0.5 0.0027 1.68E-8 1.82E-11 2.93E-14 2.73E-14
    0.6 0.0026 1.64E-8 1.98E-11 9.45E-13 2.80E-13
    0.7 0.0030 2.12E-8 2.23E-11 1.15E-12 3.83E-13
    0.8 0.0023 2.63E-8 2.61E-11 8.05E-13 2.75E-13
    0.9 0.0025 2.51E-8 2.65E-11 8.01E-13 1.01E-13
    Time 0.957s 1.014s 1.091s 2.103s 9.406s

     | Show Table
    DownLoad: CSV

    In this paper, a new numerical method is presented to solve FODEs by using the hybrid functions consisting of BPFs and Taylor polynomials. The fractional integral operational matrix of the hybrid functions is derived. The FODEs to be solved is transformed into a set of algebraic equations via the operational matrix. Finally, the numerical solutions of FODEs are get by solving the algebraic equations. Illustrative examples verify that the proposed algorithm can get more accurate numerical solutions of FODEs than BPOM method and FTOM method although it has a large amount of calculation.

    This work is partially supported by Public Welfare Project of Zhejiang Province of China (LGN20C050002).

    The authors declare no conflict of interest.

    In this section, the derivation process of the projection matrix from HBT functions to BPFs is given. Note that Eq (3.10) is the HBT function vector. First, elements h1j(t) of H(t) are considered and can be expanded into BPFs as follows:

    h1j(t)mk=1a(j+1)kψk(t) (A.1)

    where j=0,1,,(M1) and m=N×M. The coefficient a(j+1)k can be computed according to Eq (3.6)

    a(j+1)k=h1j(t),ψk(t)ψk(t),ψk(t)=T0h1j(t)ψk(t)dtT0ψk(t)ψk(t)dt (A.2)

    Since h1j(t) is nonzero when 0t<TN, and ψk(t) is defined on the interval [k1NMT,kNMT),h1j(t)ψk(t) is nonzero only when 1kM. Observing this fact, Eq (A.2) can be written as

    a(j+1)k={MNTkMNTk1MNT(NTt)jdt,1kM,0,otherwise.

    After some manipulations, it is easy to obtain

    a(j+1)k={1j+11Mj[kj+1(k1)j+1],1kM,0,otherwise.

    Define

    pk=kj+1(k1)j+1 (A.3)

    Then

    [a(j+1)1,a(j+1)2,,a(j+1)m]=1j+11Mj[p1,p2,,pM,(N1)×M terms00] (A.4)

    Thus, when 0t<1N

    [h10(t),h11(t),,h1m(t)]T=Φ1ΨT(t), (A.5)

    where

    Φ1=(a11a12a1M00a21a22a2M00aM1aM2aMM00), (A.6)

    similarly, when TNt<2NT,h2j(t),j=0,1,,(M1) can be also expanded into BPFs, and

    h2j(t)mk=1b(j+1)kψk(t) (A.7)

    where

    b(j+1)k={MNTkMNTk1MNT(NTt1)jdt,M+1k2M,0,otherwise.

    Let x=tTN, one has

    b(j+1)k={MNTkMMNTkM1MNT(NTt)jdt=a(j+1)(kM),M+1k2M,0,otherwise.

    Rewriting b(j+1)k into vector form,

    [b(j+1)1,b(j+1)2,,b(j+1)m]=1j+11Mj[M terms0,,0,p1,p2,,pM,(N2)×M terms00]. (A.8)

    Then

    Φ2=(000a11a12a1M00000a21a22a2M00000aM1aM2aMM00), (A.9)

    similarly, when N1NTt<T,

    hNj(t)mk=1d(j+1)kψk(t) (A.10)

    where j=0,1,,(m1) and m=N×M.

    And finally, take the same process as above, the matrix ΦN can be obtained, where ΦN=[d(j+1)k]M×NM

    [d(j+1)1,d(j+1)2,,d(j+1)m]=1j+11Mj[(N1)×M terms0,,0,p1,p2,,pM]. (A.11)

    Then

    ΦN=(000a11a12a1M000a21a22a2M000aM1aM2aMM), (A.12)

    Therefore

    Φ={Φ1=[aij]M×NM,0t<TN,Φ2=[bij]M×NM,TNt<2NT,ΦN=[dij]M×NM,N1NTt<T. (A.13)


    [1] V. Lakshmikantham, A. Vatsala, Basic theory of fractional differential equations, Nonlinear Anal. Theor., 69 (2008), 2677–2682. doi: 10.1016/j.na.2007.08.042
    [2] F. C. Meral, T. J. Royston, R. Magin, Fractional calculus in viscoelaticity: An experimental study, Commun. Nonlinear Sci., 15 (2010), 939–945. doi: 10.1016/j.cnsns.2009.05.004
    [3] J. Wang, Realizations of generalized warburg impedance with RC ladder networks and transmission lines, J. Electrochem. Soc., 134 (1987), 1915. doi: 10.1149/1.2100789
    [4] L. L. Huang, J. H. Park, G. C. Wu, Z. W. Mo, Variable-order fractional discrete-time recurrent neural networks, J. Comput. Appl. Math., 370 (2020), 112633. doi: 10.1016/j.cam.2019.112633
    [5] G. C. Wu, M. Luo, L. L. Huang, S. Banerjee, Short memory fractional differential equations for new memristor and neural network design, Nonlinear Dynam., 100 (2020), 3611–3623. doi: 10.1007/s11071-020-05572-z
    [6] G. C. Wu, M. Niyazi Çankaya, S. Banerjee, Fractional q-deformed chaotic maps: A weight function approach, Chaos, 30 (2020), 121106. doi: 10.1063/5.0030973
    [7] T. U. Khan, M. A. Khan, Y. M. Chu, A new generalized Hilfer-type fractional derivative with applications to space-time diffusion equation, Results Phys., 22 (2021), 103953.
    [8] M. F. El Amin, Derivation of fractional-derivative models of multiphase fluid flows in porous media, J. King Saud. Univ. Sci., 33 (2021), 101346. doi: 10.1016/j.jksus.2021.101346
    [9] H. R. Marzban, A new fractional orthogonal basis and its application in nonlinear delay fractional optimal control problems, ISA T., 2020, In press.
    [10] Y. Liu, Y. Du, H. Li, S. He, W. Gao, Finite difference/finite element method for a nonlinear time-fractional fourth-order reaction diffusion problem, Comput. Math. Appl., 70 (2015), 573–591. doi: 10.1016/j.camwa.2015.05.015
    [11] V. Daftardar-Gejji, Y. Sukale, S. Bhalekar, A new predictor-corrector method for fractional differential equations, Appl. Math. Comput., 244 (2014), 158–182. doi: 10.1016/j.amc.2014.06.097
    [12] C. Li, Y. Wang, Numerical algorithm based on adomian decomposition for fractional differential equations, Comput. Math. Appl., 57 (2009), 1672–1681. doi: 10.1016/j.camwa.2009.03.079
    [13] G. C. Wu, E. Lee, Fractional variational iteration method and its application, Phys. Lett. A, 374 (2010), 2506–2509. doi: 10.1016/j.physleta.2010.04.034
    [14] B. Ghazanfari, F. Veisi, Homotopy analysis method for the fractional nonlinear equations, J. King Saud. Univ. Sci., 23 (2011), 389–393. doi: 10.1016/j.jksus.2010.07.019
    [15] H. Dehestani, Y. Ordokhani, M. Razzaghi, Pseudo-operational matrix method for the solution of variable-order fractional partial integro-differential equations, Engineering with Computers, (2020), 1–16.
    [16] S. Najafalizadeh, R. Ezzati, A block pulse operational matrix method for solving two-dimensional nonlinear integro-differential equations of fractional order, J. Comput. Appl. Math., 326 (2017), 159–170. doi: 10.1016/j.cam.2017.05.039
    [17] M. H. Alshbool, O. Isik, I. Hashim, Fractional bernstein series solution of fractional diffusion equations with error estimate, Axioms, 10 (2021), 6. doi: 10.3390/axioms10010006
    [18] J. R. Loh, C. Phang, Numerical solution of fredholm fractional integro-differential equation with right-sided caputo's derivative using bernoulli polynomials operational matrix of fractional derivative, Mediterr J. Math., 16 (2019), 1–25. doi: 10.1007/s00009-018-1275-9
    [19] İ. Avcı, N. I. Mahmudov, Numerical solutions for multi-term fractional order differential equations with fractional Taylor operational matrix of fractional integration, Mathematics, 8 (2020), 96. doi: 10.3390/math8010096
    [20] S. Mashayekhi, M. Razzaghi, Numerical solution of nonlinear fractional integro-differential equations by hybrid functions, Eng. Anal. Bound. Elem., 56 (2015), 81–89. doi: 10.1016/j.enganabound.2015.02.002
    [21] B. Yuttanan, M. Razzaghi, Legendre wavelets approach for numerical solutions of distributed order fractional differential equations, Appl. Math. Model., 70 (2019), 350–364. doi: 10.1016/j.apm.2019.01.013
    [22] M. H. Heydari, Chebyshev cardinal wavelets for nonlinear variable-order fractional quadratic integral equations, Appl. Numer. Math., 144 (2019), 190–203. doi: 10.1016/j.apnum.2019.04.019
    [23] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Academic Press, 1998.
    [24] H. Marzban, M. Razzaghi, Solution of multi-delay systems using hybrid of block-pulse functions and Taylor series, J. Sound. Vib., 292 (2006), 954–963. doi: 10.1016/j.jsv.2005.08.007
    [25] A. Kilicman, Z. A. A. A. Zhour, Kronecker operational matrices for fractional calculus and some applications, Appl. Math. Comput., 187 (2007), 250–265 doi: 10.1016/j.amc.2006.08.122
    [26] E. Kreyszig, Introductory functional analysis with applications, New York: Wiley, 1978.
    [27] S. K. Damarla, M. Kundu, Numerical solution of multi-order fractional differential equations using generalized triangular function operational matrices, Appl. Math. Comput., 263 (2015), 189–203. doi: 10.1016/j.amc.2015.04.051
    [28] Z. Odibat, S. Momani, Modified Homotopy perturbation method: Application to quadratic Riccati differential equation of fractional order, Chaos Soliton. Fract., 36 (2008), 167–174. doi: 10.1016/j.chaos.2006.06.041
    [29] Y. Li, Solving a nonlinear fractional differential equation using chebyshev wavelets, Commun. Nonlinear Sci., 15 (2010), 2284–2292. doi: 10.1016/j.cnsns.2009.09.020
    [30] A. El-Mesiry, A. El-Sayed, H. El-Saka, Numerical methods for multi-term fractional (arbitrary) orders differential equations, Appl. Math. Comput., 160 (2005), 683–699. doi: 10.1016/j.amc.2003.11.026
  • This article has been cited by:

    1. Osama Moaaz, Ahmed E. Abouelregal, Multi-fractional-differential operators for a thermo-elastic magnetic response in an unbounded solid with a spherical hole via the DPL model, 2022, 8, 2473-6988, 5588, 10.3934/math.2023282
    2. Mohamed Karim Bouafoura, Naceur Benhadj Braiek, Suboptimal control synthesis for state and input delayed quadratic systems, 2022, 236, 0959-6518, 944, 10.1177/09596518211067476
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3427) PDF downloads(262) Cited by(2)

Figures and Tables

Figures(2)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog