Novel derivative operational matrix in Caputo sense with applications

Abstract

The main objective of this study is to present a computationally efficient numerical method for solving fractional-order differential equations with initial conditions. The proposed method is based on the newly developed generalized derivative operational matrix and generalized integral operational matrix derived from Laguerre polynomials, which belong to the class of orthogonal polynomials. Through the utilization of these operational matrices, the fractional-order problems can be transformed into a system of Sylvester-type matrix equations. This system is easily solvable using any computational software, thereby providing a practical framework for solving such equations. The results obtained are compared against various benchmarks, including an existing exact solution, Podlubny numerical techniques, analytical and numerical solvers, and reported solutions from stochastic techniques employing hybrid approaches. This comparative analysis serves to validate the accuracy of our proposed design scheme.

Keywords:

• Laguerre polynomials
• operational matrices
• orthogonal polynomials
• spectral methods
• Tau method
• fractional derivative differential equations

Mathematic Subject classifications:

• 42C10
• 65D15
• 65D25
• 34G10

1. Introduction

Fractional Calculus (FC) is a generalization of integer-order calculus used for computing integrals and derivatives of non-integer orders. Initially, it was regarded solely as an abstract mathematical idea with minimal or no practical significance. Nevertheless, the present scenario regarding FC is entirely distinct, and commonly acknowledged as one of the most valuable and crucial tools in the scientific community. Various physical and biological phenomena find precise mathematical representation in FC, such as fractional-order models for nonlocal epidemics and partial bed-load transport. These models effectively capture the combined influences of sub-diffusion and non-local behaviour. In the proposed study [1], a fractional-order model was introduced for nonlocal epidemics, with a focus on examining the stability of equations with fractional orders. The implications of these findings were anticipated to extend to diseases like SARS, avian flu, and foot-and-mouth disease. The authors in [2] employed a fractional derivative model to precisely depict the combined effects of sub-diffusion and non-local behaviours in partial bed-load transport, induced by turbulence and micro-relief within an armor layer. Another insightful contribution is presented in [3], where a comprehensive discussion on variable-order differential operators in modelling anomalous diffusion was provided. Unlike constant-order fractional diffusion models, these new models exhibit variable properties dependent on concentration, time, space, and other independent variables.

We are always interested in the development of new numerical methods that are easy to apply and capable of producing efficient numerical results for both linear and non-linear ordinary and partial fractional-order differential equations (FODEs). Spectral methods and Finite-Difference methods are the two leading numerical tools that have been used widely to solve linear and non-linear ordinary and partial FODEs. Among the other finite-difference techniques, probably the predictor-corrector technique of Adams–Bashforth–Moulton type has been widely employed to solve FODEs, see [4,5]. However, Adams–Bashforth–Moulton type technique has two main disadvantages: one is their linear convergence; the other is their development for solving only commensurate-order FODEs. On the other hand, the spectral methods have been applied successfully for solving both commensurate-order and incommensurate-order FODEs, see [6,7]. The construction of the spectral methods involved the derivative and/or integral operational matrices associated with orthogonal polynomials, see [7–11] and references therein. Using these operational matrices, FODEs are transformed into a system of readily solvable matrix equations, which can then be solved by using any computational software. In addition, the spectral methods provides the exponential convergence as compared to the finite difference methods where the convergence is linear, i.e. as more basis functions are incorporated, the error diminishes swiftly, resulting in remarkably precise outcomes when compared to alternative numerical techniques.

In our current study, we explore the practical use of derivative and integral operational matrices associated with Laguerre polynomials (LPs) in the development of the proposed numerical method (PNM), named as The Operational Matrices Technique. This method falls within the category of spectral methods as it also provides the solution as series expansion of orthogonal polynomials. The structure of the PNM relies on the newly constructed fractional-order derivative operational matrix (FODOM) of LPs together its fractional-order integral operational matrix (FOIOM) derived in [12]. Through the application of FODOM and FOIOM, we are able to transform FODEs into a system of matrix equations of Sylvester type, which can be directly solved within the MATLAB environment using its built-in “lyap()” function.

Traditionally, the operational matrices approach, along with the Tau method, has been widely accepted in scientific literature for addressing FODEs, see [12–17]. However, the use of the operational matrices approach within the framework of spectral Tau methods introduces challenges related to the computation of Residual Functions. This limitation creates difficulties for its use, potentially leading to computational and precision issues. For instance, in [12], the operational matrices approach, combined with the Tau method, is suggested for solving FODEs. Nevertheless, it is noted as a computationally expensive method, as demonstrated by the data presented in Tables 3 and 5. The results shown in these tables highlight the computational expense associated with the method outlined in [12]. Furthermore, this approach is found to be less suitable for addressing problems with nonpolynomial solutions, as illustrated in Figure 6 of Example 7.5.

Similarly, the operational matrices approach, along with the Tau method, is employed in [15–17] for FODE solutions. However, to achieve promising results, the construction of Residual functions requires the consideration of additional terms of orthogonal polynomials, including Legendre and Laguerre. Consequently, the increased number of series coefficients evaluations makes this approach more costly and less user-friendly. The findings related to these results are detailed in Tables 4 and 9.

In contrast, our PNM's framework exclusively relies on the operational matrices of LPs, eliminating the need for constructing Residual Functions as required by methodologies introduced in [12–17]. Unlike these approaches, our proposed methodology can transform FODEs into a system of algebraic equations without the necessity of computing Residual Functions, as is the case with the hybrid approach based on the operational matrices together with Tau method. This distinctive feature enables us to achieve higher accuracy in the approximate results compared to the numerical techniques that combine the operational matrices approach with the Tau method or other reported numerical solvers in the literature. To provide a clearer understanding and distinguish our proposed methodology from other operational matrices-based approaches, we outline a step-by-step procedure for solving Example 7.1. Furthermore, in order to showcase the numerical precision of our PNM, we solve multiple examples and compare the obtained results with the available exact solution and the results reported in [12,15–17], and also with other numerical solvers including, Podlubny [18], Genetic Algorithm with Pattern Search technique (GA-PS) [19], the Particle Swarm Optimisation method combined with the Pattern Search technique (PSO-PS) [20], the Neural Network Interior Point Algorithm [21], Legendre artificial neural network method (LANNM) [22], and He's variational iteration method (He'VIM) [23]

The structure of the paper is outlined as follows: In Section 2, fundamental definitions and some useful properties of fractional-order operators are presented which are used throughout the paper. In Section 3, the analytical expressions of LPs and their properties are recalled. In Section 4, some useful results are derived to construct the operational matrices of LPs. The Section 4.1 discussed the derivation of newly proposed FODOM of LPs in Caputo sense. In Section 5, the construction of an improved and efficient numerical method is discussed. In Section 7, the applicability of the PNM is assessed by solving various examples and comparing the obtained results with the results obtained by employing other numerical solvers. The final section provides a summary of the findings and outlines potential research directions for future studies in this field.

2. Basic definitions and some useful properties

In this section, some useful properties and basic definitions of fractional-order operators are listed that will be used throughout this article. Multiple definitions are presented in the literature to define the fractional-order integrals and fractional-order derivatives operators, however, they do not coincide with each other, see [24]. The most commonly used definitions presented in the sense of Riemann–Liouville and Caputo are expressed in the following way [18,24].

Definition 2.1

The Riemann–Liouville fractional-order integral extends the notion of integration to include integrals of functions with non-integer orders. It's defined as: (1) $\begin{array}{rl}{}_{\mathrm{RL}}{J}_{b^{+}}^{\alpha }f\left(x\right)& =\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_b^x(x-y{)}^{\alpha -1}f(y)\mathrm{d}y,\\ & x>b,\alpha >0.\end{array}$(1)

Consequently, the fractional-order derivatives operators in Riemann–Liouville's sense and Caputo's sense are defined in the following way: (2) $\begin{array}{rl}{}_{\mathrm{RL}}{\mathcal{D}}_{b^{+}}^{\alpha }f\left(x\right)& =D_{\mathrm{RL}}^n{J}_{b^{+}}^{n-\alpha }f\left(x\right)\\ & =\frac{1}{\mathrm{\Gamma }(n-\alpha )}\frac{{\mathrm{d}}^{\mathrm{n}}}{\mathrm{d}{x}^n}{\int }_b^x(x-y{)}^{n-\alpha -1}\\ & \times f\left(y\right)\mathrm{d}y,x>b,\\ {}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }f\left(x\right)& {=}_{\mathrm{RL}}{J}_{b^{+}}^{n-\alpha }{D}^{n}f\left(x\right)\\ & =\frac{1}{\mathrm{\Gamma }(n-\alpha )}{\int }_b^x(x-y{)}^{n-\alpha -1}\\ & \times f^{\left(n\right)}\left(y\right)\mathrm{d}y,x>b,\end{array}$(2) where $n-1<\alpha <n,n\in \mathbb{N}$, and $\alpha >0$. The following are some useful properties that will be used later in establishing the main results. (3) $\begin{array}{rl}{}_{\mathrm{RL}}{J}_{b^{+}}^{\alpha }{(x-b)}^l& =\frac{\mathrm{\Gamma }(l+1)}{\mathrm{\Gamma }(l+1+\alpha )}\\ & \times {(x-b)}^{l+\alpha },\alpha \in {\mathbb{R}}_{\mathbb{+}},\end{array}$(3) (4) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }{(x-b)}^l& =\frac{\mathrm{\Gamma }(l+1)}{\mathrm{\Gamma }(l+1-\alpha )}{(x-b)}^{l-\alpha },\\ & l\in {\mathbb{R}}_{\mathbb{+}},\mathrm{\&}l\ge \lceil \alpha \rceil ,\end{array}$(4) (5) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }\left(f\right(x)+g(x\left)\right)& {=}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }f\left(x\right){+}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }g\left(x\right).\end{array}$(5)

3. LPs and its properties

LPs, introduced by the French mathematician Edmond Laguerre in the late 19th century, hold great significance across various disciplines, including mathematics, physics, and engineering. They form a set that is orthogonal over a specific interval, typically ranging from zero to infinity. This orthogonality property plays a pivotal role in numerous mathematical and computational techniques, including numerical integration, approximation theory, and solving differential equations. They serve as a comprehensive basis for expanding functions and addressing differential equations involving exponential decay, see [25–28]. Furthermore, LPs establish connections with probability theory, particularly in the analysis of random variables and stochastic processes. They find application in studying queueing systems and other probabilistic models, see [29–31]. Additionally, LPs have significant implications in quantum mechanics. They appear in the radial component of the solution to the Schrödinger equation for a one-electron atom. Additionally, they describe the static Wigner functions of oscillator systems in quantum mechanics when considering phase space. For instance, by employing the technique of separation of variables in spherical coordinates, the Schrödinger equation for the hydrogen-like atom can be solved exactly. In this solution, the radial component of the wave function takes the form of a (generalized) LPs, see [32–34].

The LPs has found extensive application in addressing significant real-world problems. For example, in the proposed study [35], the Lane-Emden equation, Bratu's type equation, and the nonlinear space-time Burger's type equation were tackled. The authors accomplished this by representing the solution function as a truncated expansion of a newly introduced generalization of LPs and applying the collocation spectral method. In another proposed work [36], the Euler-Bernoulli beam equation was addressed using a computationally efficient numerical algorithm based on Legendre and Laguerre polynomials.

The LPs are defined through a recursive formulation as [36,37] (6) $\begin{array}{rl}{L}_{i+1}\left(t\right)& =\frac{(2i+1-t)L_i\left(t\right)-\left(i\right)L_{i-1}\left(t\right)}{i+1},\\ L_0\left(t\right)& =1,L_1\left(t\right)=1-t.\end{array}$(6) The analytical expression of Equation (6(6) $\begin{array}{rl}{L}_{i+1}\left(t\right)& =\frac{(2i+1-t)L_i\left(t\right)-\left(i\right)L_{i-1}\left(t\right)}{i+1},\\ L_0\left(t\right)& =1,L_1\left(t\right)=1-t.\end{array}$(6) ) is given as (7) $L_i\left(t\right)=\sum _{k=0}^i\frac{(-1{)}^k}{k!}(\genfrac{}{}{0ex}{}{i}{k})t^k,i=0,1,2,\dots ,n.$(7) With respect to the weight function $w\left(t\right)$=${\mathrm{e}}^{-t}$, the set of LPs $\left\{L_0\right(t),L_1(t),\dots ,L_n(t\left)\right\}$ on the interval $[0,\mathrm{\infty })$ possesses an orthonormal basis that exhibits the following property (8) ${\int }_0^{\mathrm{\infty }}{L}_i\left(t\right)L_j\left(t\right)w\left(t\right)\mathrm{d}t={\delta }_{\mathrm{ij}},i,j=0,1,\dots ,n,$(8) where ${\delta }_{\mathrm{ij}}$ is the kronecker delta function. By employing Equation (7(7) $L_i\left(t\right)=\sum _{k=0}^i\frac{(-1{)}^k}{k!}(\genfrac{}{}{0ex}{}{i}{k})t^k,i=0,1,2,\dots ,n.$(7) ), we may compute the first few LPs for $i=0,1,\dots ,4$, given below (9) $\left(\begin{array}{c}{L}_0\left(t\right)\\ L_1\left(t\right)\\ L_2\left(t\right)\\ L_3\left(t\right)\\ L_4\left(t\right)\end{array}\right)=\left(\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ -{\displaystyle \frac{t^3}{6}+{\displaystyle \frac{3t^2}{2}-3t+1}}\\ {\displaystyle \frac{t^4}{24}-{\displaystyle \frac{2t^3}{3}+3t^2-4t+1}}\end{array}\right).$(9)

3.1. Function approximation using Laguerre series

Arbitrary function $u\in L^2\left(\right[0,\mathrm{\infty });{\mathrm{e}}^{-t})$, for which the integral (10) $||u|{|}^2={\int }_0^{\mathrm{\infty }}{\mathrm{e}}^{-t}|u\left(t\right){|}^2\mathrm{d}t<\mathrm{\infty }$(10) is finite, can be expanded into the Laguerre series as (11) $u\left(t\right)=\underset{n\to \mathrm{\infty }}{lim}\sum _{i=0}^n{c}_i{L}_i\left(t\right).$(11) The series coefficients, $c_i$ can be computed by using the following expression (12) $c_i={\int }_0^{\mathrm{\infty }}u\left(t\right){\mathrm{e}}^{-t}{L}_i\left(t\right)\mathrm{d}t.$(12) By truncating the series up to the first N-terms, we can express Equation (11(11) $u\left(t\right)=\underset{n\to \mathrm{\infty }}{lim}\sum _{i=0}^n{c}_i{L}_i\left(t\right).$(11) ) as follows (13) $u\left(t\right)\approx \sum _{i=0}^N{c}_i{L}_i\left(t\right)=F^T\psi \left(t\right),$(13) where, $\begin{array}{rl}{F}^T& =[c_0,c_1,c_2,\dots ,c_N],\\ \psi \left(t\right)& =\left[L_0\right(t),L_1(t),\dots ,L_N(t){]}^T.\end{array}$Consider the cubic function, $u\left(t\right)=4t^3-1$ that has a finite square norm with weight ${\mathrm{e}}^{-t}$, i.e. $||4t^3-1|{|}^2={\int }_0^{\mathrm{\infty }}(4t^3-1{)}^2{\mathrm{e}}^{-t}\mathrm{d}t=11473.$Therefore, $4t^3-1\in L^2\left(\right[0,\mathrm{\infty });{\mathrm{e}}^{-t})$ and it can be expanded into convergent Laguerre series (which is actually a finite sum) as follows $4t^3-1=\underset{n\to \mathrm{\infty }}{lim}\sum _{i=0}^n{c}_i{L}_i\left(t\right),$where $\begin{array}{rl}{c}_0& ={\int }_0^{\mathrm{\infty }}(4t^3-1){\mathrm{e}}^{-t}\mathrm{d}t=23,\\ c_1& =-72,c_2=72,c_3=-24,\end{array}$and all other coefficients are zeroes. Consequently, we get the following identity $4t^3-1=23L_0\left(t\right)-72L_1\left(t\right)+72L_2\left(t\right)-24L_3\left(t\right).$Now we consider the rational function, $u\left(t\right)=\frac{4t^3-1}{t^2+1}$ that has a finite square norm with weight ${\mathrm{e}}^{-t}$, i.e. ${\left|\left|\frac{4t^3-1}{t^2+1}\right|\right|}^2={\int }_0^{\mathrm{\infty }}{\left(\frac{4t^3-1}{t^2+1}\right)}^2{\mathrm{e}}^{-t}\mathrm{d}t\approx \mathrm{21.3606.}$Therefore, $\frac{4t^3-1}{t^2+1}\in L^2\left(\right[0,\mathrm{\infty });{\mathrm{e}}^{-t})$ and it can be expanded into convergent Laguerre series as follows $\frac{4t^3-1}{t^2+1}=\underset{n\to \mathrm{\infty }}{lim}\sum _{i=0}^n{c}_i{L}_i\left(t\right),$where $\begin{array}{rl}{c}_0& ={\int }_0^{\mathrm{\infty }}\frac{4t^3-1}{t^2+1}{\mathrm{e}}^{-t}{L}_0\left(t\right)\mathrm{d}t\approx 22.00504,\\ c_1& \approx -4.13738,c_2\approx 0.217678,c_3\approx 0.260623,\end{array}$and so on. We build Laguerre approximation of $\frac{4t^3-1}{t^2+1}$ across various values of N, see Figure 1.

Figure 1. Approximation of a rational function $u\left(t\right)=\frac{4t^3-1}{t^2+1}$ with various orthogonal polynomials on interval $[0,20]$. (a) Laguerre approximation at N = 10. (b) Laguerre approximation at N = 15. (c) Legendre approximation at N = 10. (d) Legendre approximation at N = 15. (e) Chelyshkov approximation at N = 10. (f) Chelyshkov approximation at N = 15. (g) Chebyshev approximation at N = 10 and (h) Chebyshev approximation at N = 15.

In a similar way, we build the Laguerre approximation for the following shifted Heaviside function, $H(t-a)$ for a positive number a, (see Figure 2). $\begin{array}{r}H(t-a)=\{\begin{array}{ll}1& \mathrm{for}t>a,\\ {\displaystyle \frac{1}{2}}& \mathrm{for}t=a,\\ 0& \mathrm{for}t<a.\end{array}\end{array}$

4. Operational matrices of LPs

Lemma 4.1

The Caputo fractional-order derivative of order α of LPs defined in (7(7) $L_i\left(t\right)=\sum _{k=0}^i\frac{(-1{)}^k}{k!}(\genfrac{}{}{0ex}{}{i}{k})t^k,i=0,1,2,\dots ,n.$(7) ) can be computed as (14) $\begin{array}{r}{}_{\mathrm{C}}{\mathcal{D}}_{0^{+}}^{\alpha }{L}_i\left(t\right)=\{\begin{array}{l}{\displaystyle \sum _{k=\lceil \alpha \rceil }^i\frac{(-1{)}^k(i!\left)\mathrm{\Gamma }\right(k+1)}{(k!{)}^2(i-k)!\mathrm{\Gamma }(k-\alpha +1)}{t}^{k-\alpha },}\\ i=\lceil \alpha \rceil ,\dots ,n,\\ 0,i=0,1,\dots ,\lceil \alpha \rceil -1.\end{array}\end{array}$(14)

Figure 2. Laguerre approximation of shifted Heaviside function for a = 1. (a) N = 9. (b) N = 14. (c) N = 19 and (d) N = 29.

Proof.

Using Equations (4(4) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }{(x-b)}^l& =\frac{\mathrm{\Gamma }(l+1)}{\mathrm{\Gamma }(l+1-\alpha )}{(x-b)}^{l-\alpha },\\ & l\in {\mathbb{R}}_{\mathbb{+}},\mathrm{\&}l\ge \lceil \alpha \rceil ,\end{array}$(4) ), (5(5) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }\left(f\right(x)+g(x\left)\right)& {=}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }f\left(x\right){+}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }g\left(x\right).\end{array}$(5) ), and (7(7) $L_i\left(t\right)=\sum _{k=0}^i\frac{(-1{)}^k}{k!}(\genfrac{}{}{0ex}{}{i}{k})t^k,i=0,1,2,\dots ,n.$(7) ), the result can be proved.

Lemma 4.2

For $t^{k-\alpha }\in L^2\left(\right[0,\mathrm{\infty });{\mathrm{e}}^{-t})$, we have the following result (15) $\begin{array}{rl}& t^{k-\alpha }\approx \sum _{j=0}^M{b}_j{L}_j\left(t\right),k-\alpha \in {\mathbb{R}}_{+},\end{array}$(15) (16) $\begin{array}{rl}& b_j=\sum _{k=l}^j\frac{(-1{)}^l(j!\left)\mathrm{\Gamma }\right(k-\alpha +l+1)}{(l!{)}^2(j-l)!},j=0,1,\dots ,n.\end{array}$(16)

Proof.

Since $t^{k-\alpha }\in L^2\left(\right[0,\mathrm{\infty });{\mathrm{e}}^{-t})$ then it can be expanded into Laguerre series as (17) $t^{k-\alpha }\approx \sum _{j=0}^n{b}_j{L}_j\left(t\right).$(17) By employing Equation (12(12) $c_i={\int }_0^{\mathrm{\infty }}u\left(t\right){\mathrm{e}}^{-t}{L}_i\left(t\right)\mathrm{d}t.$(12) ), the Laguerre series coefficients $b_j$ can be computed easily, given below (18) $b_j=\sum _{l=0}^j\frac{(-1{)}^{l}j!}{{l!}^2(j-l)!}\mathrm{\Gamma }(k-\alpha +l+1).$(18) Equations (17(17) $t^{k-\alpha }\approx \sum _{j=0}^n{b}_j{L}_j\left(t\right).$(17) ) and (18(18) $b_j=\sum _{l=0}^j\frac{(-1{)}^{l}j!}{{l!}^2(j-l)!}\mathrm{\Gamma }(k-\alpha +l+1).$(18) ) prove the result.

4.1. New generalized derivative operational matrix of LPs in Caputo sense

The operational matrices associated with orthogonal polynomials play a vital role in various domains of applied mathematics, particularly in the realms of computational mathematics. By employing operational matrices associated with orthogonal polynomials, it becomes possible to represent differential and integral operators as matrix operations. As a result, FODEs can be transformed into algebraic equations, simplifying the process of finding numerical solutions significantly.

In this section, we discuss the construction of a novel operational matrix of LPs in the sense of a Caputo fractional-order derivative operator.

Theorem 4.1

Let $\psi \left(t\right)=\left[L_0\right(t),L_1(t),\dots ,L_n(t){]}^T$ defined in Equation (13(13) $u\left(t\right)\approx \sum _{i=0}^N{c}_i{L}_i\left(t\right)=F^T\psi \left(t\right),$(13) ) then its FOIOM in the sense of Riemann–Liouville has dimensions $(N+1\times N+1)$ can be computed using the following result (19) ${}_{\mathrm{RL}}{J}^{\alpha }\psi \left(t\right)\approx Q_{(N+1\times N+1)}^{\left(\alpha \right)}\psi \left(t\right).$(19) The matrix $Q_{(N+1\times N+1)}^{\left(\alpha \right)}$, expresses the FOIOM of order $\alpha \in {\mathbb{R}}_{+}$ in the sense of Riemann-Liouville fractional-order integral operator, defined as follows: (20) $Q_{(N+1\times N+1)}^{\left(\alpha \right)}=\sum _{k=0}^i{\theta }_{(i,j)},$(20) where (21) $\begin{array}{rl}{\theta }_{(i,j)}& =\sum _{k=0}^i\sum _{r=0}^j\frac{(-1{)}^{k+r}(i!)(j!)\mathrm{\Gamma }(k+\alpha +r+1)}{{(r!)}^2(j-r)!k!(i-k)!\mathrm{\Gamma }(k+\alpha +1)},\\ & i,j=0,1,\dots ,N.\end{array}$(21)

Proof.

For the proof, see [12].

Theorem 4.2

Let $\psi \left(t\right)=\left[L_0\right(t),L_1(t),\dots ,L_n(t){]}^T$ defined in Equation (13(13) $u\left(t\right)\approx \sum _{i=0}^N{c}_i{L}_i\left(t\right)=F^T\psi \left(t\right),$(13) ) then its FODOM in the sense of Caputo has dimensions $(N+1\times N+1)$ can be computed using the following result (22) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }\psi \left(t\right)\approx R_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right).$(22) The matrix $R_{(N+1,N+1)}^{\left(\alpha \right)}$ is used to compute the FODOM of order $\alpha \in {\mathbb{R}}_{+}$ in the sense of Caputo fractional-order derivative operator, defined as follows: (23) $R_{(N+1,N+1)}^{\left(\alpha \right)}=\sum _{k=\lceil \alpha \rceil }^i{\phi }_{(i,j)},$(23) where (24) $\begin{array}{rl}\phi (i,j,k)& =\sum _{l=0}^j\frac{(-1{)}^{k+l}(i!\left)\right(j!\left)\mathrm{\Gamma }\right(k-\alpha +l+1)}{(l!{)}^2(k!\left)\right(j-l)!(i-k)!\mathrm{\Gamma }(k-\alpha +1)},\\ & i=\lceil \alpha \rceil ,\dots ,N,j=0,1,\dots ,N.\end{array}$(24)

Proof.

By employing Equation (5(5) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }\left(f\right(x)+g(x\left)\right)& {=}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }f\left(x\right){+}_{\mathrm{C}}{\mathcal{D}}_{b^{+}}^{\alpha }g\left(x\right).\end{array}$(5) ) and Lemma 4.1 to Equation (7(7) $L_i\left(t\right)=\sum _{k=0}^i\frac{(-1{)}^k}{k!}(\genfrac{}{}{0ex}{}{i}{k})t^k,i=0,1,2,\dots ,n.$(7) ), we have (25) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)=\sum _{k=\lceil \alpha \rceil }^i\frac{(-1{)}^{k}i!\mathrm{\Gamma }(k+1)}{{k!}^2(i-k)!\mathrm{\Gamma }(k-\alpha +1)}{t}^{k-\alpha }.$(25) Using Lemma 4.2, Equation (25(25) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)=\sum _{k=\lceil \alpha \rceil }^i\frac{(-1{)}^{k}i!\mathrm{\Gamma }(k+1)}{{k!}^2(i-k)!\mathrm{\Gamma }(k-\alpha +1)}{t}^{k-\alpha }.$(25) ) can be written as (26) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)& \simeq \sum _{k=\lceil \alpha \rceil }^i\frac{(-1{)}^{k}i!\mathrm{\Gamma }(k+1)}{{k!}^2(i-k)!\mathrm{\Gamma }(k-\alpha +1)}\\ & \times \sum _{j=0}^N\left(\sum _{l=0}^j\frac{(-1{)}^{l}j!}{{l!}^2(j-l)!}\mathrm{\Gamma }\right(k-\alpha +l+1\left)\right)\\ & \times L_j\left(t\right).\end{array}$(26) Equation (26(26) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)& \simeq \sum _{k=\lceil \alpha \rceil }^i\frac{(-1{)}^{k}i!\mathrm{\Gamma }(k+1)}{{k!}^2(i-k)!\mathrm{\Gamma }(k-\alpha +1)}\\ & \times \sum _{j=0}^N\left(\sum _{l=0}^j\frac{(-1{)}^{l}j!}{{l!}^2(j-l)!}\mathrm{\Gamma }\right(k-\alpha +l+1\left)\right)\\ & \times L_j\left(t\right).\end{array}$(26) ) can further be written as (27) $\begin{array}{r}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)\simeq \sum _{j=0}^N\left(\sum _{k=\lceil \alpha \rceil }^i{\phi }_{(i,j,k)}\right)L_j\left(t\right),i=\lceil \alpha \rceil ,\dots ,N.\end{array}$(27) The expression for $\phi (i,j,k)$ as provided in Equation (24(24) $\begin{array}{rl}\phi (i,j,k)& =\sum _{l=0}^j\frac{(-1{)}^{k+l}(i!\left)\right(j!\left)\mathrm{\Gamma }\right(k-\alpha +l+1)}{(l!{)}^2(k!\left)\right(j-l)!(i-k)!\mathrm{\Gamma }(k-\alpha +1)},\\ & i=\lceil \alpha \rceil ,\dots ,N,j=0,1,\dots ,N.\end{array}$(24) ), is used to rephrase Equation (27(27) $\begin{array}{r}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)\simeq \sum _{j=0}^N\left(\sum _{k=\lceil \alpha \rceil }^i{\phi }_{(i,j,k)}\right)L_j\left(t\right),i=\lceil \alpha \rceil ,\dots ,N.\end{array}$(27) ) in a vectorized form, resulting in the following representation (28) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)& \simeq [\sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,0,k},\sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,1,k},\sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,2,k},,\dots ,\\ & \sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,N,k}]\psi \left(t\right),i=\lceil \alpha \rceil ,\dots ,N.\end{array}$(28) By consulting Lemma 4.1, we derive the subsequent result (29) $\begin{array}{r}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)\simeq [0,0,\dots ,0]\psi \left(t\right),i=0,1,2,\dots ,\lceil \alpha \rceil -1.\end{array}$(29) The desired result can be obtained by combining Equations (28(28) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)& \simeq [\sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,0,k},\sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,1,k},\sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,2,k},,\dots ,\\ & \sum _{k=\lceil \alpha \rceil }^i{\phi }_{i,N,k}]\psi \left(t\right),i=\lceil \alpha \rceil ,\dots ,N.\end{array}$(28) ) and (29(29) $\begin{array}{r}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }{L}_i\left(t\right)\simeq [0,0,\dots ,0]\psi \left(t\right),i=0,1,2,\dots ,\lceil \alpha \rceil -1.\end{array}$(29) ).

4.2. Motivation and innovation

In this section, our primary objective is to present more detailed information on the challenges associated with traditional methods relying on the operational matrices of orthogonal polynomials. Spectral Tau methods, in conjunction with operational matrices of orthogonal polynomials is a widely used efficient numerical approach for solving FODEs. This approach involves computing Residual Functions, which limits its applicability due to potential computational and precision challenges. This limitation arises because a large residual value indicates a significant error in the approximation.

To enhance accuracy, it is crucial to minimize the Residual Functions by adjusting approximation parameters, such as the coefficients of the basis functions and the Tau parameter. Additionally, in implementing spectral Tau methods, Residual Functions must be expanded as a series of orthogonal polynomials, and initial or boundary conditions are then applied as constraints. This complexity contributes to the less user-friendly nature of spectral Tau methods.

To reinforce our arguments, we apply the Tau methods proposed in [12,15–17,38] to solve specific examples, including Examples 7.1, 7.2, 7.4, 7.5, and 7.6 by computing their Residual Functions. Notably, achieving promising results requires minimizing Residual Functions by increasing the number of terms of orthogonal polynomials. This, however, escalates the computational cost of the algorithm, as seen in Tables 3 and 5. Similarly, minimizing the error between the exact solution and the Tau method-approximated solution demands increasing the number of terms in the series expansion based on the basis functions of orthogonal polynomials. This raises concerns about the efficiency of approaches that integrate operational matrices and the Tau method, as indicated in Tables 4 and 9.

While Tau methods combined with operational matrices are effective for FODEs with polynomial solutions, they compromise accuracy when tackling problems with non-polynomial solutions, as demonstrated in Example 7.5 and Figure 6. These observations underscore the effectiveness of our PNM when compared to methods blending Tau methods with operational matrices. Importantly, our PNM is independent of the computation choice for obtaining the system of algebraic equations through Residual Functions. Our method transforms FODEs into a system of algebraic equations differently, as exemplified in the solution of Example 7.1. Consequently, we are motivated to introduce an efficient class of operational matrices based approach by incorporating modifications into existing techniques.

5. Application of the operational matrices of LPs

Our prime focus in this section is to discuss the development of the framework of an efficient and effective numerical method to solve FODEs with homogeneous and nonhomogeneous initial conditions and diverse source terms. Our PNM relies on the utilization of FODOM and FOIOM associated with LPs. A notable benefit of the PNM is its ability to convert FODEs into a set of Sylvester-type matrix equations, simplifying the solving process with the application of any computational software.

Consider the following generalized FODEs corresponding to initial conditions (30) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)& =H(t,y(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\right(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\right(t),\dots {,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y(t\left)\right),\\ y^{\left(q\right)}\left(0\right)& =d_q,q=0,1,\dots ,N,\end{array}$(30) where, $N-1<\alpha \le N,0<{\rho }_0<{\rho }_1<\cdots <{\rho }_N<\alpha$, H in general is a nonlinear function, and ${}_{\mathrm{C}}\mathcal{D}$ represents a Caputo fractional-order derivative.

Consider the following holds true (31) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)=F^T\psi \left(t\right).$(31) By using the result of Theorem 4.1 and the initial conditions (30(30) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)& =H(t,y(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\right(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\right(t),\dots {,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y(t\left)\right),\\ y^{\left(q\right)}\left(0\right)& =d_q,q=0,1,\dots ,N,\end{array}$(30) ), we can express Equation (31(31) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)=F^T\psi \left(t\right).$(31) ) as (32) $\begin{array}{r}y\left(t\right)\approx F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+d_0+d_{1}t+\cdots +d_N{t}^N.\end{array}$(32) The expression $d_0+d_{1}t+\cdots +d_N{t}^N$ in Equation (32(32) $\begin{array}{r}y\left(t\right)\approx F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+d_0+d_{1}t+\cdots +d_N{t}^N.\end{array}$(32) ) can be approximated as (33) $d_0+d_{1}t+\cdots +d_N{t}^N\approx \sum _{p=0}^M{\omega }_p{L}_p\left(t\right)={\mathrm{\Omega }}^T\psi \left(t\right).$(33) The expression, ${\mathrm{\Omega }}^T=[{\omega }_0,{\omega }_1,\dots ,{\omega }_M]$ can be determined as (34) $\begin{array}{rl}{\omega }_p& ={\int }_0^{\mathrm{\infty }}(d_0+d_{1}t+\cdots +d_N{t}^N){\mathrm{e}}^{-t}{L}_p\left(t\right)\mathrm{d}t,\\ & p=0,1,\dots ,M.\end{array}$(34) Using Equation (33(33) $d_0+d_{1}t+\cdots +d_N{t}^N\approx \sum _{p=0}^M{\omega }_p{L}_p\left(t\right)={\mathrm{\Omega }}^T\psi \left(t\right).$(33) ) in Equation (32(32) $\begin{array}{r}y\left(t\right)\approx F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+d_0+d_{1}t+\cdots +d_N{t}^N.\end{array}$(32) ), we have (35) $y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(35) Using Theorem 4.2 and Equation (35(35) $y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(35) ), the terms involved in the right-hand side of the problem (30(30) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)& =H(t,y(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\right(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\right(t),\dots {,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y(t\left)\right),\\ y^{\left(q\right)}\left(0\right)& =d_q,q=0,1,\dots ,N,\end{array}$(30) ) can be determined as

(36) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right),\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right),\\ & ⋮⋮⋮⋮\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right),\\ & \mathrm{and}\\ & H\left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).\end{array}\}$(36) Now substituting Equations (31(31) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)=F^T\psi \left(t\right).$(31) ), (35(35) $y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(35) ), and (36(36) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right),\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right),\\ & ⋮⋮⋮⋮\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right),\\ & \mathrm{and}\\ & H\left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).\end{array}\}$(36) ) in Equation (30(30) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)& =H(t,y(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\right(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\right(t),\dots {,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y(t\left)\right),\\ y^{\left(q\right)}\left(0\right)& =d_q,q=0,1,\dots ,N,\end{array}$(30) ), we have (37) $F^T-F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\mathrm{\Phi }={\mathrm{\Omega }}^T\mathrm{\Phi }+{\mathrm{{\rm Y}}}^T,$(37) where $\mathrm{\Phi }=(R_{(N+1,N+1)}^{\left({\rho }_0\right)}+R_{(N+1,N+1)}^{\left({\rho }_1\right)}+\cdots +R_{(N+1,N+1)}^{\left({\rho }_N\right)}+I_{(N+1,N+1)})$. The matrix equation represented by Equation (37(37) $F^T-F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\mathrm{\Phi }={\mathrm{\Omega }}^T\mathrm{\Phi }+{\mathrm{{\rm Y}}}^T,$(37) ) is of the Sylvester type, where the unknown vector is $F^T$. By substituting the values of the unknown vector $F^T$ into Equation (35(35) $y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(35) ), we can determine an approximate solution to the fractional problem described in Equation (30(30) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)& =H(t,y(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\right(t\left){,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\right(t),\dots {,}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y(t\left)\right),\\ y^{\left(q\right)}\left(0\right)& =d_q,q=0,1,\dots ,N,\end{array}$(30) ).

6. Error estimation

Let $w=w\left(t\right)={\mathrm{e}}^{-t}$ represents a weighted function defined on the interval $[0,\mathrm{\infty })$, i.e. a non-negative integrable function within this range. Denote the space $L^2[0,\mathrm{\infty })$ as $L_w^2[0,\mathrm{\infty })$, defined as the set of functions $y:[0,\mathrm{\infty })\to \mathbb{R}$, such that (38) $\begin{array}{rl}{L}_w^2& =L_w^2[0,\mathrm{\infty })\\ & =\{y:[0,\mathrm{\infty })\to \mathbb{R}\mathrm{s}.\mathrm{t}{\int }_0^{\mathrm{\infty }}|y(t\left){|}^{2}w\right(t)\mathrm{d}t<\mathrm{\infty }\}.\end{array}$(38) Consider ${\psi }_N$ as a finite-dimensional and closed subspace of $L_w^2[0,\mathrm{\infty })$, generated by the Laguerre polynomials' bases, given by ${\psi }_N=\{L_0,L_1,\dots ,L_N\}$. For any $y\left(t\right)\in L^{2}w[0,\mathrm{\infty })$, there exists a unique best approximation $y_N\left(t\right)\in {\psi }_N$, such that (39) $||y\left(t\right)-y_N\left(t\right)|{|}_w\le ||y-\nu |{|}_w,\mathrm{forall}\nu \in {\psi }_N.$(39) Furthermore, if $y_N\left(t\right)\in {\psi }_N$, it can be expressed as (40) $y_N\left(t\right)=\sum _{q=0}^N{c}_q{L}_q=F^T\psi \left(t\right).$(40) The series coefficients $c_q$ can be computed using Equation (12(12) $c_i={\int }_0^{\mathrm{\infty }}u\left(t\right){\mathrm{e}}^{-t}{L}_i\left(t\right)\mathrm{d}t.$(12) ).

Lemma 6.1

Suppose ${\psi }_N=\{L_0,L_1,\dots ,L_N\}$ and $y\left(t\right)\in C^N[0,\mathrm{\infty })$, where $y:[0,\mathrm{\infty })\to \mathbb{R}$. If $y_N\left(t\right)$ is the best approximation out of ${\psi }_N,i.e.y_N\left(t\right)=\sum _{q=0}^N{c}_q{L}_q=F^T\psi \left(t\right)$ then the error bound can be computed by using the following relation (41) $||y\left(t\right)-y_N\left(t\right)|{|}_w\le \frac{K}{\mathrm{\Gamma }(N+1)}{\mathrm{e}}^b\mathrm{\Gamma }(2N+1,b).$(41) where $K=\underset{x\in \mathbb{[}0,\mathrm{\infty })}{max}\left(y^{\left(N\right)}\right(t\left)\right)$.

Proof.

Suppose that $y\left(t\right)$ is N times continuously differentiable real-valued function, then its power series representation about t = b exists, given as (42) $p\left(t\right)=\sum _{q=0}^N\frac{(t-b{)}^q}{q!}{y}^{\left(q\right)}\left(b\right),$(42) with the following error bound (43) $\begin{array}{rl}|y\left(t\right)-p\left(t\right)|& =|y^{\left(N\right)}\left(\xi \right)\frac{(t-b{)}^N}{N!}|,b<\xi <\mathrm{\infty },\\ & \le K\frac{(t-b{)}^N}{N!}.\end{array}$(43) Since $y_N\left(t\right)$ is the best approximation out of ${\psi }_N$ and $p\left(t\right)\in {\psi }_N$, then we have (44) $\begin{array}{rl}||y\left(t\right)-y_N\left(t\right)|{|}_w^2& \le ||y\left(t\right)-p\left(t\right)|{|}_w^2\\ & ={\int }_0^{\mathrm{\infty }}\left(y\right(t)-p(t\left){)}_w^{2}w\right(t)\mathrm{d}t\\ & \le \frac{K^2}{(N!{)}^2}{\int }_0^{\mathrm{\infty }}(t-b{)}^{2N}w(t)\mathrm{d}t.\end{array}$(44) Given that $w\left(t\right)={\mathrm{e}}^{-t}$ is a non-negative integrable function defined on $[0,\mathrm{\infty })$, and by employing the substitution $(t-b)=u$, the integral on the right side of inequality (44(44) $\begin{array}{rl}||y\left(t\right)-y_N\left(t\right)|{|}_w^2& \le ||y\left(t\right)-p\left(t\right)|{|}_w^2\\ & ={\int }_0^{\mathrm{\infty }}\left(y\right(t)-p(t\left){)}_w^{2}w\right(t)\mathrm{d}t\\ & \le \frac{K^2}{(N!{)}^2}{\int }_0^{\mathrm{\infty }}(t-b{)}^{2N}w(t)\mathrm{d}t.\end{array}$(44) ) takes the form of an incomplete gamma function. Consequently, (44(44) $\begin{array}{rl}||y\left(t\right)-y_N\left(t\right)|{|}_w^2& \le ||y\left(t\right)-p\left(t\right)|{|}_w^2\\ & ={\int }_0^{\mathrm{\infty }}\left(y\right(t)-p(t\left){)}_w^{2}w\right(t)\mathrm{d}t\\ & \le \frac{K^2}{(N!{)}^2}{\int }_0^{\mathrm{\infty }}(t-b{)}^{2N}w(t)\mathrm{d}t.\end{array}$(44) ) can be represented as (45) $||y\left(t\right)-y_N\left(t\right)|{|}_w^2\le \frac{K^2}{\left(\mathrm{\Gamma }\right(N+1){)}^2}{\mathrm{e}}^b\mathrm{\Gamma }(2N+1,b).$(45) Equation (45(45) $||y\left(t\right)-y_N\left(t\right)|{|}_w^2\le \frac{K^2}{\left(\mathrm{\Gamma }\right(N+1){)}^2}{\mathrm{e}}^b\mathrm{\Gamma }(2N+1,b).$(45) ) can further be expressed as (46) $||y\left(t\right)-y_N\left(t\right)|{|}_w\le \frac{K}{\mathrm{\Gamma }(N+1)}{\mathrm{e}}^b\mathrm{\Gamma }(2N+1,b).$(46)

7. Examples

In this section, we assess the effectiveness of the PNM by solving several illustrative examples. The results obtained are then compared with those previously reported in the literature. We present these comparative findings using both tables and plots.

Example 7.1

Consider the following FODEs of Bagley–Torvik type associated with the initial conditions [15,19–21] (47) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(47) The exact solution of Example 7.1 at $\alpha =2$, $d_0=1$, $d_1=1$, and $H\left(t\right)=1+t$ exists, given as (48) $y\left(t\right)=1+t.$(48) By applying the algorithm described in Section 5, with N = 2, we have (49) ${}_{\mathrm{C}}{\mathcal{D}}^{2}y\left(t\right)=F^T\psi \left(t\right).$(49) Using Equation (32(32) $\begin{array}{r}y\left(t\right)\approx F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)+d_0+d_{1}t+\cdots +d_N{t}^N.\end{array}$(32) ), we can express Equation (49(49) ${}_{\mathrm{C}}{\mathcal{D}}^{2}y\left(t\right)=F^T\psi \left(t\right).$(49) ) in the following way (50) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+1+t.$(50) The expression, 1 + t in Equation (50(50) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+1+t.$(50) ) can be approximated as (51) $1+t\approx \sum _{p=0}^2{\omega }_p{L}_p\left(t\right)={\mathrm{\Omega }}^T\psi \left(t\right).$(51) The series coefficients, ${\mathrm{\Omega }}^T=[{\omega }_0,{\omega }_1,{\omega }_2]$ can be determined as (52) ${\omega }_p={\int }_0^{\mathrm{\infty }}(1+t){\mathrm{e}}^{-t}{L}_p\left(t\right)\mathrm{d}t,p=0,1,2.$(52) Using Equation (51(51) $1+t\approx \sum _{p=0}^2{\omega }_p{L}_p\left(t\right)={\mathrm{\Omega }}^T\psi \left(t\right).$(51) ), we can express Equation (50(50) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+1+t.$(50) ) as (53) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(53) Using Equation (36(36) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right),\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right),\\ & ⋮⋮⋮⋮\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right),\\ & \mathrm{and}\\ & H\left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).\end{array}\}$(36) ), we can represent the term ${}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)$ in Equation (47(47) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(47) ) as (54) ${}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0=\frac{3}{2}}y\left(t\right)=F^T{Q}_{(3,3)}^{\left(2\right)}{R}_{(3,3)}^{\left(\frac{3}{2}\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(3,3)}^{\left(\frac{3}{2}\right)}\psi \left(t\right).$(54) The expression on the right-hand side, $H\left(t\right)=1+t$, in Equation (47(47) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(47) ) can be approximated as (55) $1+t\approx \sum _{q=0}^2{\upsilon }_q{L}_q\left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).$(55) The series coefficients, ${\mathrm{{\rm Y}}}^T=[{\upsilon }_0,{\upsilon }_1,{\upsilon }_2]$ can be determined as (56) ${\upsilon }_q={\int }_0^{\mathrm{\infty }}(1+t){\mathrm{e}}^{-t}{L}_q\left(t\right)\mathrm{d}t,q=0,1,2.$(56) Now substituting, Equations (49(49) ${}_{\mathrm{C}}{\mathcal{D}}^{2}y\left(t\right)=F^T\psi \left(t\right).$(49) ), (53(53) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(53) ), (54(54) ${}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0=\frac{3}{2}}y\left(t\right)=F^T{Q}_{(3,3)}^{\left(2\right)}{R}_{(3,3)}^{\left(\frac{3}{2}\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(3,3)}^{\left(\frac{3}{2}\right)}\psi \left(t\right).$(54) ), and (55(55) $1+t\approx \sum _{q=0}^2{\upsilon }_q{L}_q\left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).$(55) ) into Equation (47(47) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(47) ), we have (57) $\begin{array}{rl}& F^T\psi \left(t\right)+F^T{Q}_{(3,3)}^{\left(2\right)}{R}_{(3,3)}^{\left(\frac{3}{2}\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(3,3}^{\left(\frac{3}{2}\right)}\psi \left(t\right)\\ & +F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).\end{array}$(57) Following simplifications, we can express Equation (57(57) $\begin{array}{rl}& F^T\psi \left(t\right)+F^T{Q}_{(3,3)}^{\left(2\right)}{R}_{(3,3)}^{\left(\frac{3}{2}\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(3,3}^{\left(\frac{3}{2}\right)}\psi \left(t\right)\\ & +F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).\end{array}$(57) ) as: (58) $\begin{array}{rl}& F^T+F^T{Q}_{(3,3)}^{\left(2\right)}(R_{(3,3)}^{\left(\frac{3}{2}\right)}+I_{(3,3)})\\ & ={\mathrm{{\rm Y}}}^T-{\mathrm{\Omega }}^T(R_{(3,3)}^{\left(\frac{3}{2}\right)}+I_{(3,3)}).\end{array}$(58) Equation (58(58) $\begin{array}{rl}& F^T+F^T{Q}_{(3,3)}^{\left(2\right)}(R_{(3,3)}^{\left(\frac{3}{2}\right)}+I_{(3,3)})\\ & ={\mathrm{{\rm Y}}}^T-{\mathrm{\Omega }}^T(R_{(3,3)}^{\left(\frac{3}{2}\right)}+I_{(3,3)}).\end{array}$(58) ) represents a Sylvester-type equation, and the unknown vector $F^T$ can be easily computed in MATLAB by employing the following expressions. (59) $\begin{array}{rl}{Q}_{(3,3)}^{\left(2\right)}& =\left[\begin{array}{ccc}1& -2& 1\\ 0& 1& -2\\ 0& 0& 1\end{array}\right],\\ R_{(3,3)}^{\left(\frac{3}{2}\right)}& =\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 1.0000& -0.5000& -0.1250\end{array}\right],\\ {\mathrm{\Omega }}^T& =\left[\begin{array}{ccc}2& -1& 0\end{array}\right],\\ {\mathrm{{\rm Y}}}^T& =\left[\begin{array}{ccc}2& -1& 0\end{array}\right].\end{array}$(59) By incorporating the results from Equation (59(59) $\begin{array}{rl}{Q}_{(3,3)}^{\left(2\right)}& =\left[\begin{array}{ccc}1& -2& 1\\ 0& 1& -2\\ 0& 0& 1\end{array}\right],\\ R_{(3,3)}^{\left(\frac{3}{2}\right)}& =\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 1.0000& -0.5000& -0.1250\end{array}\right],\\ {\mathrm{\Omega }}^T& =\left[\begin{array}{ccc}2& -1& 0\end{array}\right],\\ {\mathrm{{\rm Y}}}^T& =\left[\begin{array}{ccc}2& -1& 0\end{array}\right].\end{array}$(59) ) into Equation (58(58) $\begin{array}{rl}& F^T+F^T{Q}_{(3,3)}^{\left(2\right)}(R_{(3,3)}^{\left(\frac{3}{2}\right)}+I_{(3,3)})\\ & ={\mathrm{{\rm Y}}}^T-{\mathrm{\Omega }}^T(R_{(3,3)}^{\left(\frac{3}{2}\right)}+I_{(3,3)}).\end{array}$(58) ), we obtain $F^T=[0,0,0]$. Subsequently, applying this result in Equation (53(53) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(53) ) yields the solution for Example 7.1 as follows: (60) $\begin{array}{rl}y\left(t\right)& =\left[\begin{array}{ccc}2& -1& 0\end{array}\right]\left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\end{array}\right]\\ & =1+t.\end{array}$(60) Because of the non-homogeneous initial conditions in Example 7.1, we can't use Podlubny's method to find its numerical solution. To show how our PNM compares to other methods, we solve Example 7.1 for two different choices of ${\rho }_0:$ 1.5 and 1.75, as seen in Tables 1 and 2. For ${\rho }_0=1.5$, comparing the maximum absolute errors, GA-PS [19], PSO-PS [20], and NN-IPA [21] report values of $4.62\times {10}^{-2}$, $1.95\times {10}^{-2}$, and $2.90\times {10}^{-9}$, respectively. In our case, the maximum absolute error is 0, highlighting the effectiveness of our proposed approach, as shown in Table 1. For ${\rho }_0=1.75$, our PNM also gives good results. In Table 2, we present the exact solution values of Example 7.1 for $0.1\le t\le 1$, along with absolute error values reported in NN–IPA [21] and computed using our PNM with N = 5, N = 10, and N = 20. NN–IPA [21] reports a maximum absolute error of $6.59\times {10}^{-9}$, with error values ranging from ${10}^{-9}$ to ${10}^{-10}$. In our case, the maximum absolute error is 0, showing that our PNM aligns perfectly with the exact solution at different values of N while keeping ${\rho }_0=1.75$ constant, as shown in Table 2. Table 3 compares the execution times (in seconds) of the PNM and an algorithm from [12] for Example 7.1 with ${\rho }_0=1.5$. The table includes approximate solutions and execution times for different values of N (5, 10, 15, and 20). The results highlight the computational efficiency of the PNM in comparison to the algorithm proposed by [12]. Overall, our PNM yields significantly better results compared to other numerical solvers.

Table 1. Comparing the absolute errors of the results, for Example 7.1 computed using PNM with those obtained from the literature, at ${\rho }_0=1.5$.

t	$y\left(t\right)$	GA-PS [19]	PSO-PS [20]	NN-IPA [21]	PNM, N = 2
0.0	1.0	$1.60\times {10}^{-2}$	$2.84\times {10}^{-3}$	$4.90\times {10}^{-10}$	0.00
0.1	1.1	$4.73\times {10}^{-3}$	$3.40\times {10}^{-3}$	$7.97\times {10}^{-10}$	0.00
0.2	1.2	$1.95\times {10}^{-4}$	$4.52\times {10}^{-3}$	$7.42\times {10}^{-10}$	0.00
0.3	1.3	$6.66\times {10}^{-4}$	$2.02\times {10}^{-3}$	$6.70\times {10}^{-10}$	0.00
0.4	1.4	$1.62\times {10}^{-3}$	$2.55\times {10}^{-3}$	$8.56\times {10}^{-10}$	0.00
0.5	1.5	$4.97\times {10}^{-3}$	$7.59\times {10}^{-3}$	$1.37\times {10}^{-9}$	0.00
0.6	1.6	$7.42\times {10}^{-3}$	$2.84\times {10}^{-2}$	$2.08\times {10}^{-9}$	0.00
0.7	1.7	$6.70\times {10}^{-3}$	$2.84\times {10}^{-2}$	$2.69\times {10}^{-9}$	0.00
0.8	1.8	$1.27\times {10}^{-5}$	$2.84\times {10}^{-3}$	$2.90\times {10}^{-9}$	0.00
0.9	1.9	$1.62\times {10}^{-2}$	$2.84\times {10}^{-3}$	$2.61\times {10}^{-9}$	0.00
1.0	2.0	$4.62\times {10}^{-2}$	$2.84\times {10}^{-2}$	$2.14\times {10}^{-9}$	0.00

Table 2. Comparing the absolute errors of the results, for Example 7.1 computed using PNM with those obtained from the literature, at ${\rho }_0=1.75$.

t	$y\left(t\right)$	NN-IPA [21]	PNM, N = 5	PNM, N = 10	PNM, N = 20
0.1	1.1	$3.74\times {10}^{-10}$	0.00	0.00	0.00
0.2	1.2	$1.74\times {10}^{-9}$	0.00	0.00	0.00
0.3	1.3	$2.38\times {10}^{-9}$	0.00	0.00	0.00
0.4	1.4	$2.49\times {10}^{-9}$	0.00	0.00	0.00
0.5	1.5	$2.76\times {10}^{-9}$	0.00	0.00	0.00
0.6	1.6	$3.58\times {10}^{-9}$	0.00	0.00	0.00
0.7	1.7	$4.85\times {10}^{-9}$	0.00	0.00	0.00
0.8	1.8	$6.04\times {10}^{-9}$	0.00	0.00	0.00
0.9	1.9	$6.59\times {10}^{-9}$	0.00	0.00	0.00
1.0	2.0	$6.55\times {10}^{-9}$	0.00	0.00	0.00

Example 7.2

Consider the following FODEs associated with the homogeneous initial conditions [16] (61) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\frac{3}{2}}y\left(t\right)+3y\left(t\right)=H\left(t\right),t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(61) The exact solution of Example 7.2 at $d_0=0=d_1$ and $H\left(t\right)=3t^3+\frac{8}{\mathrm{\Gamma }\left(0.5\right)}{t}^{1.5}$ exists, given as (62) $y\left(t\right)=t^3.$(62) We apply both the PNM and the operational matrices approach along with the Tau method to solve Example 7.2, as proposed in [15–17]. We find that our PNM gives an approximate solution that coincides with the exact solution for different values of N, resulting in a maximum absolute error of 0, see Table 4. On the other hand, when we use the operational matrices approach along with the Tau method, the maximum errors are $1.19808\times {10}^{-17}$, $9.95\times {10}^{-18}$, and $9.82\times {10}^{-18}$ at N = 3, N = 5, and N = 8, respectively, see Table 4. This clearly shows that our PNM provides significantly better results compared to the operational matrices approach along with the Tau method for Example 7.2. Table 5 compares the execution times (in seconds) of our PNM and an algorithm from [12]. The table shows approximate solutions and execution times for different values of N (3, 5, 10, 13, and 15). The results indicating that our PNM is computationally more efficient than the algorithm proposed by [12]. In summary, our PNM produces significantly better results compared to other numerical methods.

Example 7.3

Consider the following FODEs of Bagley–Torvik type associated with the homogeneous initial conditions and the generalized source term that appeared when studying the dynamics of a thin, large plate immersed in a Newtonian fluid [17,39] (63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) The FODEs of Bagley-Torvik type described by Equation (63(63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) ) were first introduced in [39] and thoroughly investigated in [18, Section 8.3.2]. In his book, Podlubny computed the analytical solution for these equations under homogeneous initial conditions by solving the following integral equation (64) $y\left(t\right)={\int }_0^{t}G(t-z)H\left(z\right)\mathrm{d}z.$(64) To find the analytical solution to the problem (63(63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) ), we need to compute the convolution integral that includes Green's function and the source term H. However, if H is a general function, this integral is difficult to compute (see [18] and [24]) because it involves the multivariate Mittag-Leffler functions. Consequently, we prefer to solve the problem using numerical methods.

Using our PNM, we solve the problem (63(63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) ) with various source terms denoted as $H\left(t\right)$. The absolute error of the problem (63(63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) ) is computed using our PNM, Podlubny's approach [18], VIM [23], GA–PS [19], PSO–PS [20], LeNN [22], and NN–IPA [21]. The computed absolute error values are listed in Tables 6 and 7. Additionally, Table 8 presents the results obtained by employing the PNM with various $H\left(t\right)$. Upon analysis, we observe that the maximum absolute error obtained with our PNM for N = 3, N = 10, ${\rho }_0=1.5$, and the source term $H\left(t\right)=t^2+4\sqrt{\frac{t}{\pi }}+2$ is $9.62\times {10}^{-15}$ and $3.34\times {10}^{-16}$, respectively. In contrast, when using Podlubny's approach [18], VIM [23], GA-PS [19], PSO-PS [20], LeNN [22], and NN-IPA [21], the maximum absolute error values are $9.12\times {10}^{-5}$, $2.25\times {10}^{-1}$, $3.43\times {10}^{-2}$, $2.12\times {10}^{-2}$, $5.74\times {10}^{-7}$, and $2.03\times {10}^{-5}$, respectively. Promising results are also obtained for the problem (63(63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) ) when computed with the following source terms: $H\left(t\right)=t^3+6t+\frac{128}{15\mathrm{\Gamma }(\frac{1}{4})}{t}^{\frac{9}{4}}$, $H\left(t\right)=t^3+6t+\frac{8}{\mathrm{\Gamma }\left(0.5\right)}{t}^{\frac{3}{2}}$, and $H\left(t\right)=C(1+t)$, where $C\in \mathbb{R}$, see Table 8. Furthermore, the effectiveness of the PNM is showcased through the solution of problem (63(63) $\begin{array}{rl}& a_1{}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+a_2{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+a_{3}y\left(t\right)=H\left(t\right),\\ & 1<{\rho }_0<2,t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(63) ) across different source terms and various values of N, as illustrated in Figure 3. Therefore, in general, the utilization of our PNM yields highly accurate results compared to other numerical algorithms.

Figure 3. Plots showing the approximate solution and exact solution for Example 7.3 are displayed for different values of N and various source terms.

Table 3. Comparing the results of Example 7.1 based on Algorithm Execution Time(s), as measured by both the PNM and the algorithm proposed in [12], with ${\rho }_0=1.5$.

	Approximate Solution		Execution Time (s)
N	PNM	Algorithm [12]	PNM	Algorithm [12]
5	1 + t	1 + t	9.42	24.08
10	1 + t	1 + t	31.16	261.19
15	1 + t	1 + t	62.21	1427.78
20	1 + t	1 + t	1314.67	5467.58

Table 4. Comparing the absolute errors of the results, for Example 7.2 computed using PNM and the operational matrices approach along with the Tau method proposed in [15–17].

	PNM		Tau Method [15–17]
t	N = 3, 5, 8	N = 3	N = 5	N = 8
0.0	0	0.0000	0.0000	0.0000
0.2	0	$5.04\times {10}^{-19}$	$2.61\times {10}^{-19}$	$9.26\times {10}^{-20}$
0.4	0	$1.99\times {10}^{-18}$	$1.22\times {10}^{-18}$	$7.88\times {10}^{-19}$
0.6	0	$4.42\times {10}^{-18}$	$3.07\times {10}^{-18}$	$2.52\times {10}^{-18}$
0.8	0	$7.76\times {10}^{-18}$	$5.95\times {10}^{-18}$	$5.51\times {10}^{-18}$
1.0	0	$1.19808\times {10}^{-17}$	$9.95\times {10}^{-18}$	$9.82\times {10}^{-18}$

Table 5. Comparing the results of Example 7.2 based on Algorithm Execution Time(s), as measured by both the PNM and the algorithm proposed in [12].

	Approximate Solution		Execution Time (s)
N	PNM	Algorithm [12]	PNM	Algorithm [12]
3	$t^3$	$t^3$	4.87	7.34
5	$t^3$	$t^3$	7.49	21.16
10	$t^3$	$t^3$	19.91	244.98
13	$t^3$	$t^3$	32.65	759.52
15	$t^3$	$t^3$	43.72	1384.75

Table 6. Comparing the absolute error of Example 7.3 outcomes, where the source term is given by $H\left(t\right)=t^2+4\sqrt{\frac{t}{\pi }}+2$ and ${\rho }_0$ is 1.5.

t	$y\left(t\right)=t^2$	Podlubny [18]	VIM [23]	PNM, N = 3
0.0	0.00	0.000000	0.000000	0.000000
0.1	0.01	$5.76\times {10}^{-5}$	$5.48\times {10}^{-5}$	$1.31\times {10}^{-15}$
0.2	0.04	$8.29\times {10}^{-5}$	$6.31\times {10}^{-4}$	$8.40\times {10}^{-16}$
0.3	0.09	$9.12\times {10}^{-5}$	$2.66\times {10}^{-3}$	$4.22\times {10}^{-16}$
0.4	0.16	$8.74\times {10}^{-5}$	$7.48\times {10}^{-3}$	$5.40\times {10}^{-17}$
0.5	0.25	$7.42\times {10}^{-5}$	$1.67\times {10}^{-2}$	$2.67\times {10}^{-16}$
0.6	0.36	$5.36\times {10}^{-5}$	$3.22\times {10}^{-2}$	$5.42\times {10}^{-16}$
0.7	0.49	$2.68\times {10}^{-5}$	$5.80\times {10}^{-2}$	$7.73\times {10}^{-16}$
0.8	0.64	$5.07\times {10}^{-6}$	$9.58\times {10}^{-2}$	$9.62\times {10}^{-15}$
0.9	0.81	$4.12\times {10}^{-5}$	$1.50\times {10}^{-1}$	$1.11\times {10}^{-15}$
1.0	1.00	$8.08\times {10}^{-5}$	$2.25\times {10}^{-1}$	$1.22\times {10}^{-15}$

Table 7. Comparing the absolute error of Example 7.3 outcomes, where the source term is given by $H\left(t\right)=t^2+4\sqrt{\frac{t}{\pi }}+2$ and ${\rho }_0$ is 1.5.

t	GA-PS [19]	PSO-PS [20]	LeNN [22]	NN-IPA [21]	PNM, N = 10
0.0	$3.34\times {10}^{-2}$	$2.12\times {10}^{-2}$	0.000000	$7.63\times {10}^{-6}$	0.000000
0.1	$3.43\times {10}^{-2}$	$2.20\times {10}^{-3}$	$7.59\times {10}^{-8}$	$8.73\times {10}^{-6}$	$4.83\times {10}^{-17}$
0.2	$3.33\times {10}^{-2}$	$2.63\times {10}^{-3}$	$1.57\times {10}^{-7}$	$1.12\times {10}^{-5}$	$2.12\times {10}^{-16}$
0.3	$3.04\times {10}^{-2}$	$2.98\times {10}^{-3}$	$2.09\times {10}^{-7}$	$1.08\times {10}^{-5}$	$3.24\times {10}^{-16}$
0.4	$2.57\times {10}^{-2}$	$2.97\times {10}^{-3}$	$2.74\times {10}^{-7}$	$8.19\times {10}^{-6}$	$3.34\times {10}^{-16}$
0.5	$1.96\times {10}^{-2}$	$2.46\times {10}^{-3}$	$3.61\times {10}^{-7}$	$7.06\times {10}^{-6}$	$2.74\times {10}^{-16}$
0.6	$1.26\times {10}^{-2}$	$1.49\times {10}^{-3}$	$4.37\times {10}^{-7}$	$1.01\times {10}^{-5}$	$1.73\times {10}^{-16}$
0.7	$5.49\times {10}^{-3}$	$2.67\times {10}^{-4}$	$4.68\times {10}^{-7}$	$1.60\times {10}^{-5}$	$5.09\times {10}^{-17}$
0.8	$8.80\times {10}^{-4}$	$8.34\times {10}^{-4}$	$4.64\times {10}^{-7}$	$2.03\times {10}^{-5}$	$7.42\times {10}^{-17}$
0.9	$5.42\times {10}^{-3}$	$1.27\times {10}^{-3}$	$4.87\times {10}^{-7}$	$1.86\times {10}^{-5}$	$1.91\times {10}^{-16}$
1.0	$6.91\times {10}^{-3}$	$3.05\times {10}^{-4}$	$5.74\times {10}^{-7}$	$1.24\times {10}^{-5}$	$2.89\times {10}^{-16}$

Table 8. Absolute error values of Example 7.3 at various source terms, $H\left(t\right)$.

	$y\left(t\right)=t^3$, ${\rho }_0=\frac{3}{4}$, N = 15	$y\left(t\right)=t^3$, ${\rho }_0=\frac{3}{2}$, N = 15	$y\left(t\right)=C(1+t)$, ${\rho }_0=\frac{7}{4}$, N = 25
t	$H\left(t\right)=t^3+6t+\frac{128}{15\mathrm{\Gamma }\left(\frac{1}{4}\right)}{t}^{\frac{9}{4}}$	$H\left(t\right)=t^3+6t+\frac{8}{\mathrm{\Gamma }\left(0.5\right)}{t}^{\frac{3}{2}}$	$H\left(t\right)=C(1+t)$, $C\in \mathbb{R}$
0.0	0.00000	0.00000	0.00000
0.2	$2.08\times {10}^{-14}$	$1.05\times {10}^{-15}$	0.00000
0.4	$9.21\times {10}^{-15}$	$1.81\times {10}^{-15}$	0.00000
0.6	$1.28\times {10}^{-14}$	$2.86\times {10}^{-15}$	0.00000
0.8	$2.35\times {10}^{-14}$	$2.85\times {10}^{-15}$	0.00000
1.0	$1.99\times {10}^{-14}$	$2.32\times {10}^{-15}$	0.00000

Table 9. Comparing the absolute errors of the results, for Example 7.4 computed using PNM and the operational matrices approach, along with the Tau method proposed in [15–17].

	PNM		Tau Method
t	N = 7	N = 7	N = 10	N = 12
0.0	$3.46\times {10}^{-11}$	$3.21\times {10}^{-10}$	$2.08\times {10}^{-9}$	$3.1\times {10}^{-9}$
0.2	$4.49\times {10}^{-11}$	$1.09\times {10}^{-10}$	$1.47\times {10}^{-10}$	$6.00\times {10}^{-11}$
0.4	$4.84\times {10}^{-11}$	$3.28\times {10}^{-10}$	$9.86\times {10}^{-10}$	$1.33\times {10}^{-9}$
0.6	$3.29\times {10}^{-11}$	$1.19\times {10}^{-10}$	$9.69\times {10}^{-10}$	$1.60\times {10}^{-9}$
0.8	$5.08\times {10}^{-12}$	$1.59\times {10}^{-10}$	$4.98\times {10}^{-10}$	$1.11\times {10}^{-9}$
1.0	$2.91\times {10}^{-11}$	$1.65\times {10}^{-10}$	$1.36\times {10}^{-10}$	$4.18\times {10}^{-9}$

Example 7.4

Consider the following multi-order FODEs associated with the homogeneous initial conditions [16] (65) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right){+}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)-2_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\left(t\right)+y\left(t\right)\\ & =H\left(t\right),t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(65) The exact solution of Example 7.4 at $\alpha =2$, ${\rho }_0=\frac{1}{2}$, ${\rho }_1=1$, $d_0=0=d_1$, and $H\left(t\right)=4t-t^2-\frac{6775690661951975t^{\frac{3}{2}}}{4503599627370496}+42t^5-14t^6+t^7+\frac{1516238469807435t^{\frac{13}{2}}}{562949953421312}-2$ exists, given as (66) $y\left(t\right)=t^7-t^2.$(66) We solve Example 7.4 by applying both the PNM and the operational matrices approach, along with the Tau method as proposed in [15–17]. We observe that the approximate solution obtained by using our PNM at N = 7 yields a maximum absolute error of $4.84\times {10}^{-11}$, see Table 9. In contrast, when we employ the operational matrices approach, along with the Tau method, we encounter maximum absolute errors of $3.28\times {10}^{-10}$, $2.08\times {10}^{-9}$, and $4.18\times {10}^{-9}$ at N = 7, N = 10, and N = 12, respectively (refer to Table 9). The effectiveness of our PNM can also be seen in its ability to solve Example 7.4 across different values of N. We notice a great resemblance between the approximate solution and the exact solution at all these values, as shown in Figure 4. These findings and observations clearly emphasize the superiority of our PNM over the operational matrices approach, along with the Tau method for Example 7.4. The error bounds and the weighted norms are also computed for various values of N. The error bounds decrease consistently and approach zero, illustrating a similar trend for the weighted norm. Notably, these results highlight the effectiveness and efficiency of the PNM, as it consistently converges towards exact solution of the problem, see Figure 5.

Figure 4. Plots showing the approximate solution and exact solution for Example 7.4 are displayed for different values of N.

Figure 5. Error estimation analysis for Example 7.4 at different values of N. (a) Error bound evaluation for various N and (b) Weighted norm evaluation for various N.

Figure 6. Plots showing the approximate solution and exact solution for Example 7.5 are displayed for different values of N. (a) PNM and (b) The method in [12].

Example 7.5

Consider the following multi-order FODEs associated with the non-homogeneous initial conditions (67) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+4_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+3y\left(t\right)=H\left(t\right),t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(67) The exact solution of Example 7.5 at $\alpha =2$, $0<{\rho }_0\le 1$, $d_0=3$, $d_1=4$, and $H\left(t\right)=0$ exists, given as (68) $y\left(t\right)=\frac{13{\mathrm{e}}^{-t}}{2}-\frac{7{\mathrm{e}}^{-3t}}{2}.$(68) We address Example (7.5) using both the PNM and the operational matrices approach in conjunction with the Tau method, as proposed in [12]. To assess the effectiveness of our PNM over the method proposed in [12], we solve Example (7.5) for various values of N. It is observed that as N increases, the approximate solution consistently approaches the exact solution of the problem. In contrast, as N increases, the approximate solution exhibits a divergent behavior in [12], as illustrated in Figure 6(b). The error bounds and the weighted norms are also computed for various values of N. The error bounds decrease consistently and approach zero, illustrating a similar trend for the weighted norm. Notably, these results highlight the effectiveness and efficiency of the PNM, as it consistently converges towards exact solution of the problem, see Figure 7.

Figure 7. Error estimation analysis for Example 7.5 at different values of N. (a) Error bound evaluation for various N and (b) Weighted norm evaluation for various N.

Example 7.6

Consider the following FODEs of Bagley–Torvik type associated with the initial conditions [38] (69) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)+{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right)=H\left(t\right),t\in [0,\mathrm{\infty }),\\ & y\left(0\right)=d_0,y^{\prime }\left(0\right)=d_1.\end{array}$(69) The exact solution of Example 7.6 at $\alpha =2$, $d_0=0$, $d_1=0$, ${\rho }_0=\frac{3}{2}$, and $H\left(t\right)=t^5-t^4-12t^2-\frac{4065414397171185t^{\frac{5}{2}}}{562949953421312}+\frac{1451933713275423t^{\frac{7}{2}}}{140737488355328}$ exists, given as (70) $y\left(t\right)=t^5-t^4,$(70) and the exact solution for $H\left(t\right)=6t+t^3+\frac{5081767996463981t^{\frac{3}{2}}}{1125899906842624}$, given as (71) $y\left(t\right)=t^3.$(71) Applying Theorems 4.1 and 4.2, the algorithm outlined in Section 5, and selecting $H\left(t\right)=6t+t^3+\frac{5081767996463981t^{\frac{3}{2}}}{1125899906842624}$ with N = 3, the result is as follows (72) $\begin{array}{rl}{Q}_{(4,4)}^{\left(2\right)}& =\left[\begin{array}{cccc}1& -2& 1& 0\\ 0& 1& -2& 1\\ 0& 0& 1& -2\\ 0& 0& 0& 1\end{array}\right],\end{array}$(72) (73) $\begin{array}{rl}{R}_{(4,4)}^{\left(1.5\right)}& =\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1.0000& -0.5000& -0.1250& -0.0625\\ 2.0000& 0& -0.7500& -0.2500\end{array}\right],\end{array}$(73) (74) $\begin{array}{rl}{c}_0& =6,c_1=-6,c_2=0,c_3=0,\end{array}$(74) (75) $\begin{array}{rl}{\omega }_0& ={\omega }_1={\omega }_2={\omega }_3=0,\end{array}$(75) and (76) $\psi \left(t\right)=\left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ {\displaystyle -\frac{t^3}{6}+\frac{3t^2}{2}-3t+1}\end{array}\right].$(76) By substituting the outputs of Equations (72(72) $\begin{array}{rl}{Q}_{(4,4)}^{\left(2\right)}& =\left[\begin{array}{cccc}1& -2& 1& 0\\ 0& 1& -2& 1\\ 0& 0& 1& -2\\ 0& 0& 0& 1\end{array}\right],\end{array}$(72) ), (74(74) $\begin{array}{rl}{c}_0& =6,c_1=-6,c_2=0,c_3=0,\end{array}$(74) ), (75(75) $\begin{array}{rl}{\omega }_0& ={\omega }_1={\omega }_2={\omega }_3=0,\end{array}$(75) ), and (76(76) $\psi \left(t\right)=\left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ {\displaystyle -\frac{t^3}{6}+\frac{3t^2}{2}-3t+1}\end{array}\right].$(76) ) into Equation (53(53) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(53) ), we have

(77) $\begin{array}{rl}y\left(t\right)& =\left(\left[\begin{array}{cccc}6& -6& 0& 0\end{array}\right]\left[\begin{array}{cccc}1& -2& 1& 0\\ 0& 1& -2& 1\\ 0& 0& 1& -2\\ 0& 0& 0& 1\end{array}\right]\right)\\ & \times \left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ {\displaystyle -\frac{t^3}{6}+\frac{3t^2}{2}-3t+1}\end{array}\right]\\ & =\left[\begin{array}{cccc}6& -18& 18& -6\end{array}\right]\left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ {\displaystyle -\frac{t^3}{6}+\frac{3t^2}{2}-3t+1}\end{array}\right]\\ & =t^3,\end{array}$(77) which coincides with the exact solution of Example 7.6. For the case, $H\left(t\right)=t^5-t^4-12t^2-\frac{4065414397171185t^{\frac{5}{2}}}{562949953421312}+\frac{1451933713275423t^{\frac{7}{2}}}{140737488355328}$ with N = 5, we have the following result (78) $\begin{array}{rl}{Q}_{(6,6)}^{\left(2\right)}& =\left[\begin{array}{cccccc}1& -2& 1& 0& 0& 0\\ 0& 1& -2& 1& 0& 0\\ 0& 0& 1& -2& 1& 0\\ 0& 0& 0& 1& -2& 1\\ 0& 0& 0& 0& 1& -2\\ 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}$(78) (79) $\begin{array}{rl}{R}_{(6,6)}^{\left(1.5\right)}& =[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 1.0000& -0.5000& -0.1250\\ 2.0000& 0& -0.7500\\ 3.0000& 0.5000& -0.3750\\ 4.0000& 1.0000& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ -0.0625& -0.0391& -0.0273\\ -0.2500& -0.1406& -0.0938\\ -0.9375& -0.3672& -0.2227\\ -0.6250& -1.0938& -0.4766\end{array}],\end{array}$(79) (80) $\begin{array}{rl}{c}_0& =96,c_1=-312,c_2=336,c_3=-120,\\ c_4& =0,c_5=0,\end{array}$(80) (81) $\begin{array}{rl}{\omega }_0& ={\omega }_1={\omega }_2={\omega }_3={\omega }_4={\omega }_5=0,\end{array}$(81) and (82) $\psi \left(t\right)=\left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ {\displaystyle -\frac{t^3}{6}+\frac{3t^2}{2}-3t+1}\\ {\displaystyle \frac{t^4}{24}-\frac{2t^3}{3}+3t^2-4t+1}\\ {\displaystyle -\frac{t^5}{120}+\frac{5t^4}{24}-\frac{5t^3}{3}+5t^2-5t+1}\end{array}\right].$(82) Substituting the outputs of Equations (78(78) $\begin{array}{rl}{Q}_{(6,6)}^{\left(2\right)}& =\left[\begin{array}{cccccc}1& -2& 1& 0& 0& 0\\ 0& 1& -2& 1& 0& 0\\ 0& 0& 1& -2& 1& 0\\ 0& 0& 0& 1& -2& 1\\ 0& 0& 0& 0& 1& -2\\ 0& 0& 0& 0& 0& 0\end{array}\right]\end{array}$(78) ), (80(80) $\begin{array}{rl}{c}_0& =96,c_1=-312,c_2=336,c_3=-120,\\ c_4& =0,c_5=0,\end{array}$(80) ), (81(81) $\begin{array}{rl}{\omega }_0& ={\omega }_1={\omega }_2={\omega }_3={\omega }_4={\omega }_5=0,\end{array}$(81) ), and (82(82) $\psi \left(t\right)=\left[\begin{array}{c}1\\ 1-t\\ {\displaystyle \frac{t^2}{2}-2t+1}\\ {\displaystyle -\frac{t^3}{6}+\frac{3t^2}{2}-3t+1}\\ {\displaystyle \frac{t^4}{24}-\frac{2t^3}{3}+3t^2-4t+1}\\ {\displaystyle -\frac{t^5}{120}+\frac{5t^4}{24}-\frac{5t^3}{3}+5t^2-5t+1}\end{array}\right].$(82) ) into Equation (53(53) $y\left(t\right)\approx F^T{Q}_{(3,3)}^{\left(2\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T\psi \left(t\right).$(53) ), we have $y\left(t\right)=t^5-t^4$, which coincides with the exact solution of Example 7.6.

It's important to mention that the Example 7.6 with the identical source terms have been previously addressed in [38]. However, our approach differs in the sense that we employ PNM, and the approximate solution coincides with the exact solutions of Example 7.6. Specifically, for the case of $H\left(t\right)=6t+t^3+\frac{5081767996463981t^{\frac{3}{2}}}{1125899906842624}$, our PNM exactly matches with the exact solution at N = 3, representing the initial four terms of LPs. In contrast, the corresponding result in [38] is obtained at N = 6, representing the initial six terms of Fermat polynomials.

Similarly, for the case of $H\left(t\right)=t^5-t^4-12t^2-\frac{4065414397171185t^{\frac{5}{2}}}{562949953421312}+\frac{1451933713275423t^{\frac{7}{2}}}{140737488355328}$, our PNM coincides the exact solution at N = 5, representing the first six terms of LPs. This corresponds to the result in [38], where the approximate solution aligns with the exact solution at N = 6, corresponding to the initial six terms of Fermat polynomials. The comparison highlights the effectiveness of our proposed approach and the approach addressed in [38] in achieving accurate solutions with a reduced number of terms in the polynomial expansion. In conclusion, both the approaches are fit to solve the multi-order FODEs. The error analysis part and the numerical algorithm proposed in [38] are worth to study.

8. Comparative analysis of the results

LPs are particularly well-suited for addressing challenges associated with exponential weight functions and semi-infinite domains. In scenarios where the problems at hand possess these specific characteristics, Laguerre polynomials may demonstrate superior convergence properties when compared to Legendre, Chebyshev, and Chelyshkov polynomials, each designed for distinct weight functions or domains. As an illustration, in a previous study [15], the authors utilized the Legendre basis for addressing fractional-order differential equations. However, the results presented in Table 4 indicate that the solution approximated with the Laguerre basis exhibits greater efficiency compared to the one achieved using the Legendre basis in [15]. Similarly, promising outcomes are evident in Table 9, where the Laguerre basis performs better than the Legendre basis.

Moreover, when extending the approximation to a larger domain, Legendre, Chebyshev, Chelyshkov, and other orthogonal polynomials (excluding Hermite polynomials) compromise accuracy. In contrast, the Laguerre basis consistently provides accurate approximations on larger domain sets. To support this assertion, we present the approximation of a rational function using various terms of Legendre polynomials, Chelyshkov polynomials, Chebyshev polynomials, and Laguerre polynomials. Our observations, detailed in Figure 1 confirm that Laguerre polynomials deliver better approximation compared to others.

To illustrate the effectiveness of LPs, we solve Examples 7.1 and 7.6 over a semi-infinite region and find that the approximate solution attained on $[0,\mathrm{\infty })$ aligns with the exact solution of Examples 7.1 and 7.6. Likewise, promising outcomes are achieved for Examples 7.2–7.5 when solved over $[0,\mathrm{\infty })$. For comparative analysis with the results found in the existing literature, we calculate the absolute errors for Examples 7.1–7.4 over $[0,1]$. Furthermore, we extend the findings of Examples 7.3–7.5 to wider domains beyond $[0,1]$, namely $[0,5)$, $[0,10)$, and $[0,5)$, respectively.

In [19], a novel stochastic approach was employed to address the Bagley–Torvik equation for fractional order systems using feed-forward artificial neural networks. These networks were trained with evolutionary computational intelligence, incorporating genetic algorithms and pattern search techniques. The authors computed the absolute errors for Example 7.1 at N = 10, representing the number of neurons, i.e. the initial ten terms of the series involving the exponential function. They considered a fractional order parameter ${\rho }_0=1.5$ with a source term $H\left(t\right)=1+t$.

In our study, we applied our PNM to tackle the same problem with identical parameters and source term for N = 2 (representing the initial three terms of the series involving LPs). The obtained results are summarized in Table 1, where we compare absolute errors with the results obtained in [19]. The notable findings from the table highlight that our PNM, with a minimum number of terms (N = 2), consistently yields zero absolute errors across different time points. In comparison to the well-established method in [19], our proposed technique exhibits superior accuracy, achieving zero errors for each considered time step. These results emphasize the effectiveness and precision of our PNM in approximating solutions for FODEs.

In [21], an efficient computational method was presented for solving fractional order systems with initial value problems in the Bagley–Torvik equations. This approach utilized the fractional neural networks optimized through interior point algorithms (IPAs). The authors calculated the absolute errors for Example 7.1 with N = 10, representing the number of neurons. They considered a fractional order parameter ${\rho }_0=1.5$ with a source term $H\left(t\right)=1+t$. In our investigation, we employed our PNM to address the same problem, maintaining identical parameters and a consistent source term across various values of N. The results are presented in Table 2, allowing a comparison of absolute errors with those obtained in [21]. Notably, our PNM consistently yields zero absolute errors across different time points for various values of N. In contrast to the well-established method in [21], our proposed technique demonstrates superior accuracy, achieving zero errors for each considered time step. These findings underscore the effectiveness and precision of our PNM in approximating solutions for FODEs.

In [15], a hybrid approach, utilizing the fractional-order derivative operational matrix of Legendre polynomials in the Caputo sense along with the Tau method, was applied to solve FODEs. The absolute error for Examples 7.2 and 7.4 was computed using various terms of Legendre polynomials (LPs). The results were compared between the method proposed in [15] and our PNM. Across different values of N (number of terms of Legendre polynomials or LPs), our PNM consistently demonstrated higher accuracy compared to the method in [15]. These comparative results are presented in Tables 4 and 9.

Similarly, in [16], the same approach proposed in [15], but employing the fractional-order integral operational matrix of LPs in the Riemann–Liouville sense, was utilized to address FODEs. The absolute error for Examples 7.2 and 7.4 was computed using different terms of LPs, and the results were compared between the method proposed in [16] and our PNM. Across various values of N (number of terms of LPs), our PNM consistently outperformed the method in [16]. These comparative findings are detailed in Tables 4 and 9.

Furthermore, in [17], a modification of the approach presented in [15] was implemented, using the fractional-order derivative operational matrix of LPs in the Caputo sense to solve FODEs. The absolute error for Examples 7.2 and 7.4 was computed using different terms of LPs, and the results were compared between the method proposed in [17] and our PNM. Across various values of N (number of terms of LPs), our PNM consistently demonstrated superior accuracy compared to the method in [17]. These comparative findings are outlined in Tables 4 and 9.

In Table 3, we compare the execution times of Example 7.1 obtained using our PNM and the algorithm from [12]. The analysis focuses on various values of N, with a fixed parameter ${\rho }_0=1.5$. The algorithm in [12] employed a hybrid approach combining the fractional-order integral operational matrix of LPs in the Riemann–Liouville sense with the Tau method. As illustrated in Table 3, the execution time results for both methods demonstrate distinct performance characteristics. For instance, when N = 5, the PNM exhibits an execution time of 9.42 seconds, while the algorithm in [12] takes 24.08 seconds. As N increases to 20, the execution time of the PNM remains relatively stable, whereas the algorithm in [12] experiences a more substantial increase. Similar comparison is shown for the results of Example 7.2, see Table 5. The results emphasize the effectiveness of our method in addressing FODEs with computational efficiency.

To address the effectiveness of our PNM in handling problems with exponential functions in their solutions, we apply it to Example 7.5. We utilized both the PNM and the operational matrices approach with the Tau method, as introduced in [12]. We solved Example 7.5 for various values of N, where N represents the number of terms of LPs. Notably, as N increases, our PNM consistently converges toward the exact solution, as depicted in Figure 6(b).

In contrast, the method proposed in [12] exhibits a divergent behavior as N increases, as illustrated in the same figure. Additionally, we compute error bounds and weighted norms for various values of N. The error bounds consistently decrease and approach zero, reflecting a similar trend observed in the weighted norm. These comprehensive results highlight the efficiency of the PNM, as it consistently converges towards the exact solution of the problem, as visually represented in Figure 7.

In [22], the authors employed the Legendre artificial neural network to solve the Bagley–Torvik equation. This involved employing the Caputo fractional-order derivative to address FODEs of Bagley–Torvik type. The optimal weights of the network were determined through a simulated annealing optimization technique. The authors computed the absolute error for Example 7.3 at N = 10, where N represents the number of terms in the Legendre artificial neural network, and considered a fractional order parameter ${\rho }_0=1.5$ with a source term $H\left(t\right)=t^2+4\sqrt{\frac{t}{\pi }}+2$. In our study, we tackled the same problem with identical parameters using our PNM. The obtained results are summarized in Table 7, where we compare absolute errors with the method proposed in [22] and our PNM with N = 10 terms. The noteworthy findings from the table highlight that our PNM, with a modest number of terms (N = 10, representing the number of terms in LPs), consistently yields minimal absolute errors across different time points. In comparison to other methods, including the well-established [22], our proposed technique exhibits superior accuracy, achieving nearly negligible errors for each considered time step. These results emphasize the effectiveness and precision of our PNM in approximating solutions for FODEs of Bagley–Torvik type.

In the work by [23], the variational iteration method was introduced for solving nth-order semi-differential equations. This method relies on an iterative formula that incorporates the integrand, including the source terms $H\left(t\right)$, derivative terms, and the Lagrange multiplier λ. To derive the suitable solution, it is necessary to split the source terms into two appropriate functions, and the optimal computation of λ is achieved through variational theory. The authors applied this method to solve Example 7.3, calculating its absolute error with a fractional order parameter ${\rho }_0=1.5$ and a source term $H\left(t\right)=t^2+4\sqrt{\frac{t}{\pi }}+2$ for N = 3, representing the number of iterations of the proposed iterative method of [23]. We also solved the same problem with our PNM with identical fractional-order parameter and the source term at N = 3, denoting the number of Lps. The results are presented in Table 6, providing a comparison of absolute errors with those obtained in [23]. Notably, our PNM for N = 3 consistently yields minimum absolute errors across different time points. In contrast to the established method in [23], our proposed technique demonstrates superior accuracy, achieving minimum errors for each considered time step. These findings underscore the effectiveness and precision of our PNM in approximating solutions for FODEs.

In [38], the author introduced a novel operational matrix for fractional-order derivatives (in the sense of Caputo) applied to Fermat polynomials. This matrix serves as a crucial tool in solving the fractional Bagley-Torvik equation using the Tau spectral method. The fundamental strategy of the algorithm involves transforming the FODEs, along with its initial (boundary) conditions, into a system of algebraic equations with unknown expansion coefficients. The author thoroughly discussed the convergence and error analysis of the proposed expansion by introducing new inequalities, including the modified Bessel function of the first kind. To validate the proposed methodology, the study in [38] addressed Example 7.6 with a fractional-order parameter of ${\rho }_0=1.5$, considering two distinct source terms as outlined in Example 7.6. We applied our PNM to solve the same problem with an identical fractional-order parameter and the source terms at N = 3 and N = 5, indicating the number of LPs. Specifically, when dealing with $H\left(t\right)=6t+t^3+\frac{5081767996463981t^{\frac{3}{2}}}{1125899906842624}$, our PNM precisely matched the exact solution at N = 3, representing the initial four terms of LPs. In contrast, the results in [38] were obtained at N = 6, representing the initial six terms of Fermat polynomials.

Similarly, for the scenario where $H\left(t\right)=t^5-t^4-12t^2-\frac{4065414397171185t^{\frac{5}{2}}}{562949953421312}+\frac{1451933713275423t^{\frac{7}{2}}}{140737488355328}$, our PNM coincided with the exact solution at N = 5, representing the first six terms of LPs. This aligns with the results in [38], where the approximate solution matched the exact solution at N = 6, corresponding to the initial six terms of Fermat polynomials. This comparison underscores the effectiveness of our proposed approach and the method presented in [38] in achieving accurate solutions with a reduced number of terms in the polynomial expansion. In conclusion, both approaches prove suitable for solving multi-order FODEs. It is recommended to delve into the error analysis and numerical algorithm proposed in [38].

9. Conclusion

Based on the numerical findings achieved in the preceding sections, we can derive the following conclusions.

• This study proposed an efficient numerical method to solve FODEs with both homogeneous and non-homogeneous initial conditions and diverse source terms.

• The construction of the PNM is based on the newly proposed fractional-order derivative operational matrix of LPs in Caputo sense and the fractional-order integral operational matrix of LPs in Riemann–Liouville sense, proposed in [12]. These operational matrices are then used to transform FODEs into Sylvester-type matrix equations which are solved by using the computational software MATLAB for finding the unknown vector, $F^T$.

• Upon comparing the results obtained through the use of the PNM with those presented in Podlubny [18], GA-PS [19], PSO-PS [20], NN-IPA [21], LeNN [22], and He'VIM [23], we observe that the approximate solution obtained by employing the PNM was more similar to the exact solution of the problems, as opposed to the solutions obtained by using other numerical solvers. Our findings indicate that our PNM provides the most precise outcomes compared to other hybrid computational techniques, such as GA-PS, PSO-PS, NN-IPA, and LeNN.

• The effectiveness of the PNM is also highlighted by solving Examples 7.2 and 7.4 for various values of N with PNM and Tau method proposed in [15–17]. We note that our PNM exhibits significantly higher precision and faster convergence toward the exact solution of Examples 7.2 and 7.4. As an illustration, the maximum absolute error computed for Example 7.2 by using the Tau method at N = 3, N = 5, and N = 8 is $1.198\times {10}^{-17}$, $9.94\times {10}^{-18}$, and $9.82\times {10}^{-18}$, respectively, see Table 4. In contrast, the amount of maximum absolute error computed by using PNM is 0, i.e. the exact solution and the approximate solution of Example 7.2 coincide at N = 3, N = 5, and N = 8. Similarly, the maximum absolute error computed for Example 7.4 by using the Tau method at N = 7, N = 10, and N = 12 is $3.28\times {10}^{-10}$, $2.08\times {10}^{-9}$, and $4.18\times {10}^{-9}$, respectively. However, in our case, the amount of maximum absolute error is $4.84\times {10}^{-11}$, see Table 9. Generally, the PNM provides the most precise results as compared to the Tau method for both Examples 7.2 and 7.4.

• To assess the efficiency of the PNM, we explore its capabilities in solving the Bagley-Torvik FODEs with diverse source terms, denoted as $H\left(t\right)$, where computing the analytical solution (64(64) $y\left(t\right)={\int }_0^{t}G(t-z)H\left(z\right)\mathrm{d}z.$(64) ) becomes challenging. We observe that the PNM provides significantly high precise results by solving the Example 7.3 across diverse source terms, see Tables 6, 7, 8, and Figure 3. Additionally, the rational functions, and shifted Heaviside functions are approximated on the half line by truncating the Laguerre series, see Figures 1 and 2.

Moving forward, it would be worthwhile to broaden the scope of the PNM to address the fractional-order partial differential equations and their coupled systems. Additionally, the effectiveness of the hybrid computational techniques previously introduced in [19–21] could be enhanced by approximating the network terms using Laguerre basis and training the optimal weight of networks through the interior point algorithm. This approach could result in the development of a more efficient hybrid computational algorithm that combines the operational matrices of orthogonal polynomials and interior point algorithms.

Our PNM is well-suited for solving problems involving FODEs with initial conditions. However, its effectiveness compromises when applied to FODEs that include boundary conditions. Nevertheless, we can enhance its capabilities by introducing some new operational matrices to approximate the terms of the type, $t{F}_{\mathrm{RL}}^T{J}_{0^{+}}^{\alpha }\psi \left(t\right)$. To illustrate, let's consider the following problem with boundary conditions. (83) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)={}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right),1<\alpha <2,t\in [0,1],\\ & y\left(0\right)=d_0,y\left(1\right)=d_1.\end{array}$(83) Consider the following approximation holds true (84) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)=F^T\psi \left(t\right).$(84) By using the result of Theorem 4.1 and the boundary conditions (83(83) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)={}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right),1<\alpha <2,t\in [0,1],\\ & y\left(0\right)=d_0,y\left(1\right)=d_1.\end{array}$(83) ), we can express Equation (84(84) ${}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)=F^T\psi \left(t\right).$(84) ) as $\begin{array}{rl}y\left(t\right)& \approx F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(t\right)-t{F}^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(1\right)\\ & +d_0+(d_1-d_0)t.\\ & =F^T\left(Q_{(N+1,N+1)}^{\left(\alpha \right)}\psi \right(t)-t{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi (1\left)\right)\\ & +{\mathrm{\Omega }}^T\psi \left(t\right),\end{array}$where, $d_0+(d_1-d_0)t={\mathrm{\Omega }}^T\psi \left(t\right)$. The term, ${}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)$ of Equation (83(83) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)={}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right),1<\alpha <2,t\in [0,1],\\ & y\left(0\right)=d_0,y\left(1\right)=d_1.\end{array}$(83) ) can be expressed as (85) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)& =F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)\\ & -{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}\left(t{F}^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \right(1\left)\right)\\ & +{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right).\end{array}$(85) The primary challenge arises from the fact that the expression in Equation (85(85) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)& =F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)\\ & -{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}\left(t{F}^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \right(1\left)\right)\\ & +{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right).\end{array}$(85) ) does not precisely match the one in Equation (36(36) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right),\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_1}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_1\right)}\psi \left(t\right),\\ & ⋮⋮⋮⋮\\ & {}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_N}y\left(t\right)=F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right)+{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_N\right)}\psi \left(t\right),\\ & \mathrm{and}\\ & H\left(t\right)={\mathrm{{\rm Y}}}^T\psi \left(t\right).\end{array}\}$(36) ) of Section 5. As a result, the problem (83(83) $\begin{array}{rl}& {}_{\mathrm{C}}{\mathcal{D}}^{\alpha }y\left(t\right)={}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)+y\left(t\right),1<\alpha <2,t\in [0,1],\\ & y\left(0\right)=d_0,y\left(1\right)=d_1.\end{array}$(83) ) cannot be simplified into an algebraic system of Sylvester-type equations. To transform them into a system of algebraic equations, it is necessary to approximate the term $t{F}^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \left(1\right)$ in Equation (85(85) $\begin{array}{rl}{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}y\left(t\right)& =F^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right)\\ & -{}_{\mathrm{C}}{\mathcal{D}}^{{\rho }_0}\left(t{F}^T{Q}_{(N+1,N+1)}^{\left(\alpha \right)}\psi \right(1\left)\right)\\ & +{\mathrm{\Omega }}^T{R}_{(N+1,N+1)}^{\left({\rho }_0\right)}\psi \left(t\right).\end{array}$(85) ) by introducing a new operational matrix.

Acknowledgements

We would like to express our sincere gratitude to the anonymous reviewers whose insightful comments and constructive feedback significantly contributed to the refinement and improvement of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254).