J. Nig. Soc. Phys. Sci. 5 (2023) 1221 Journal of the Nigerian Society of Physical Sciences Solving fractional variable-order differential equations of the non-singular derivative using Jacobi operational matrix M. Basima, N. Senua,b,∗, A. Ahmadianc, Z. B. Ibrahimb, S. Salahshourd aInstitute for Mathematical Research, Universiti Putra Malaysia, 43400 UPM, Serdang, Malaysia bDepartment of Mathematics and statistics Universiti Putra Malaysia, 43400 UPM, Serdang, Malaysia cDecision Lab, Mediterranea University of Reggio Calabria, Reggio Calabria, Italy dFaculty of Engineering and Natural Sciences, Bahcesehir University, Istanbul, Turkey Abstract This research derives the shifted Jacobi operational matrix (JOM) with respect to fractional derivatives, implemented with the spectral tau method for the numerical solution of the Atangana-Baleanu Caputo (ABC) derivative. The major aspect of this method is that it considerably simplifies problems by reducing them to ones that can be solved by solving a set of algebraic equations. The main advantage of this method is its high robustness and accuracy gained by a small number of Jacobi functions. The suggested approaches are applied in solving non-linear and linear ABC problems according to initial conditions, and the efficiency and applicability of the proposed method are proved by several test examples. A lot of focus is placed on contrasting the numerical outcomes discovered by the new algorithm together with those discovered by previously well-known methods. Keywords: Fractional differential equations; Atangana-Baleanu Caputo; Variable order; Operational matrix. Article History : Received: 23 November 2022 Received in revised form: 31 January 2023 Accepted for publication: 14 February 2023 Published: 05 May 2023 c© 2023 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license (https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Communicated by: O. Adeyeye 1. Introduction Over the past three decades, the most focus was placed on fractional calculus study, together with its countless implemen- tations in the fields of engineering and physics. Fractional dif- ferential equations (FDEs) can effectively represent the imple- mentations of fractional calculus employed in a variety of fields, which includes signal processing, optics, statistics and proba- bility, electrochemistry of corrosion, control theory of dynam- ical systems, electrical networks, as well as chemical physics. There have been a number of important early papers on frac- tional derivatives and FDEs, as may be seen in [1, 2]. These ∗Corresponding author tel. no: +23480xxxx572 Email address: norazak@upm.edu.my (N. Senu ) publications offer a systematic explanation of fractional calcu- lus, including its uniqueness and existence, and are regarded as the introduction with respect to the FDEs and fractional deriva- tive theory. Numerous other scholars have recently focused on the findings of the initial value problem (IVP) and boundary value problem (BVP) solutions for FDEs, which can be further read in [3-5]. It is crucial to determine approximate or exact solutions to FDEs. We have trouble locating their analytical solutions for any but a small subset of these equations. Many different types of differential equations in a variety of fields in science, engineering physical and natural applications can be solved, which are extremely effective [6-11]. Other than that, many authors have been inspired to adopt these approaches for various equations because of their high accuracy and simplic- ity of usage. The collocation, Galerkin, as well as tau methods 1 FULAFIA MIS Typewriter DOI:10.46481/jnsps.2023.1221 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 2 are particular spectral methods that are more suitable and fre- quently employed. When Saadatmandi and Dehghan [12] used spectral meth- ods in solving multi-term linear as well as non-linear FDEs nu- merically, they instigated the shifted Legendre operational ma- trix with respect to fractional derivatives. To find approximate solutions for multi-term linear and non-linear FDEs, Doha et al. [13] developed a novel formula that explicitly expresses any fractional-order derivatives of shifted Chebyshev polynomials of any degree with respect to the shifted Chebyshev polynomi- als themselves. They combined this formula with tau and col- location spectral methods. Recently, Bhrawy et al. [14] treated multi-term linear FDEs having variable coefficients employ- ing a quadrature-shifted Legendre tau technique. The shifted Chebyshev operational matrix was recently presented by Doha et al. [15] and used in conjunction with spectral methods to solve multi-term linear and non-linear FDEs according to IVPs and BVPs conditions. Additionally, in [16, 17], the authors introduced the spectral tau method for the numerical solution of a few FDEs, while in [18], Pedas and Tamme created the spline collocation methods for solving FDEs. In recent years, Esmaeili and Shamsi [19] in- stigated a direct solution method for solving a particular family of fractional IVPs employing a pseudospectral method. More- over, Esmaeili et al. [20] introduced a computational method with regard to the Müntz polynomials and collocation method for the FDEs solution. Furthermore, the algorithms utilized in the current work are associated with those applied by Saa- datmandi and Dehghan [12], Doha et al. [13-15], as well as Bhrawy et al. [14] to create accurate algorithms for a variety of uses. The classical Jacobi polynomials, represented by Ju,vi (x)(i ≥ 0, u, v > −1) [21] are crucial to the study and use of spectral methods and have been widely employed in math- ematical analysis and real-world applications. The benefit of using general Jacobi polynomials is that they may be used to calculate solutions using the Jacobi parameters a and b (refer to [22]). Therefore, to generalize, it is beneficial to perform a systematic study with regards to Jacobi polynomials (u, v > −1) having general indexes. This may then be immediately imple- mented in other contexts rather than generating approximation findings for each specific indices pair. This is the reason we in- troduce the Jacobi polynomials family having indexes u, v > −1 in this work. To solve numerically linear fractional and variable or- der problems with initial conditions, this work introduces the shifted Jacobi operational matrix (JOM) of fractional deriva- tive. This method relies on the Jacobi tau method. Additionally, we present an appropriate technique for approximating the non- linear fractional and variable order IVPs on the interval [0, L] using the spectral shifted Jacobi collocation approach relying on JOM in order to determine the solution y(t). At (N − m + 1) points, the non-linear fractional and variable orders collocate. These equations produce (N + 1) non-linear algebraic equations that may be resolved employing Newton’s iterative method af- ter being combined with m initial conditions. Finally, test prob- lems are used to show how accurate the suggested algorithms are. We point out that Saadatmandi and Dehghan [12] and Doha et al. [15] introduced the two shifted Legendre and Chebyshev operational matrices, correspondingly, and that several more extremely interesting situations may be produced directly as special cases emerging from the shifted JOM [23]. As a re- sult, we were inspired to pursue the shifted Jacobi polynomials because it is the most generalized of the orthogonal polynomi- als. This paper is structured accordingly: First, we commence by going over several fundamental information about Jacobi polynomials and fractional calculus theory that is necessary for supporting our findings in Section 2. The JOM for Atangana- Baleanu Caputo (ABC) is obtained in Section 3. The spectral tau, JOM of ABC derivative, as well as collocation methods are all applied in Section 4 to solve general linear and non- linear ABC. The variable-order ABC-derivative JOM is found in section 5. The suggested methods are used in several cases in Section 6. In addition, Section 7 provides a conclusion. 2. Basic concepts and notations This section defines the Caputo derivative, CF-derivative, as well as Atangana-Baleanu Caputo (ABC)-derivative in which fractional order and variable order are concisely highlighted. 2.1. Fractional derivatives Caputo fractional-order differential equation is expressed by [1] as: CDαy(t) = 1 Γ(1 −α) ∫ t 0 y′(s)d s (t − s)α , (1) where 0 < α < 1. The most popular fractional derivative is the Caputo derivative, which is frequently used in engineering and science domains. Definition 2.1.1. To 0 < α < 1, y(t) ∈ H1(ı, ),  > ı, the fractional-order of the CF-derivative is defined by [29]: CFDαy(t) = M(α) 1 −α ∫ t 0 y′(s)e −α(t−s) 1−α d s, (2) in which M(α) denotes a normalization function. Definition 2.1.2. For 0 < α < 1, y(t) ∈ H1(ı, ),  > ı, the fractional-order of the ABC-derivative is represented by [30]: ABCDαy(t) = M(α) 1 −α ∫ t 0 y ′ (s)Eα [ −α(t − s)α 1 −α ] d s, (3) in which 0 < α < 1, M(α) resembles a normalization function, while Eα denotes Mittag-Leffler function. To broaden the ABC-derivative for the case of n < α < n + 1 having y(s)(a) = 0 f or s = 1, 2, ..., n, we have ABCDαy(t) = ABCDα(Dny(t)) = M(α) 1 −α ∫ χ 0 y(n+1)(s)Eα [ −α(t − s)α 1 −α ] d s, y(n+1)(s) = D(n+1)y(s) = Ddαey(s), (4) in which dαe represents ceil α. 2 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 3 Definition 2.1.3. The fractional integral of the ABC-derivative may be expressed by [30]: AB 0 I α t {y(t)} = 1 −α M(α) y(t)+ α M(α)Γ(α) ∫ t 0 y(s)(t−s)α−1d s.(5) Definition 2.1.4. The ABC-derivative having variable-order α(t), 0 < α(t) < 1 of function y(t) may be expressed by [31]: ABCDα(t)y(t) = M(α(t)) 1 −α(t) ∫ t 0 y ′ (s)Eα(t) [ −α(t) 1 −α(t) (t − s)α(t) ] d s, (6) where Eα(t) is the Mittag Leffer function. Theorem 2.1. Suppose β > 0. Then, the variable-order ABC-derivative may be expressed as [31]: ABCDα(t)tβ = M(α(t)) 1 −α(t) Γ(β+1)tβEα(t),β+1 [ −α(t) 1 −α(t) tα(t) ] .(7) 2.2. Some properties of SJPs The following recurrence formula can be used to create the famous Jacobi polynomials, which are specified in the interval [−1, 1]: J(u,v)i (t) = (u + v + 2i − 1)[u2 − v2 + t(u + v + 2i)(u + v + 2i − 2)] 2i(u + v + i)(u + v + 2i − 2) J(u,v)i−1 (t) − (u + i − 1)(v + i − 1)(u + v + 2i) i(u + v + 2i)(u + v + 2i − 2) J(u,v)i−2 , for i = 2, 3, ..., where J(u,v)0 (t) = 1 and J (u,v) 1 (t) = u + v + 2 2 t + u − v 2 . We now define the shifted Jacobi polynomials by in- stigating the change of variable t = 2tL − 1 to apply these polynomials with respect to the interval t ∈ [0, L]. Let the shifted Jacobi polynomials J(u,v)i ( 2t L − 1) be ex- pressed by J(u,v)L,i (t). Then, J (u,v) L,i (t) can be generated from: J(u,v)L,i (t) = (u + v + 2i − 1)[u2 − v2 + ( 2tL − 1)(u + v + 2i)(u + v + 2i − 2) 2i(u + v + i)(u + v + 2i − 2) J(u,v)L,i−1(t) − (u + i − 1)(v + i − 1)(u + v + 2i) i(u + v + i)(u + v + 2i − 2) J(u,v)L,i−2(t), i = 2, 3, ..., (8) where J(u,v)L,0 (t) = 1 and J (u,v) L,1 (t) = u+v+2 2 ( 2t L − 1) + u−v 2 . The analytic form with regards to the shifted Jacobi poly- nomials J(u,v)L,i (t) of degree i may be expressed as J(u,v)L,i (t) = i∑ k=0 (−1)i−k Γ(i + v + 1)Γ(i + k + u + v + 1) Γ(k + v + 1)Γ(i + u + v + 1)(i − k)!k!Lk tk, (9) where J(u,v)L,i (0) = (−1) i Γ(i+v+1) Γ(v+1)i! , J (u,v) L,i (L) = Γ(i+u+1) Γ(u+1)i! . Of these polynomials, the most commonly utilized are the shifted Chebyshev polynomials with respect to the first kind TL,i(t), the shifted Legendre polynomials PL,i(t), as well as the shifted Chebyshev polynomials with respect to the second kind UL,i(t). Moreover, for the non-symmetric shifted Jacobi polynomials, two essential special cases with regards to shifted Chebyshev polynomials of third and fourth kinds VL,i(t) and WL,i(t) are also taken into account. The following relations connect these orthogonal polynomials with respect to the shifted Jacobi polynomials. TL,i(t) = i!Γ(0.5) Γ(i+0.5) J (−0.5,−0.5) L,i (t), PL,i(t) = J (0,0) L,i (t), UL,i(t) = (i+1)!Γ(0.5) Γ(i+1.5) J 0.5,0.5 L,i (t), VL,i(t) = (2i)!! (2i−1)!! J (0.5,−0.5) L,i (t), WL,i(t) = (2i)!! (2i−1)!! J (−0.5,0.5) L,i (t). The orthogonality condition with respect to the shifted Jacobi polynomials is expressed as ∫ L 0 J(u,v)L, j (t)J (u,v) L,k (t)W (u,v) L (t)dt = hk, (10) where W (u,v)L (t) = t v(L − t)u and hk =  Lu+v+1 Γ(k+u+1)Γ(k+v+1) (2k+u+v+1)k!Γ(k+u+v+1) , i = j, 0, i , j. Let y(t) refers to a polynomial with degree n. Now, we may write these in terms of shifted Jacobi polynomials given by y(t) = N∑ j=0 c j J (u,v) L, j (t) = c T , (11) in which the coefficients c j are provided as follows c j = 1 h j ∫ L 0 W (u,v)L (t)y(t)J (u,v) L, j (t)dt, j = 0, 1, .... (12) Suppose the shifted Jacobi coefficient vector C, as well as the shifted Jacobi vector φ(t), is expressed as CT = [c0, c1, ..., cN ], (13) 3 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 4 and φ(t) = [J(u,v)L,0 (t), J (u,v) L,1 (t), ..., J (u,v) L,N (t)] T , (14) accordingly. Therefore, the first-order derivative with respect to the vector φ(t) may be written as dφ(t) dt = D(1)φ(t), (15) in which D(1) denotes the (N + 1) × (N + 1) operational matrix of derivative written as D(1) = (di j) = C1(i, j), i > j,0 otherwise, in which C1(i, j) = Lu + v(i + u + v + 1)(i + u + v + 2) j( j + u + 2)i− j−1Γ( j + u + v + 1) (i − j − 1)!Γ(2 j + u + v + 1) × 3F2  −i + 1 + j, i + j + u + v + 2, j + u + 1 ; 1 j + u + 2, 2 j + u + v + 2  (The proof can be found in [24], and the general definitions of a generalized hypergeometric series, as well as special 3F2, may be found in [25], accordingly on pp. 41 and 103–104). For instance, for even N, we obtain D(1) =  0 0 0 . . . 0 0 C1(1, 0) 0 0 . . . 0 0 C1(2, 0) C1(2, 1) 0 . . . 0 0 C1(3, 0) C1(3, 1) C1(3, 2) . . . 0 0 ... ... ... . . . ... ... C1(N, 0) C1(N, 1) C1(N, 2) . . . C1(N, N − 1) 0.  (16) 3. Generalized SJPs operational matrix to fractional calculus This section’s major goal is to expand the Jacobi operational matrix (JOM) of derivatives for Atangana-Baleanu Caputo (ABC) [32, 33]. Theorem 3.1. Suppose Ψ(t) vector be SJPs defined in Eq.(11) such that α > 0. Then ABCDαΨ(t) ' ABCD(α)Ψ(t), (17) in which Dα denotes the operational matrix (m + 1) × (m + 1) that may be expressed as: Let φ(t) vector be SJPs defined in Eq.(14). Here, suppose α > 0, then the ϑ in OM Eq.(18) is obtained using SJPs as follows ϑi, j,k = M(α) (1 −α) j∑ ł=0 (−1)i+ j−k+łΓ(i + v + 1)Γ(i + k + u + v + 1)Γ( j + v + 1)Γ( j + ł + u + v + 1) h jΓ(k + v + 1)Γ(i + u + v + 1)(i − k)!LkΓ( j + v + 1)Γ( j + u + v + 1)( j − ł)!ł!Lł a j,ł. (19) Proof. ABCDαtv = M(α) 1 −α ∫ t 0 Y (n+1)(s)Eα [ −α(t − s)α 1 −α ] d s, v > 1, v > dαe = M(α) 1 −α ∫ t 0 Γ(v + 1) Γ(v − n) S v−n−1 Eα [ −α(t − s)α 1 −α ] d s I f 0 < α < 1, bαc = n = 0, ψ(t) is solve f or ∫ t 0 S v−1 Eα [ −α(t − s)α 1 −α ] d s ABCDαJu,vL,i (t) = i∑ k=dαe (−1)i− jΓ(i + v + 1)Γ(i + k + u + v + 1) Γ(k + v + 1)Γ(i + u + v + 1)(i − k)!k!Lk ABCDαtk = i∑ k=dαe M(α) 1 −α ψ(t) (−1)i−kΓ(i + v + 1)Γ(i + k + u + v + 1) Γ(k + v + 1)Γ(i + u + v + 1)(i − k)!LkΓ(k)4 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 5 where ψ(t) = µ∑ j=0 ak, j J u,v L, j(t) ak, j = 1 h j j∑ ł=0 Γ( j + v + 1)Γ( j + ł + u + v + 1) Γ( j + v + 1)Γ( j + u + v + 1)( j − ł)!ł!Lł ∫ 1 0 ψ(t)W (u,v)L t łdt ABCDαJu,vL,i (t) = i∑ k=dαe µ∑ j=0 (−1)i−kΓ(i + v + 1)Γ(i + k + u + v + 1) Γ(k + v + 1)Γ(i + u + v + 1)(i − k)!k!Lk ak, j J u,v L, j(t) = m∑ j=0 ( i∑ k=dαe ϑi, j,k)J u,v L, j(t). Corollary 1. In the case of u = v = 0, it is apparent that the JOM of derivatives for integer calculus aligns with Legendre operational matrix of derivatives with respect to integer calculus as gained by Saadatmandi and Dehghan (refer [12] Eq. (11)). Corollary 2. In the case of u = v = −0.5, it is evident that the JOM of derivatives for integer calculus aligns with Chebyshev’s operational matrix of derivatives with respect to integer calculus as gained by Doha et al. (refer [15] Eq. (3.2)). 4. Applications of the operational matrix of fractional derivative This section solves an FDE in order to demonstrate the great significance of an operational matrix based on SLPs of frac- tional derivatives. We follow the same steps when using SJPs. 4.1. Linear FDEs Consider the linear FDEs ABCDαy(t) = b1 D βky(t) + ... + bk D β1y(t) + bk+1y(t) + bk+2q(t), for k = 1, 2, .... (20) The initial conditions are y(v)0 = dv, v = 0, ..., n, (21) in which bk denotes real constant coefficients with n < α ≤ n+1, 0 < β1 < β2 < ... < βk < α, in which ABCDβ refers to the ABC- derivative of order β. To solve Eq.(20), we present an approximation of the function q(t) ' m∑ i=0 qi J (u,v) L,i (t) = Q T Ψ(t), (22) ABCDαy(t) ' CT ABCD(α)Ψ(t), (23) where Q = [q0, ..., qm]T is a known vector. Employing Eqs.(22),(23) and (13), the residual Rm(t) for Eq.(20) can be written as Rm(t) ' ( CT ABCD(α) − bk+1C T − bk+2 Q T ) Ψ(t). (24) We now establish m − n linear equations as in a typical tau method by applying 〈Rm(t),p j(t)〉 = ∫ 1 0 Rm(t)pi(t)dt = 0, i = 0, 1, ..., m − n − 1. (25) Moreover, by substituting Eq.(13) and Eq.(14) with Eq.(21), we get y0 = C T Ψ(0) = d0, y(1)0 = C T D(1)Ψ(0) = d1, ... y(n)0 = C T D(n)Ψ(0) = dn. (26) Eqs. (25) and (26) generate (m−n) linear equations, which may be solved using arbitrary coefficients with respect to the vector C. 4.2. Nonlinear FDEs Consider the non-linear FDEs ABCDαy(t) = F(t, y(t), Dβ1 y(t), ...Dβk y(t)). (27) The initial conditions are y(v)0 = dv, v = 0, ..., n, (28) in which n < α ≤ n + 1, 0 < β1 < β2 < ... < βr < α, as well as ABCDα resembles the ABC-derivative of order α. We put in mind that F can be non-linear in general. CT ABCD(α)Ψ(t) ' F(t, CT Ψ(t), CT D(β1 )Ψ(t), ..., CT D(βk )Ψ(t)). (29) Furthermore, upon substituting Eq.(13) and Eq.(14) with Eq.(28), we now have Y0 = C T Ψ(0) = d0, Y (v)0 = C T D(v)Ψ(0) = dv, v = 1, 2, ..., n. (30) To discover the solution y(t), we collocate Eq.(29) by employ- ing first (m−n) points shifted Legendre roots of P̄m+1(t). These equations along with Eq.(30) establish (m + 1) non-linear equa- tions, which can be resolved by employing Newton’s iterative method. 5 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 6 5. JOM of variable-order ABC-derivative This section is devoted to tackling the problem using the SJPs OM of variable order. λ(t) = [1, t, t2, ..., tn]T . (31) Thus, the vector Ψ(t) can be expressed as: Ψ(t) = A(u,v)λ(t), (32) in which A(u,v) is (n + 1) × (n + 1) denotes a square matrix that specified by: (ai, j)0≤i, j≤n =  (−1)n−i Γ(n+β+1)Γ(n+i+α+β+1) Γ(i+β+1)Γ(n+α+β+1)Γ(n−i+1)Γ(i+1)li , i ≥ j, 0, otherwise. (33) Hence, by employing Eq.(32), we now have λ(t) = A−1Ψ(t). (34) Using the OM for variable-order fractional differential operator Dα(t)Ψ(t), and Eq.(32), we now have Dα(t)Ψ(t) = Dα(t)(Aλ(t)) = ADα(t)[1, t, t2, ..., tι]T . (35) Here, the Atangana-Baleanu Caputo (ABC) derivative with re- spect to the variable order provided in Eq.(4) may be employed. Then, we may obtain Eq.(35) as given below: Dα(t)Ψ(t) = [0, Γ(2) Γ(1 −α(t)) t ∞∑ k=0 ( −α(t)1−α(t) t α(t))k Γ(kα(t) + 2) , Γ(3) Γ(1 −α(t)) t2 ∞∑ k=0 ( −α(t)1−α(t) t α(t))k Γ(kα(t) + 3) , ..., Γ(ι + 1) Γ(1 −α(t)) tι ∞∑ k=0 ( −α(t)1−α(t) t α(t))k Γ(kα(t) + ι + 1) ]T , Dα(t)Ψ(t) = AB(t)λ(t), (36) in which B(t) =  0 0 0 . . . 0 0 Γ(2) Γ(1−α(t)) ∑ ∞ k=0 ( −α(t)1−α(t) t α(t) )k Γ(kα(t)+2) 0 . . . 0 0 0 Γ(3) Γ(1−α(t)) ∑ ∞ k=0 ( −α(t)1−α(t) t α(t) )k Γ(kα(t)+3) . . . 0 ... ... ... ... ... 0 0 0 . . . Γ(ι+1) Γ(1−α(t)) ∑ ∞ k=0 ( −α(t)1−α(t) t α(t) )k Γ(kα(t)+ι+1)  (37) Substituting Eq.(34) into Eq.(36), we get Dα(t)Ψ(t) = AB(t)A−1Ψ(t), (38) in which AB(t)A−1 denotes the OM of the variable-order ABC- derivative Dα(t)Ψ(t). Here, the approximate solution may be given as Dα(t)y(t) ' Dα(t)(CT Ψ(t)) = CT Dα(t)Ψ(t) = CT AB(t)A−1Ψ(t),(39) CT AB(t)A−1Ψ(t) = F[t, CT Ψ(t), CT AD(1) A−1Ψ(t), ..., CT AD(n) A−1Ψ(t)], 0 ≤ t ≤ 1. (40) Here, we employ the collocation points, tu = 2u+1 2n+2 , u = 0, 1, ..., n, in converting the system of equations given in Eq.(40) into an algebraic equations system as follows: CT AB(tu)A −1 Ψ(tu) = F[tu, C T Ψ(tu), C T AD(1) A−1Ψ(tu), ..., CT AD(n) A−1Ψ(tu)], CT Ψ(0) = y0 (41) Ultimately, the arbitrary vector C in Eq.(9) may be gained by solving the algebraic equations system provided in Eq.(41). 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 0 0.05 0.1 y( t) a approximate exact 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 0 0.05 y( t) b approximate exact 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 0 0.05 y( t) c approximate exact Figure 1: Comparison between exact and approximate solution for (a)α = 0.95, (b)α = 0.9 and (c)α = 0.85 for Example 6.2. Example 5.1. Suppose the following[36] ABCDαy(t) = y2(t) − 2(t + 1)−2, 6 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 7 α Method t = 0.1 t = 0.3 t = 0.5 t = 0.7 t = 0.9 0.85 LOM 2.19158e-3 2.89152e-3 4.26507e-3 4.40040e-4 3.13683e-3 PRCO 2.36581e-3 4.24947e-3 6.00249e-3 7.71664e-3 9.42102e-3 MTSLP 1.08027e-2 5.75464e-3 3.35770e-3 1.33299e-3 5.46216e-4 JOM(0,0.5) 3.00081e-3 5.85134e-4 2.97585e-3 2.07844e-3 1.57973e-4 JOM(0.5,0) 3.28839e-3 1.50328e-3 5.29734e-3 2.26751e-3 1.67970e-3 0.9 LOM 1.26158e-3 2.08766e-3 2.57667e-3 5.91395e-4 1.94142e-3 PRCO 1.18803e-3 2.39605e-3 3.56382e-3 4.72981e-3 5.90684e-3 MTSLP 1.10428e-2 7.41747e-3 5.88865e-3 5.57577e-3 3.34146e-3 JOM(0,0.5) 2.51233e-3 1.28073e-4 2.28507e-3 1.93179e-3 7.25019e-4 JOM(0.5,0) 2.01129e-3 1.21534e-3 3.35330e-3 1.02220e-3 1.41226e-3 0.95 LOM 4.88352e-4 1.16544e-3 1.10247e-3 5.09878e-4 8.67245e-4 PRCO 4.09639e-4 9.83473e-4 1.55920e-3 2.14596e-3 2.74713e-3 MTSLP 9.83430e-3 8.58128e-3 7.92322e-3 7.35794e-3 6.82360e-3 JOM(0,0.5) 1.77326e-3 6.45112e-4 1.54473e-3 1.53799e-3 1.01582e-3 JOM(0.5,0) 8.75527e-4 7.70414e-4 1.55386e-3 1.52607e-4 9.24438e-4 Table 1: The absolute error obtained by employing various values of α for Example 6.2. α Method t = 0.1 t = 0.3 t = 0.5 t = 0.7 t = 0.9 0.85 LOM 8.14287e-4 1.92419e-3 1.00572e-2 4.67247e-2 1.15441e-1 PRCO 1.02731e-3 7.40800e-3 1.19908e-2 5.84856e-2 1.50751e-1 MTSLP 4e-3 2.74688e-3 1.97787e-3 1.91297e-2 7.45086e-2 JOM(0,0.5) 8.01348e-4 1.93026e-3 1.00587e-2 4.67324e-2 1.15433e-1 JOM(0.5,0) 8.16107e-4 1.91932e-3 1.00499e-2 4.67267e-2 1.15445e-1 0.9 LOM 6.70490e-4 1.42673e-3 8.20714e-3 3.66835e-2 8.89397e-2 PRCO 8.64541e-4 1.70507e-3 9.89488e-3 4.55878e-2 1.14652e-1 MTSLP 4e-3 2.06924e-3 2.35421e-3 1.18382e-2 3.19227e-1 JOM(0,0.5) 6.62878e-4 1.42923e-3 8.20792e-3 3.66879e-2 8.89366e-2 JOM(0.5,0) 6.71528e-4 1.42406e-3 8.20355e-3 3.66839e-2 8.89426e-2 0.95 LOM 4.13897e-4 7.86685e-4 5.01340e-3 2.15783e-2 5.13325e-2 PRCO 5.90418e-4 8.64405e-4 6.50201e-3 2.77896e-2 2.90733e-1 MTSLP 4e-3 6.57765e-3 3.02391e-3 1.94353e-3 1.90389e-2 JOM(0,0.5) 4.11399e-4 7.86925e-4 5.01366e-3 2.15797e-2 5.13320e-2 JOM(0.5,0) 2.94860e-4 3.13253e-4 3.61277e-3 1.25339e-2 2.53848e-2 Table 2: The absolute error obtained employing various values of α for Example 6.3. For y0 = −2 and the exact solution y(t) = −2 (t+1) in case of α = t 0, Figure 3 displays the approximate values of α = 0.85, 0.9, 0.95 and m = 6. A good approximation that is comparable to the exact answer can be obtained via an operation matrix based on SJPs. 6. Numerical examples The numerical examples of linear and non-linear fractional- order and variable-order scenarios will acquire some solutions in this section. Our computational findings will measure the difference between the exact and approximate solutions using absolute error. The MATLAB R2020b software is used to code and perform all of the numerical programs, whereas the CPU is for the next. • JOM Jacobi Operational matrix method derived in this study. • CPSKOM Chebyshev polynomials with respect to the second kind Operational matrix method derived in this study. • LOM Legendre Operational matrix method [26]. • PRCO Predictor-Corrector method provided in [27]. • MTSLP Mixture two-step Lagrange polynomial as well 7 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 8 α Method t = 0.1 t = 0.3 t = 0.5 t = 0.7 t = 0.9 0.85 LOM -1.74986 -1.50024 -1.35680 -1.23342 -1.13618 PRCO -1.60319 -1.52292 -1.41361 -1.32635 -1.25359 MTSLP -1.67179 -1.45348 -1.31389 -1.22725 -1.17437 JOM(0,0.5) -1.76123 -1.50623 -1.35787 -1.23557 -1.13710 JOM(0.5,0) -1.74359 -1.49912 -1.35666 -1.23182 -1.13669 0.9 LOM -1.78420 -1.52494 -1.36023 -1.22657 -1.11789 PRCO -1.65250 -1.54279 -1.42052 -1.32416 -1.24500 MTSLP -1.78773 -1.58943 -1.45396 -1.35011 -1.26552 JOM(0,0.5) -1.79249 -1.53108 -1.36213 -1.22826 -1.11925 JOM(0.5,0) -1.77931 -1.52306 -1.36002 -1.22533 -1.11773 0.95 LOM -1.82307 -1.54990 -1.35819 -1.21279 -1.09254 PRCO -1.71332 -1.56366 -1.42509 -1.31793 -1.23183 MTSLP -1.89701 -1.64139 -1.47105 -1.34753 -1.25202 JOM(0,0.5) -1.82560 -1.55437 -1.36040 -1.21334 -1.09386 JOM(0.5,0) -1.82039 -1.54740 -1.35792 -1.21227 -1.09156 1 Exact -1.81818 -1.53846 -1.33333 -1.17647 -1.05263 Table 3: The approximate solutions obtained employing various values of α for Example 6.4. α Method t = 0.1 t = 0.3 t = 0.5 t = 0.7 t = 0.9 0.85 LOM 0.13534 0.20345 0.24998 0.28275 0.31040 PRCO 0.19782 0.22782 0.27041 0.30270 0.32864 MTSLP 0.164103 0.22650 0.25520 0.28130 0.30738 JOM(0,0.5) 0.12887 0.20271 0.24847 0.28210 0.30963 JOM(0.5,0) 0.13737 0.20354 0.25049 0.28286 0.31040 0.9 LOM 0.11432 0.19514 0.24890 0.28804 0.31966 PRCO 0.16976 0.21614 0.26616 0.30394 0.33403 MTSLP 0.10613 0.18237 0.23502 0.27718 0.31138 JOM(0,0.5) 0.10899 0.19386 0.24741 0.28717 0.31881 JOM(0.5,0) 0.11608 0.19533 0.24944 0.28817 0.31989 0.95 LOM 0.08994 0.18737 0.24973 0.29612 0.33192 PRCO 0.13780 0.20414 0.26292 0.30695 0.34158 MTSLP 0.051491 0.15966 0.23010 0.28087 0.32011 JOM(0,0.5) 0.08644 0.18564 0.24848 0.29510 0.33113 JOM(0.5,0) 0.09129 0.18767 0.25026 0.29628 0.33233 1 Exact 0.08375 0.19264 0.26323 0.31422 0.35351 Table 4: The approximate solutions obtained using different values of α for Example 6.5. as the fundamental theorem with respect to fractional cal- culus stated in [28]. Example 6.1. We now consider the Bagley–Torvik equation governing the motion of a rigid plate immersed in the Newto- nian fluid given as follows ABCD1.5y(t) + D2y(t) + y(t) = t + 1. Here, y0 = t0, y ′ 0 = t 0 and y(t) = t + 1 denotes the exact solution. Using SLPs, the approximate solution for m = 3 is y(t) = [1.5 0.5 0 0]Ψ(t) = t + 1, which equals to the exact solution. Using SJPs(0.5,0), the approximate solution for m = 3 is y(t) = [1.4 0.4 0 0]φ(t) = t + 1, which equals to the exact solution. For SJPs(0,0.5), the approximate solution for m = 3 is y(t) = [1.6 0.4 0 0]φ(t) = t + 1, which equals to the exact solution. 8 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 9 t Error Error Error Error Error Error LOM PRCO MTSLP CPSKOM JOM(0,0.5) JOM(0.5,0) 0.1 7.71674e-12 1.09802e-4 1.1e-2 2.09480e-9 8.56178e-11 2.50289e-11 0.2 7.77495e-12 4.53577e-4 2.24064e-2 3.26627e-9 8.53458e-11 3.76301e-11 0.3 8.81307e-12 1.09934e-3 3.46062e-2 3.66969e-9 8.97760e-11 3.88097e-11 0.4 2.77689e-11 2.12757e-3 4.60297e-2 3.46035e-9 2.89525e-10 3.28610e-11 0.5 3.48143e-11 3.61364e-3 5.77915e-2 2.79354e-9 3.63677e-10 2.40773e-11 0.6 1.56708e-11 5.63079e-3 6.74063e-2 1.82454e-9 1.62011e-10 1.67516e-11 0.7 4.39398e-11 8.25300e-3 7.82267e-2 7.08644e-10 4.65697e-10 1.51773e-11 0.8 1.58296e-10 1.15567e-2 8.50985e-2 3.98864e-10 1.66967e-9 2.36475e-11 0.9 3.41676e-10 1.56209e-2 9.75044e-2 1.34269e-9 3.60013e-9 4.64555e-11 Table 5: The absolute error for Example 6.6 for m = 4 t Error Error Error Error Error Erorr LOM PRCO MTSLP CPSKOM JOM(0,0.5) JOM(0.5,0) 0.1 1.75338e-4 6.54123e-4 e-2 1.75338e-4 1.75338e-4 1.75338e-4 0.2 3.38732e-3 3.95957e-3 2.70108e-2 3.38732e-3 3.38732e-3 3.38732e-3 0.3 9.96953e-3 9.35838e-3 3.72260e-2 9.96953e-3 9.96953e-3 9.96953e-3 0.4 1.88528e-2 1.27167e-2 4.22493e-2 1.88528e-2 1.88528e-2 1.88528e-2 0.5 2.93188e-2 8.50431e-3 4.70844e-2 2.93188e-2 2.93188e-2 2.93188e-2 0.6 4.06491e-2 8.72007e-3 5.81601e-2 4.06490e-2 4.06490e-2 4.06490e-2 0.7 5.21251e-2 4.36191e-2 8.11371e-2 5.21250e-2 5.21250e-2 5.21250e-2 0.8 6.30285e-2 1.00115e-1 1.20052e-1 6.30284e-2 6.30284e-2 6.30284e-2 0.9 7.26408e-2 1.82050e-1 1.77509e-1 7.26406e-2 7.26406e-2 7.26406e-2 Table 6: The absolute error for Example 6.7 for m = 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 y( t) =0.85 =0.9 =0.95 exact Figure 2: Comparison between approximate solutions for α = 0.85, α = 0.9 and α = 0.95 with the exact solution for Example 6.3. Example 6.2. Suppose we have the following model[34]: ABCDαy(t) = −K(1 − y(t)), 0 < α < 1, and y0 = 0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 y( t) =0.85 =0.9 =0.95 exact Figure 3: Comparison of approximate solutions for α = 0.85, α = 0.9 and α = 0.95 with the exact solution in case α = t0 for Example 6.4. The exact solution is as follows: y(t) = −K(1 −α) M(α) − K(1 −α) Eα ( Kα M(α) − K(1 −α) tα ) + [ 1 − Eα ( Kα M(α) − K(1 −α) tα )] + M(α)y(0) M(α) − K(1 −α) Eα ( Kα M(α) − K(1 −α) tα ) . 9 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 y( t) =0.85 =0.9 =0.95 exact Figure 4: Comparison between approximate solutions for α = 0.85, α = 0.9 α = 0.95 with exact solution in case of α = t0 for Example 6.5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 y( t) Exact LOM MTSLP PRCO Figure 5: The approximate solution and exact solution for Example 6.6 for m = 4 Table 1 and Figure 1 compare the absolute error for the two methods for K is constant K = 0.1,α = 0.85, 0.9, 0.95, and m = 4 respectively. An enhanced approximate solution comparable to the exact solution can be obtained using an operation matrix based on SJPs. Example 6.3. Suppose we have the following[35] ABCDαy(t) = −y(t) + t4 − 0.5t3 − 3 Γ(4 −α) t3−α + 24 Γ(5 −α) t4−α, for y0 = 0, while the exact solution is expressed by y(t) = t4 − 0.5t3. Table 2 and Figure 2 compare the absolute error for the two approaches for α = 0.85, 0.9, 0.95, and m = 6. A good approx- imate solution that is comparable to the exact answer can be obtained using an operation matrix based on SJPs. Example 6.4. Suppose we have the following[37] ABCDαy(t) = (1 − y(t))4, 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 y( t) Exact LOM MTSLP PRCO Figure 6: The approximate solution and exact solution for Example 6.7 for m = 4 The exact solution in case of α = t0 is expressed by y(t) = 1+3t−(1+6t+9t2 )1/3 1+3t and y0 = 0. The problem is solved with m = 8, and the numerical results are shown in Figure 4. Example 6.5. Consider the following linear variable-order FDEs[38]: ABCDα(t)y(t) + ety(t) = et(t2 + t3 + 1) + M(α(t)) 1 −α(t) 2t2 Eα(t),3(− α(t) 1 −α(t) tα(t)) + M(α(t)) 1 −α(t) 6t3 Eα(t),4(− α(t) 1 −α(t) tα(t)), where α(t) = 0.5t + 0.1, y0 = t0 and the exact solution is pro- vided by y(t) = t2 + t3 + 1. The absolute error for m = 4 is shown in Table 3. An enhanced approximate solution that is comparable to the exact answer can be obtained using an operation matrix based on SJPs. Example 6.6. We now consider the following non-linear variable-order FDEs given by[39]: ABCDα(t)y(t) + y2(t) = t2 + t4 + 2t2−sin(t) Γ(3 − sin(t)) . where α(t) = 0.5t + 0.6, y0 = 0 and the exact solution is ex- pressed by y(t) = t2. The absolute error for m = 4 is shown in Figure 5. Operation matrix relying on SJPs may give an enhanced approximate so- lution that is comparable with the exact solution. 7. Conclusion We came up with a general formulation for the fractional and variable order Jacobi operational matrix (JOM), which is utilized to approximate Atangana-Baleanu Caputo (ABC) derivatives in numerical solutions. The shifting Jacobi tau and 10 Basim et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1221 11 collocation approaches served as the foundation for our strat- egy. Since the ABC fractional derivative enables the inclusion of traditional initial conditions in the formulation of the issue, the fractional derivatives are described in the ABC sense. Here, the findings presented in the preceding section show how ac- curate these algorithms are. Additionally, only a few shifted Jacobi polynomials are required to produce a good outcome. We can recommend some future studies for this paper: solve system of fractional and variable order differential equations by use Jacobi and more orthogonal polynomials. References [1] K. S. Miller, & B. Ross, “An introduction to the fractional calculus and fractional differential equations”, Wiley, (1993). [2] K. Oldham, & J. Spanier, “The fractional calculus theory and applications of differentiation and integration to arbitrary order”, Elsevier, (1974). [3] M. Amairi, M. Aoun, S. Najar & M. N. Abdelkrim, “A constant enclo- sure method for validating existence and uniqueness of the solution of an initial value problem for a fractional differential equation”, Applied Mathematics and Computation 217 (2010) 2162. [4] J. Deng, & L. Ma, “Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations”, Applied Mathe- matics Letters 23 (2010) 676. [5] L. K. Alzaki, & H. K. Jassim, “Time-Fractional Differential Equations with an Approximate Solution”, Journal of the Nigerian Society of Phys- ical Sciences (2022) 818. [6] O. A. Uwaheren, A. F. Adebisi, & O. A. Taiwo, “Perturbed collocation method for solving singular multi-order fractional differential equations of Lane-Emden type”, Journal of the Nigerian Society of Physical Sci- ences (2020) 141. [7] K. M. Owolabi, & A. Atangana, “On the formulation of Adams-Bashforth scheme with Atangana-Baleanu-Caputo fractional derivative to model chaotic problems”, Chaos: An Interdisciplinary Journal of Nonlinear Sci- ence 29 (2019) 023111. [8] S. Arshad, I. Saleem, O. Defterli, Y. Tang, & D. Baleanu, “Simpson’s method for fractional differential equations with a non-singular kernel ap- plied to a chaotic tumor model”, Physica Scripta 96 (2021) 124019. [9] S. Qureshi, & A. Yusuf, “Modeling chickenpox disease with fractional derivatives: From caputo to atangana-baleanu.” Chaos, Solitons & Frac- tals 122 (2019) 111. [10] I. Ahmed, E. F. D. Goufo, A. Yusuf, P. Kumam, P. Chaipanya, & K. Nonlaopon, “An epidemic prediction from analysis of a combined HIV- COVID-19 co-infection model via ABC-fractional operator”, Alexandria Engineering Journal 60 (2021) 2979. [11] S. Djennadi, N. Shawagfeh, M. S. Osman, J. F. Gómez-Aguilar, & O. A. Arqub, “The Tikhonov regularization method for the inverse source problem of time fractional heat equation in the view of ABC-fractional technique”, Physica Scripta 96 (2021) 094006. [12] A. Saadatmandi, & M. Dehghan, “A new operational matrix for solving fractional-order differential equations”, Computers & mathematics with applications 59 (2010) 1326. [13] E. H. Doha, A. H. Bhrawy, & S. S. Ezz-Eldien, “Efficient Cheby- shev spectral methods for solving multi-term fractional orders differential equations”, Applied Mathematical Modelling 35 (2011) 5662. [14] A. H. Bhrawy, A. S. Alofi, & S. S. Ezz-Eldien, “A quadrature tau method for fractional differential equations with variable coefficients”, Applied Mathematics Letters 24 (2011) 2146. [15] E. H. Doha, A. H. Bhrawy, & S. S. Ezz-Eldien, “A Chebyshev spectral method based on operational matrix for initial and boundary value prob- lems of fractional order”, Computers & Mathematics with Applications 62 (2011) 2364. [16] F. Ghoreishi, & S. Yazdani, “An extension of the spectral Tau method for numerical solution of multi-order fractional differential equations with convergence analysis”, Computers & Mathematics with Applications 61 (2011) 30. [17] S. K. Vanani, & A. Aminataei, “Tau approximate solution of fractional partial differential equations”, Computers & Mathematics with Applica- tions 62 (2011) 1075. [18] A. Pedas, & E. Tamme, “On the convergence of spline collocation meth- ods for solving fractional differential equations”, Journal of Computa- tional and Applied Mathematics 235 (2011) 3502. [19] S. Esmaeili, & M. Shamsi, “A pseudo-spectral scheme for the approxi- mate solution of a family of fractional differential equations”, Communi- cations in Nonlinear Science and Numerical Simulation 16 (2011) 3646. [20] S. Esmaeili, M. Shamsi, & Y. Luchko, “Numerical solution of fractional differential equations with a collocation method based on Müntz polyno- mials”, Computers & Mathematics with Applications 62 (2011) 918. [21] G. Szegö, “Am. Math. Soc. Colloq. Pub. 23”, Orthogonal Polynomials (1985). [22] E. H. Doha, A. H. Bhrawy, & R. M. Hafez “A Jacobi dual-Petrov-Galerkin method for solving some odd-order ordinary differential equations”, Ab- stract and Applied Analysis 2011 (2011). [23] E. H. Doha, A. H. Bhrawy, & S. S. Ezz-Eldien, “A new Jacobi operational matrix: an application for solving fractional differential equations”, Ap- plied Mathematical Modelling 36 (2012) 4931. [24] E. H. Doha, “On the construction of recurrence relations for the expansion and connection coefficients in series of Jacobi polynomials”, Journal of Physics A: Mathematical and General 37 (2004) 657. [25] Y. L. Luke, “The special functions and their approximations”, 53 (1969). [26] M. Basim, N. Senu, Z. B. Ibrahim, A. Ahmadian, & S. Salahshour, “a Robust Operational Matrix of Nonsingular Derivative to Solve Fractional Variable-Order Differential Equations”, Fractals 30 (2022) 2240041. [27] J. D.Djida, A. Atangana, & I. Area, “Numerical computation of a frac- tional derivative with non-local and non-singular kernel”, Mathematical Modelling of Natural Phenomena 12 (2017) 4. [28] M. Toufik, & A. Atangana, “New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models”, The European Physical Journal Plus 132 (2017) 1. [29] M.Caputo, & M. Fabrizio, “A new definition of fractional derivative with- out singular kernel”, Progress in Fractional Differentiation & Applica- tions 1 (2015) 73. [30] A. Atangana, & D. Baleanu, “New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model”, arXiv preprint arXiv:1602.03408 (2016). [31] X. Li, Y. Gao, & B. Wu, “Approximate solutions of Atangana-Baleanu variable order fractional problems”, AIMS Mathematics 5 (2020) 2285. [32] L. J. Rong, & P. Chang, “Jacobi wavelet operational matrix of fractional integration for solving fractional integro-differential equation”, Journal of Physics: Conference Series 693 (2016) 012002. [33] C. Phang, Y. T. Toh, & F. S. Md Nasrudin, “An operational matrix method based on poly-Bernoulli polynomials for solving fractional delay differ- ential equations”, Computation 8 (2020) 82. [34] S. Salahshour, A. Ahmadian, M. Salimi, M. Ferrara, & D. Baleanu, “Asymptotic solutions of fractional interval differential equations with nonsingular kernel derivative”, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (2019) 083110. [35] A. Al-Rabtah, S. Momani, & M. A. Ramadan, “Solving linear and non- linear fractional differential equations using spline functions”, Abstract and Applied Analysis 2012 (2012). [36] Z. M. Odibat, & S. Momani, “An algorithm for the numerical solution of differential equations of fractional order”, Journal of Applied Mathemat- ics & Informatics 26 (2008) 15. [37] A. Al-Rabtah, S. Momani, & M. A. Ramadan, “Solving linear and non- linear fractional differential equations using spline functions”, Abstract and Applied Analysis 2012 (2012). [38] X. Li, Y. Gao, & B. Wu, “Approximate solutions of Atangana-Baleanu variable order fractional problems”, AIMS Mathematics 5 (2020) 2285. [39] D. Baleanu, “Approximate solutions for solving nonlinear variable-order fractional Riccati differential equations”, Inst Mathematics & Informatics (2019). 11