Numerical Solutions to Some Families of Fractional Order Differential Equations by Laguerre Polynomials

Abstract

This article is devoted to compute numerical solutions of some classes and families of fractional order differential equations (FODEs). For the required numerical analysis, we utilize Laguerre polynomials and establish some operational matrices regarding to fractional order derivatives and integrals without discretizing the data. Further corresponding to boundary value problems (BVPs), we establish a new operational matrix which is used to compute numerical solutions of boundary value problems (BVPs) of FODEs. Based on these operational matrices (OMs), we convert the proposed (FODEs) or their system to corresponding algebraic equation of Sylvester type or system of Sylvester type. The resulting algebraic equations are solved by MATLAB® using Gauss elimination method for the unknown coefficient matrix. To demonstrate the suggested scheme for numerical solution, many suitable examples are provided.

1. Introduction

The theory of integrals as well as derivatives of arbitrary order is known by the special name “fractional calculus.” It has an old history just like classical calculus. The chronicle of fractional calculus and encyclopedic book can be studied in [1, 2]. Researchers have now necessitated the use of fractional calculus due to its diverse applications in different fields, specially in electrical networks, signal and image processing and optics, etc. For conspicuous work on FODEs in the fields of dynamical systems, electrochemistry, advanced techniques of microorganisms culturing, weather forecasting, as well as statistics, we refer to peruse [3, 4]. Fractional derivatives show valid results in most cases where ordinary derivatives do not. Also annotating that fractional order derivatives as well as fractional integrals are global operators, while ordinary derivatives are local operators. Fractional order derivative provides greater degree of freedom. Therefore from different aspects, the aforesaid areas were investigated. For instance, many researchers have provide understanding to existence and uniqueness results about FODEs, for few results, we refer [5, 6, 7], and many others have actualized the instinctive framework of fractional differential equations in various problems [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19] with many references included in them.

Often it is very difficult to obtain the exact solution due to global nature of fractional derivatives in differential equations. Contrarily approximate solutions are obtained by numerical methods assorted in [20, 21, 22]. Various new numerical methods have been developed, among them is one famous method called “spectral method” which is used to solve problems in various realms [23]. In this method operational matrices are obtained by using orthogonal polynomials [24]. Many authors have successfully developed operational matrices by using Legender, Jacobi, and various other polynomials [25, 26]. For delay differential and various other related equations, Laguerre spectral methods have been used [27, 28, 29, 30, 31, 32]. Bernstein polynomials and various classes of other polynomials were also used to obtain operational matrices corresponding to fractional integrals and derivatives [33, 34, 35, 36, 37, 38, 39, 40]. Apart from them, operational matrices were also developed with the collocation method (see Refs. [41, 42, 43]). Since spectral methods are powerful tools to compute numerical solutions of both ODEs and FODEs. Therefore, we bring out numerical analysis via using Laguerre polynomials of some families and coupled systems of FODEs under initial as well as boundary conditions. In this regard we investigate the numerical solutions to the given families under initial conditions

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)\pm z\left(t\right)=0,0<\gamma \le 1,& \\ z\left(0\right)=z_0,z_0\in R,& \end{array}\right.$$E1

and subject to boundary conditions

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)\pm z\left(t\right)=0,1<\gamma \le 2,& \\ z\left(0\right)=z_0,z\left(1\right)=z_1,z_0,z_1\in R.& \end{array}\right.$$E2

By similar numerical techniques, we also investigate the numerical solutions to the following systems with fractional order derivatives under initial and boundary conditions as

$$\begin{array}{c}\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)+\mathit{az}\left(t\right)+\mathit{by}\left(t\right)=f\left(t\right),& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)+\mathit{cy}\left(t\right)+\mathit{dz}\left(t\right)=g\left(t\right),& \end{array}\right.\\ z\left(0\right)=z_0,y\left(0\right)=y_0\end{array}$$E3

for $0<\gamma \le 1$ and

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)\mathit{az}\left(t\right)+\mathit{by}\left(t\right)=f\left(t\right),& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)+\mathit{cy}\left(t\right)+\mathit{dz}\left(t\right)=g\left(t\right),& \\ z\left(0\right)=z_0,y\left(0\right)=y_0,z\left(1\right)=z_1,y\left(1\right)=y_1,& \end{array}\right.$$E4

for $1<\gamma \le 2$ where $f,g:\left[0,1\right]\times R^2\to R$ and $z_0,y_0,z_1,y_1\in R.$ We first obtain OMs for fractional derivatives and integrals by using Laguerre polynomials. Also corresponding to boundary conditions, we construct an operational matrix which is needed in numerical analysis of BVPs. With the help of the OMs we convert the considered problem of FODEs under initial/boundary conditions to Sylvester-type algebraic equations. Solving the mentioned matrix equations by using MATLAB®, we compute the numerical solutions of the considered problems.

2. Preliminaries

Here we recall some basic definition results that are needed in this work onward, keeping in mind that throughout the paper we use fractional derivative in Caputo sense.

Definition 1. The fractional integral of order $\gamma >0$ of a function $z:\left(0,\infty \right)\to R$ is defined by

$${}_0{I}_t^{\gamma }z\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\gamma \right)}{\int }_0^t\frac{z\left(s\right)}{{\left(t-s\right)}^{1-\gamma }}\mathit{ds},$$

provided the integral converges at the right sides. Further a simple and important property of ${}_0{I}_t^{\gamma }$ is given by

$${}_0{I}_t^{\gamma }{t}^{\delta }=\frac{\mathrm{\Gamma }\left(\delta +1\right)}{\mathrm{\Gamma }\left(\delta +\gamma +1\right)}{t}^{\gamma +\delta }.$$

Definition 2. Caputo fractional derivative is defined as

$${}_0^c{D}_t^{\gamma }f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(n-\gamma \right)}{\int }_0^t{\left(t-s\right)}^{n-\gamma -1}{f}^{\left(n\right)}\left(s\right)\mathit{ds},$$

where $n$ is a positive integer with the property that $n-1<\gamma \le n.$ For example, if $0<\gamma \le 1$, then Caputo fractional derivative becomes

$${}_0^c{D}_t^{\gamma }f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(1-\gamma \right)}{\int }_0^t{\left(t-s\right)}^{-\gamma -1}{f}^{\text{'}}\left(s\right)\mathit{ds}.$$

Theorem 1. The FODE given by

$${}_0^c{D}_t^{\gamma }f\left(t\right)=0$$

has a unique solution, such that

$$f\left(t\right)=d_0+d_{1}t+d_2{t}^2+\dots +d_{n-1}{t}^{n-1},n=\left[\gamma \right]+1.$$

Lemma 1. Therefore in view of this result, if $h\in L^n\left[0,T\right],$ then the unique solution of nonhomogenous FODE

$${}_0^c{D}_t^{\gamma }f\left(t\right)=h\left(t\right),n-1<\gamma \le n$$

is written as

$$f\left(t\right)=d_0+d_{1}t+d_2{t}^2+\dots +d_{n-1}{t}^{n-1}{+}_0{I}_t^{\gamma }h\left(t\right),$$

where $d_i$ for $i=\mathrm{0,1,2,3}\dots n-1$ are real constants.

The above lemma is also stated as

$$f\left(t\right){=}_0{I}_t^{\gamma }h\left(t\right)+\sum _{i=0}^{n-1}\frac{f^i\left(0\right)}{i!}{t}^i.$$

Definition 3. The famous Laguerre polynomials are represented by $L_i^{\gamma }\left(t\right)$ and defined as

$$L_i^{\gamma }\left(t\right)=\sum _{k=0}^i\frac{{\left(-1\right)}^k\mathrm{\Gamma }\left(i+\gamma +1\right)}{\mathrm{\Gamma }\left(k+1+\gamma \right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+1\right)}{t}^k.$$

They are orthogonal on $\left[0,\infty \right].$ If $L_i^{\gamma }\left(t\right)$ and $L_j^{\gamma }\left(t\right)$ are Laguerre polynomials, then the orthogonality condition is given as

$${\int }_0^{\infty }{L}_i^{\gamma }\left(t\right)L_j^{\gamma }\left(t\right)W^{\gamma }\left(t\right)\mathit{dt}={\delta }_{i,j}{U}_k,$$

where

$$W^{\gamma }\left(t\right)=t^{\gamma }{e}^{-t},$$

is the weight function and

$$U_k=\left\{\begin{array}{ll}\frac{\mathrm{\Gamma }\left(1+\gamma +k\right)}{\mathrm{\Gamma }\left(1+k\right)},& i=j\\ 0& i\ne j.\end{array}\right.$$

Now let $Z\left(t\right)$ be any function, defined on the interval $\left[0,\infty \right].$ We express the function in terms of Laguerre polynomials as

$$\begin{array}{c}Z\left(t\right)=\sum _{i=0}^n{c}_i{L}_i^{\gamma }\left(t\right).\\ =c_0{L}_0^{\gamma }\left(t\right)+c_1{L}_1^{\gamma }\left(t\right)+\dots +c_N{L}_N^{\gamma }\left(t\right)\\ =\left[c_0{c}_1\dots c_N\right]\left[\begin{array}{c}{L}_0^{\gamma }\left(t\right)\\ ⋮\\ L_n^{\gamma }\left(t\right)\\ \end{array}\right].\end{array}$$E5

We set the above two vectors into their inner product and represent the column matrix by $\mathrm{\Psi }\left(t\right),$ so that

$$Z\left(t\right)=c^t\mathrm{\Psi }\left(t\right).$$

Again as

$$\begin{array}{c}Z\left(t\right)=\sum _{i=0}^n{c}_i{L}_i^{\gamma }\left(t\right),\\ {\int }_0^{L}Z\left(t\right)W^{\gamma }\left(t\right)L_j^{\gamma }\left(t\right)\mathit{dt}={\int }_0^L\sum _{i=0}^n{c}_i{L}_i^{\gamma }\left(t\right)L_j^{\gamma }\left(t\right)W^{\gamma }\left(t\right)\mathit{dt},\end{array}$$

which is written as

$$\sum _{i=0}^n{c}_i{\int }_0^L{L}_i^{\gamma }\left(t\right)L_j^{\gamma }\left(t\right)W^{\gamma }\left(t\right)\mathit{dt}.$$

We call $h_i$ to the general term of integration

$${\int }_0^{L}Z\left(t\right)W^{\gamma }\left(t\right)L_j^{\gamma }\left(t\right)\mathit{dt}=\sum _{i=0}^n{c}_i{h}_i.$$

Hence the coefficient $c_i$ is

$$c_i=\frac{1}{h_i}{\int }_0^{L}Z\left(t\right)W^{\gamma }\left(t\right)L_j^{\gamma }\left(t\right)\mathit{dt}.$$

In vector form we can write Eq. (5) as

$$Z\left(t\right)=c_M^t{\mathrm{\Psi }}_M\left(t\right).$$

where $M$ = $m+1,$$c_M$ is the $M$ terms coefficient vector and ${\mathrm{\Psi }}_M\left(t\right)$ is the $M$ terms function vector.

2.1 Representation of Laguerre polynomial with Caputo fractional order derivative

If the Caputo fractional order derivative is applied to Laguerre polynomial, by considering whole function constant except $t^k.$ We use the definition of Caputo fractional order derivative for $t^k$ to obtain (6) as

$${}_0^c{D}_t^{\gamma }{L}_i^{\gamma }\left(t\right)=\sum _{k=0}^i\left(t^{k-\gamma }\right)\frac{{\left(-1\right)}^k\mathrm{\Gamma }\left(i+\gamma +1\right)}{\mathrm{\Gamma }\left(k+1+\gamma \right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(1+k-\gamma \right)}.$$E6

2.2 Error analysis

The proof of the following results can be found with details in [20].

Lemma 2. Let $L_i^{\beta }\left(t\right)$ be given; then

$${}_0^c{D}_t^{\gamma }{L}_i^{\beta }\left(t\right)=0,i=\mathrm{0,1,2},\cdots ,\left[\beta \right]-1,\gamma >0.$$

Theorem 2. For error analysis, we state the theorem such that, $a$ be any integer and $0\le s\le a$, and then

$$\parallel P_{M,\mathit{az}}-z\left(t\right)\parallel A_{\alpha }^s,\mathrm{\Lambda }\le {\mathit{cM}}^{\frac{s-a}{2}}\mid z\left(t\right)\mid A_{\alpha }^a,\mathrm{\Lambda },\forall z\left(t\right)\mathrm{ϵ}{A}_{\alpha }^a\left(\mathrm{\Lambda }\right),$$

where $A_{\alpha }^a$$=$$\{$$z/z$ is measurable on $\mathrm{\Lambda }$ and $\parallel z\parallel$$A_{\alpha }^a,\left(\mathrm{\Lambda }\right)$$<$$\infty$$\}$ and

$$\begin{array}{c}\mid z\mid A_{\alpha }^a,\left(\mathrm{\Lambda }\right)=\parallel {\partial }_p^{a}z{\parallel }_{\mathit{w\alpha }+a,\mathrm{\Lambda }},\\ \parallel z\parallel A_{\alpha }^a,\left(\mathrm{\Lambda }\right)={\left(\sum _{k=0}^a{\left|z\right|}_{A_{\alpha }^a,\left(\mathrm{\Lambda }\right)}^2\right)}^{\frac{1}{2}}.\end{array}$$

Now let $\mathrm{\Lambda }=\mathrm{\varrho }/0<\mathrm{\varrho }<\infty$ with $\chi \left(\mathrm{\varrho }\right)$ be a weight function. Then

$L_{\chi }^2\left(\mathrm{\Lambda }\right)$$=$$\{\kappa$$/$$\kappa$ is measurable on $\mathrm{\Lambda }$ and $\parallel u{\parallel }_{L_{\chi }^2},\mathrm{\Lambda }$$<$$\infty \}$.

with the following inner product and norm

$${\left(u,v\right)}_{\chi ,\mathrm{\Lambda }}={\int }_{\mathrm{\Lambda }}u\left(\mathrm{\varrho }\right)v\left(\mathrm{\varrho }\right)d\mathrm{\varrho },\parallel v\parallel \chi ,\mathrm{\Lambda }={\sqrt{⟨u,v⟩}}_{\chi ,\mathrm{\Lambda }}.$$

3. Operational matrices corresponding to fractional derivatives and integrals

Here in this section, we provide the required OMs via Laguerre polynomials of fractional derivatives and integrals.

Lemma 3. Let ${\mathrm{\Psi }}_M\left(t\right)$ be a function vector; the fractional integral of order $\gamma$ for the function ${\mathrm{\Psi }}_M\left(t\right)$ can be generalized as

$${}_0{I}_t^{\gamma }{\mathrm{\Psi }}_M\left(t\right)\approx G_{N\times N}^{\gamma }{\mathrm{\Psi }}_M\left(t\right),$$

where $G_{N\times N}^{\gamma }$ is the OM of integration of fractional order $\gamma$ and given by

$$\left[\begin{array}{cccccc}{\mathcal{\daleth }}_{0,0,k,r}^{\gamma }& {\mathcal{\daleth }}_{0,1,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{0,j,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{0,m,k,r}^{\gamma }\\ {\mathcal{\daleth }}_{1,0,k,r}^{\gamma }& {\mathcal{\daleth }}_{1,i,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{1,j,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{1,m,k,r}^{\gamma }\\ ⋮& ⋮& ⋮& ⋮& \cdots & ⋮\\ {\mathcal{\daleth }}_{i,0,k,r}^{\gamma }& {\mathcal{\daleth }}_{i,1,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{i,j,k,r}^{\gamma }& ⋮& {\mathcal{\daleth }}_{i,m,k,r}^{\gamma }\\ ⋮& ⋮& ⋮& ⋮& \cdots & ⋮\\ {\mathcal{\daleth }}_{m,0,k,r}^{\gamma }& {\mathcal{\daleth }}_{m,1,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{m,j,k,r}^{\gamma }& \cdots & {\mathcal{\daleth }}_{m,m,k,r}^{\gamma }\end{array}\right],$$

where

$$\begin{array}{ll}& {\mathcal{\daleth }}_{i,j,k,r}^{\gamma }=\sum _{k=0}^i\sum _{r=0}^i\frac{{\left(-1\right)}^{k+r}\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(i+\gamma +1\right)\mathrm{\Gamma }\left(k+\gamma +\alpha +r+1\right)}{\mathrm{\Gamma }\left(j-r+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(r+1\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(k+\alpha +1\right)\mathrm{\Gamma }\left(\gamma +r+1\right)}.\end{array}$$

Proof. We apply the fractional order integral of order $\gamma$ to the Laguerre polynomials

$${}_0^c{I}_t^{\gamma }{L}_i^{\gamma }\left(t\right)=\sum _{k=0}^i{\frac{\mathrm{\Gamma }\left(i+\gamma +1\right)}{\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(k+1\right)}}_0^c{I}_t^{\gamma }{t}^k.$$E7

Since from (7), we have

$${}_0^c{I}_t^{\gamma }{t}^k=\frac{\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(1+k+\alpha \right)}{t}^{k+\gamma }.$$

Therefore Eq. (7) implies that

$${}_0^c{I}_t^{\gamma }{L}_i^{\gamma }\left(t\right)=\sum _{k=0}^i{t}^{k+\gamma }\frac{\mathrm{\Gamma }\left(i+\gamma +1\right)}{\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(k+1\right)}\frac{\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(1+k+\alpha \right)},$$

which is equal to

$${}_0^c{I}_t^{\gamma }{L}_i^{\gamma }\left(t\right)=\sum _{k=0}^i{\left(-1\right)}^k\frac{\mathrm{\Gamma }\left(i+\gamma +1\right)}{\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(1+k-\gamma \right)}{t}^{k+\gamma }.$$E8

We approximate $t^{k+\gamma }$ in (8) with Laguerre polynomials, i.e.

$$t^{k+\gamma }\approx \sum _{j=0}^n{H}_j{L}_j^{\gamma }\left(t\right).$$

By using the relation of orthogonality, we can find coefficients

$$H_j=\sum _{r=0}^j{\left(-1\right)}^k\frac{\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(k+\alpha +r+\gamma +1\right)}{\mathrm{\Gamma }\left(1+j-r\right)\mathrm{\Gamma }\left(1+r\right)\mathrm{\Gamma }\left(1+r+\gamma \right)}.$$

So Eq. (8) implies

$$\begin{array}{ll}& {}_0^c{I}_t^{\gamma }{L}_i^{\gamma }\left(t\right)=\sum _{k=0}^i{\left(-1\right)}^k\frac{\mathrm{\Gamma }\left(i+\gamma +1\right)}{\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(1+k-\gamma \right)}\\ & \times \sum _{r=0}^j{\left(-1\right)}^r\frac{\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(k+\alpha +r+\gamma +1\right)}{\mathrm{\Gamma }\left(j-r+1\right)\mathrm{\Gamma }\left(r+1\right)\mathrm{\Gamma }\left(r+\gamma +1\right).}\end{array}$$

$$\begin{array}{ll}& {}_0^c{I}_t^{\gamma }{L}_i^{\gamma }\left(t\right)=\\ & \sum _{k=0}^i\sum _{r=0}^j{\left(-1\right)}^{k+r}\frac{\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(i+\gamma +1\right)\mathrm{\Gamma }\left(k+\alpha +r+\gamma +1\right)}{\mathrm{\Gamma }\left(1-k+i\right)\mathrm{\Gamma }\left(j-\gamma +1\right)\mathrm{\Gamma }\left(\gamma +1\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(k+\alpha +1\right)\mathrm{\Gamma }\left(\gamma +r+1\right)}.\end{array}$$

which is the desired result.

Lemma 4. Let ${\mathrm{\Psi }}_M\left(t\right)$ be a function vector; then the fractional derivative of order $\gamma$ for ${\mathrm{\Psi }}_M\left(t\right)$ is generalized as

$${}_0^c{D}_t^{\gamma }{\mathrm{\Psi }}_M\left(t\right)\approx {\mathbf{W}}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M\left(t\right),$$

where ${\mathbf{W}}_{M\times M}^{\gamma }$ is the OM of derivative of order $\gamma$, defined as in (9)

$$\begin{array}{c}{\mathbf{W}}_{M\times M}^{\gamma }=\left[\begin{array}{cccccc}0& 0& 0& 0& \cdots 0& \\ ⋮& ⋮& ⋮& ⋮& ⋮& \\ {\mathrm{\Theta }}_{⌈\gamma ⌉,0,k,\alpha }^{\gamma }& {\mathrm{\Theta }}_{⌈\gamma ⌉,1,k,\alpha }^{\gamma }& \cdots {\mathrm{\Theta }}_{⌈\gamma ⌉,j,k,\alpha }^{\gamma }& \cdots & \cdots & {\mathrm{\Theta }}_{⌈\gamma ⌉,n,k,\alpha }^{\gamma }\\ {\mathrm{\Theta }}_{i,0,k,\alpha }^{\gamma }& {\mathrm{\Theta }}_{i,1,k,\alpha }^{\gamma }& {\mathrm{\Theta }}_{i,j,k,\alpha }^{\gamma }& \cdots & \cdots & {\mathrm{\Theta }}_{i,n,k,\alpha }^{\gamma }\\ ⋮& ⋮& ⋮& ⋮& ⋮& \\ {\mathrm{\Theta }}_{n,0,k,\alpha }^{\gamma }& {\mathrm{\Theta }}_{n,1,k,\alpha }^{\gamma }& {\mathrm{\Theta }}_{n,j,k,\alpha }^{\gamma }& \cdots & \cdots & {\mathrm{\Theta }}_{n,n,k,\alpha }^{\gamma }\end{array}\right],\end{array}$$E9

where

$$\begin{array}{ll}& {\mathrm{\Theta }}_{i,j,k,\alpha }^{\gamma }=\\ & \sum _{k=\gamma }^i\sum _{r=0}^i\frac{{\left(-1\right)}^{\gamma +k}\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(i+\alpha +1\right)\mathrm{\Gamma }\left(k+\alpha -r+\gamma +1\right)}{\mathrm{\Gamma }\left(j-r+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(r+1\right)\mathrm{\Gamma }\left(k+\alpha +1\right)\mathrm{\Gamma }\left(k-\gamma +1\right)\mathrm{\Gamma }\left(\alpha +\gamma +1\right)}.\end{array}$$

Proof. Leaving the proof as it is very similar to the proof of the above lemma.

Lemma 5. We consider a function $Z\left(t\right)$ defined on $\left[0,\infty \right]$ and $y\left(t\right)=K_M{\mathrm{\Psi }}_M^T\left(t\right)$; then

$$Z\left(t\right)\left[{}_0{I}_t^{\gamma }y\left(t\right)\right]=K_M{Q}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M\left(t\right),$$

where $Q_{M\times M}^{\gamma }$ is the operational matrix, given by

$$\left[\begin{array}{cccccc}{C}_{0,0,}& C_{0,1}& \cdots & C_{0,j}& \cdots & C_{0,m}\\ C_{1,0}& C_{1,1}& \cdots & C_{1,j}& \cdots & C_{1,m}\\ ⋮& ⋮& ⋮& ⋮& \cdots & ⋮\\ C_{i,0}& C_{i,1}& \cdots & C_{i,j}& ⋮& C_{i,m}\\ ⋮& ⋮& ⋮& ⋮& \cdots & ⋮\\ C_{m,0}& C_{m,1}& \cdots & C_{m,j}& \cdots & C_{m,m}\end{array}\right],$$

where

$$C_{i,j}=\frac{1}{h_i}{\int }_0^1{\mathrm{\Delta }}_{i,\gamma ,k}Z\left(t\right)L_j^{\gamma }\left(t\right)\mathit{dt},$$

with

$$w_i=\sum _{k=0}^i\frac{{\left(-1\right)}^{i+1}\mathrm{\Gamma }\left(i+1+\gamma \right)}{\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(1-k+i\right)\mathrm{\Gamma }\left(k+\gamma \right)}.$$

Proof. By considering the general term of ${\mathrm{\Psi }}_M\left(t\right)$

$$\begin{array}{c}{}_0{I}_1^{\gamma }{L}_i\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\gamma \right)}{\int }_0^1{\left(1-s\right)}^{\gamma -1}{L}_i\left(s\right)\mathit{ds}.\\ {}_0{I}_1^{\gamma }{L}_i\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\gamma \right)}{\int }_0^1{\left(1-s\right)}^{\gamma -1}\sum _{k=0}^i{\left(s\right)}^k\frac{{\left(-1\right)}^k\mathrm{\Gamma }\left(i+1+\gamma \right)}{\mathrm{\Gamma }\left(-k+1+i\right)\mathrm{\Gamma }\left(k+1+\gamma \right)\mathrm{\Gamma }\left(1+k\right)}\mathit{ds}.\\ {}_0{I}_1^{\gamma }{L}_i\left(t\right)=\sum _{k=0}^i\frac{{\left(-1\right)}^k\mathrm{\Gamma }\left(i+1+\gamma \right)}{\mathrm{\Gamma }\left(\gamma \right)\mathrm{\Gamma }\left(-k+1+i\right)\mathrm{\Gamma }\left(k+1+\gamma \right)\mathrm{\Gamma }\left(1+k\right)}{\int }_0^1{\left(1-s\right)}^{\gamma -1}{\left(s\right)}^k\mathit{ds}.\end{array}$$E10

Using the famous Laplace transform, we have from (10)

$$\begin{array}{c}\mathit{\text{£}}({\int }_0^1{\left(1-s\right)}^{\gamma -1}{s}^k\mathit{ds}=\frac{\mathrm{\Gamma }\left(\gamma \right)\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(\gamma +k\right)}.\\ {}_0{I}_1^{\gamma }{L}_i\left(t\right)=\sum _{k=0}^i\frac{{\left(-1\right)}^k\mathrm{\Gamma }\left(i+1+\gamma \right)}{\mathrm{\Gamma }\left(\gamma \right)\mathrm{\Gamma }\left(-k+1+i\right)\mathrm{\Gamma }\left(k+1+\gamma \right)\mathrm{\Gamma }\left(1+k\right)}\frac{\mathrm{\Gamma }\left(\gamma \right)\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(\gamma +k\right)}.\\ \sum _{k=0}^i\frac{{\left(-1\right)}^k\mathrm{\Gamma }\left(i+1+\gamma \right)}{\mathrm{\Gamma }\left(-k+1+i\right)\mathrm{\Gamma }\left(k+\gamma +1\right)\mathrm{\Gamma }\left(1+k\right)}={\mathrm{\Delta }}_{i,\gamma ,k}.\end{array}$$

Now using Laguerre polynomials, we have

$${\mathrm{\Delta }}_{i,\gamma ,k}z\left(t\right)=\sum _{j=0}^m{C}_{i,j}{L}_i\left(t\right),$$

where $C_{i,j}$ is calculated by using orthogonality as

$$C_{i,j}=\frac{1}{\mathit{hi}}{\int }_0^1{\mathrm{\Delta }}_{i,\gamma ,k}z\left(t\right)L_j^{\gamma }\left(t\right)\mathit{dt}.$$E11

To get the desired result, we evaluate the above (11) relation for $i=0,1,\dots ,m$ and $j=0,1,\dots ,m$.

4. Main result

In this section, we discuss some cases of FODEs with initial condition as well as boundary conditions. The approximate solution obtained through desired method is compared with the exact solution. Similarly we investigate numerical solutions to various coupled systems under some initial conditions as well as boundary conditions.

4.1 Treatment of FODEs under initial and boundary conditions

Here we discuss different cases.

Case 1. In the first case, we consider the fractional order differential equation

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)\pm z\left(t\right)=0,0<\gamma ⩽1,& \\ z\left(0\right)=z_0,z_0\in R& \end{array}\right.$$E12

we see that

$${}_0^c{D}_t^{\gamma }z\left(t\right)={\mathrm{Ł}}_M{\psi }_M^T\left(t\right).$$

and applying ${}_0{I}_t^{\gamma }$ by the Lemma 1, on (12) we write

$$z\left(t\right)=e_0+{}_0{I}_t^{\gamma }\left[{\mathrm{Ł}}_M{\psi }_M^T\left(t\right)\right],$$

Using the initial condition to get $e_0=z_0$ and approximate $z_0$ as $z_0\approx F_M{\psi }_M^T\left(t\right)$, Eq. (12) implies

$${\mathrm{Ł}}_M{\psi }_M^T\left(t\right)+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)+F_M{\psi }_M^T\left(t\right)=0.$$

Finally the Sylvester-type algebraic equation is obtained as

$${\mathrm{Ł}}_M+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)+F_M=0.$$

Solving the Sylvester matrix for ${\mathrm{Ł}}_M,$ we get the numerical value for $z\left(t\right)$.

Example 1.

$$\left\{\begin{array}{l}{}_0^c{D}_t^{\gamma }z\left(t\right)\pm z\left(t\right)=0,0<\gamma \le 1,\\ z\left(0\right)=1,z_0\in R.\end{array}\right.$$

Since the exact solution is given by

$$z\left(t\right)=E_{\gamma }\left(-t^{\gamma }\right),$$

where $E_{\gamma }$ is the Mittag-Leffler representation, and at $\gamma =1,$$z\left(t\right)$$=e^{-t}.$

Approximating the solution through the proposed method and plotting the exact as well as numerical solution by using scale $M=8$ corresponding to $\gamma =1$ in Figure 1, we see that the proposed method works very well.

Figure 1.

Plots of both approximate and exact solution for the Example 1 for Case 1.

Case 2.

$$\left\{\begin{array}{ll}& {}_0^c{D}_t^{\gamma }z\left(t\right)+z\left(t\right)=0,1<\gamma ⩽2,\\ & z\left(0\right)=z_0,z\left(1\right)=z_1,z_0,z_1\in R.\end{array}\right.$$E13

We take

$${}_0^c{D}_t^{\gamma }z\left(t\right)=K_M{\psi }_M^T\left(t\right).$$E14

Applying Lemma 1 to Eq. (14), we get

$$z\left(t\right)=e_0+e_1\left(t\right){+}_0{I}_t^{\gamma }{K}_M{\psi }_M^T\left(t\right).$$E15

Using the conditions by putting $t=0$ and $t=1$ to get $e_0=z_0$ and

$$e_1=z_1-z_0-K_{M0}{I}_1^{\gamma }{\psi }_M^T\left(t\right){/}_{t=1}.$$

Equation (15) implies

$$z\left(t\right)=z_0+\left(z_1-z_0\right)t-{\mathit{tK}}_{M0}{I}_1^{\gamma }{\psi }_M^T\left(t\right){/}_{t=1}{+}_0{I}_t^{\gamma }{K}_M{\psi }_M^T\left(t\right),$$

where $z_0+\left(z_1-z_0\right)t$ is the smooth function of $t$ and constants; we approximate it as

$$z_0+\left(z_1-z_0\right)t\approx G_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)$$

and

$${\mathit{tK}}_{M0}{I}_1^{\gamma }{\psi }_M^T\left(1\right)\approx K_M{Q}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right).$$

Hence

$$z\left(t\right)=G_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)-K_M{Q}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)$$

So Eq. (13) implies

$$K_M{\psi }_M^T\left(t\right)+G_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)-K_M{Q}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)=0$$

which is further solved for $K_M$ to get the required numerical solution.

For Case 2, we give the following example.

Example 2.

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)+z\left(t\right)=0,0<\gamma \le 2,& \\ z\left(0\right)=-1,z\left(1\right)=1.& \end{array}\right.$$E16

At $\gamma =2,$ we get the exact solution as of (16) as given by (17)

$$z\left(t\right)=114.58sin\left(x\right)-cos\left(x\right)$$E17

Upon using the suggested method, we see from the subplot at the left of Figure 2 that exact and numerical solutions are very close to each other for very low scale level. Also, the absolute error is given in subplot at the right of Figure 2.

Figure 2.

The plot of exact and approximate solution for Example 2 for Case 2.

4.2 Coupled systems of linear FODEs under initial and boundary conditions

In this subsection, we consider different forms of coupled systems of FODEs with the initials as well as boundary conditions.

Case 1. First we take the coupled system of FODEs as

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)+\mathit{az}\left(t\right)+\mathit{by}\left(t\right)=f\left(t\right)& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)+\mathit{cy}\left(t\right)+\mathit{dz}\left(t\right)=g\left(t\right),& \end{array}\right.$$E18

with the conditions

$$z\left(0\right)=z_0,y\left(0\right)=y_0,z_0,y_0\in R.$$E19

Let

$${}_0^c{D}_t^{\gamma }z\left(t\right)={\mathrm{Ł}}_M{\psi }_M^T\left(t\right),{}_0^c{D}_t^{\gamma }y\left(t\right)=K_M{\psi }_M^T\left(t\right).$$E20

Applying Lemma 1 to Eq. (20), we get

$$\left\{\begin{array}{ll}z\left(t\right)=e_0+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right),& \\ y\left(t\right)=d_0+K_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right).& \end{array}\right.$$E21

Using the initial conditions given in Eq. (19), from Eq. (21), we get

$$\left\{\begin{array}{ll}z\left(t\right)=F_M^1{\psi }_M^T\left(t\right)+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right),& \\ y\left(t\right)=y_0\approx F_M^2{\psi }_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right).& \end{array}\right.$$E22

We take approximation as

$$z_0\approx F_M^1{\psi }_M^T\left(t\right),$$

and

$$y_0\approx F_M^2{\psi }_M^T\left(t\right),$$

while source functions are approximated as

$$f\left(t\right)\approx F_M^3{\mathrm{\Psi }}_M^T\left(t\right),$$

and

$$g\left(t\right)\approx F_M^4{\mathrm{\Psi }}_M^T\left(t\right).$$

Therefore the consider system on using (19)–(22), (18) becomes

$$\begin{array}{l}\left\{\begin{array}{ll}{\mathrm{Ł}}_M{\psi }_M^T+a\left(F_M^1{\psi }_M^T\left(t\right)+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)\right)& \\ +b(F_M^2{\psi }_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)=F_M^3{\psi }_M^T\left(t\right).& \end{array}\right.\\ \left\{\begin{array}{ll}{K}_M{\psi }_M^T+c\left(F_M^2{\psi }_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)\right)& \\ +d(F_M^1{\psi }_M^T\left(t\right)+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\psi }_M^T\left(t\right)=F_M^4{\psi }_M^T\left(t\right).& \end{array}\right.\end{array}$$

On further rearrangement we have

$$\left\{\begin{array}{ll}{\mathrm{Ł}}_M+a\left(F_M^1+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }\right)+b(F_M^2+K_M{G}_{M\times M}^{\gamma }=F_M^3& \\ K_M+c\left(F_M^2+K_M{G}_{M\times M}^{\gamma }\right)+d(F_M^1+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }=F_M^4.& \end{array}\right.$$

which further can be written as

$$\left\{\begin{array}{ll}{\mathrm{Ł}}_M\left(I_{M\times M}+{\mathit{aG}}_{M\times M}^{\gamma }\right)+K_M\left({\mathit{bG}}_{M\times M}^{\gamma }\right)+\left({\mathit{aF}}_M^1+{\mathit{bF}}_M^2-F_M^3\right)=0& \\ K_M\left(I_{M\times M}+{\mathit{cG}}_{M\times M}^{\gamma }\right)+{\mathrm{Ł}}_M\left({\mathit{dG}}_{M\times M}^{\gamma }\right)+\left({\mathit{cF}}_M^2+{\mathit{dF}}_M^1-F_M^4\right)=0.& \end{array}\right.$$

In matrix form we write as

$$\begin{array}{c}\left[{\mathrm{Ł}}_M{K}_M\right]\left[\begin{array}{cc}{I}_{M\times M}+{\mathit{aG}}_{M\times M}^{\gamma }& 0\\ 0& I_{M\times M}+{\mathit{cG}}_{M\times M}^{\gamma }\end{array}\right]+\left[\begin{array}{cc}{\mathrm{Ł}}_M& K_M\\ & \end{array}\right]\left[\begin{array}{cc}0& {\mathit{dG}}_{M\times M}^{\gamma }\\ {\mathit{bG}}_{M\times M}^{\gamma }& 0\end{array}\right]\\ +\left[\begin{array}{c}{\mathit{aF}}_M^1+{\mathit{bF}}_M^2-F_M^3\\ {\mathit{cF}}_M^2+{\mathit{dF}}_M^1-F_M^4\end{array}\right]=0.\end{array}$$

We solve this system of matrix equation for $\left[{\mathrm{Ł}}_M{K}_M\right]$ by using Gaussian’s elimination method. The considered system is in the form of $X\overline{A}+X\overline{B}+\overline{C}=0,$.

where $X=\left[{\mathrm{Ł}}_M{K}_M\right]$$\overline{A}=\left[\begin{array}{cc}{I}_{M\times M}+{\mathit{aG}}_{M\times M}^{\gamma }& 0\\ 0& I_{M\times M}+{\mathit{cG}}_{M\times M}^{\gamma }\end{array}\right],$.

$\overline{B}=\left[\begin{array}{cc}0& {\mathit{dG}}_{M\times M}^{\gamma }\\ {\mathit{bG}}_{M\times M}^{\gamma }& 0\end{array}\right]$ and $\overline{C}=\left[\begin{array}{c}{\mathit{aF}}_M^1+{\mathit{bF}}_M^2-F_M^3\\ {\mathit{cF}}_M^2+{\mathit{dF}}_M^1-F_M^4.\end{array}\right].$

Upon computation of matrices ${\mathrm{Ł}}_M,K_M$ by using MATLAB®, we put these matrices in Eq. (22) to find $z_{\mathit{app}}$ and $y_{\mathit{app}}$, respectively.

Example 3. We now provide its example by considering the system of FODEs:

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)+z\left(t\right)+y\left(t\right)=f\left(t\right)& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)+y\left(t\right)+z\left(t\right)=g\left(t\right),& \\ z\left(0\right)=2,y\left(0\right)=1.& \end{array}\right.$$

By taking $\gamma =1,$ the exact solution is obtained as

$$z\left(t\right)=cos\left(t\right)+e^t,y=sin\left(t\right)+e^{-t},$$

where the external source functions are given by $f\left(t\right)=cos\left(t\right)+e^{-t}+2e^t$ and $g\left(t\right)=$$e^{-t}+sin\left(t\right)+2cos\left(t\right).$ The exact solution $z_{\mathit{ex}},y_{\mathit{ex}}$ can be computed by any method of ODEs. Approximating the problem by the considered method, we see that the computed numerical and exact solutions have close agreement at very small-scale level. The corresponding accuracy has been recorded in Table 1. Further the comparison between exact and numerical solution and the results about absolute error have been demonstrated in Figures 3 and 4, respectively. In Figure 3 we are given the comparison between exact solution and approximate solutions by using proposed method. Similarly the absolute errors have been described in Figure 4.

$t$	CPU time (s)	Absolute error $\parallel z_{\mathit{app}}-z_{\mathit{ex}}\parallel$	Absolute error $\parallel y_{\mathit{app}}-y_{\mathit{ex}}\parallel$	CPU time (s)
0	30.5	0.00003	0.000006	32.5
0.15	32.7	0.000016	0.000034	33.3
o.35	35.8	0.000013	0.00003	33.9
0.65	33.6	0.000012	0.00003	35.6
0.87	34.8	0.000018	0.000036	36.5
1	35.9	0.00003	0.000006	36.8

Table 1.

Absolute error at $M=5,\gamma =0.9$, for different values of $t$ in Example 3.

Figure 3.

Plots of exact and approximate solution of Example 3.

Figure 4.

Plots of absolute error of Example 3.

By comparing the exact and numerical solution through the proposed method, we observe that our numerical solution does not show any disagreement with the exact solution as can be seen in Figure 3. The absolute errors $\parallel z_{\mathit{app}}-z_{\mathit{ex}}\parallel$ and $\parallel y_{\mathit{app}}-y_{\mathit{ex}}\parallel$ plotted at the scale $M=5$ are very low as given in Figure 4, which describes the efficiency of the proposed method.

Case 2. Similarly for the coupled system of FODEs with boundary conditions, we consider

$$\begin{array}{c}\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)+\mathit{az}\left(t\right)+\mathit{by}\left(t\right)=f\left(t\right),& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)+\mathit{cy}\left(t\right)+\mathit{dz}\left(t\right)=g\left(t\right),& \end{array}\right.\\ z\left(0\right)=z_0,y\left(0\right)=y_0,z\left(1\right)=z_1,y\left(1\right)=y_1.\end{array}$$E23

Let us assume

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)={\mathrm{Ł}}_M{\psi }_M^T\left(t\right),& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)=K_M{\psi }_M^T\left(t\right).& \end{array}\right.$$E24

Applying Lemma 1 to Eq. (24), we get

$$\left\{\begin{array}{ll}z\left(t\right)=e_0+e_1\left(t\right)+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)& \\ y\left(t\right)=d_0+d_1\left(t\right)+K_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right),& \end{array}\right.$$E25

where $d_0,d_1,e_0,e_1\in R.$ Using the initial conditions in Eq. (25), we have $e_0=z_0,$$d_0=y_0.$ On using boundary conditions, we have from Eq. (25)

$$\begin{array}{l}{\left.z\left(1\right)=z_0+e_1+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right|}_{t=1},\\ {\left.z\left(1\right)-z_0-{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right|}_{t=1}=e_1.\end{array}$$

Similarly

$$\begin{array}{l}{\left.y\left(1\right)=y_0+d_1+K_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right|}_{t=1},\\ {\left.y\left(1\right)-y_0-K_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right|}_{t=1}=d_1.\end{array}$$

Equation (25) implies that

$$\left\{\begin{array}{ll}{\left.z\left(t\right)=z_0+t\left(z_1-z_0\right)-t(L_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right|}_{t=1})+L_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)& \\ {\left.y\left(t\right)=y_0+t\left(y_1-y_0\right)-t(K_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right|}_{t=1})+K_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right).& \end{array}\right.$$E26

Let $z_0+t\left(z_1-z_0\right)\approx F_M^1{\mathrm{\Psi }}_M^T\left(t\right)$ and $y_0+t\left(y_1-y_0\right)\approx F_M^2{\psi }_M^T\left(t\right),$ with

$$\begin{array}{ll}& {\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)={\mathrm{Ł}}_M{Q}_{M\times M}^{\gamma ,z}{\mathrm{\Psi }}_M^T\left(t\right)\\ & {\mathit{tK}}_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right)=K_M{Q}_{M\times M}^{\gamma ,y}{\mathrm{\Psi }}_M^T\left(t\right).\end{array}$$E27

Hence Eq. (26) implies

$$\left\{\begin{array}{ll}z\left(t\right)=F_M^1{\mathrm{\Psi }}_M^T\left(t\right)-L_M{Q}_{M\times M}^{\gamma ,z}{\mathrm{\Psi }}_M^T\left(t\right)+L_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)& \\ y\left(t\right)=F_M^2{\mathrm{\Psi }}_M^T\left(t\right)-K_M{Q}_{M\times M}^{\gamma ,y}{\mathrm{\Psi }}_M^T\left(t\right)+K_M{G}_{M\times M}^{\star \gamma }{\mathrm{\Psi }}_M^T\left(t\right).& \end{array}\right.$$E28

approximating $f\left(t\right)$ and $g\left(t\right)$ such that

$$\left\{\begin{array}{ll}f\left(t\right)\approx F_M^3{\mathrm{\Psi }}_M^T\left(t\right)& \\ g\left(t\right)\approx F_M^4{\mathrm{\Psi }}_M^T\left(t\right).& \end{array}\right.$$E29

On using (24)–(29), system (23) can be written as

$$\left\{\begin{array}{ll}{L}_M{\mathrm{\Psi }}_M^T\left(t\right)+a\left(F_M^1{\mathrm{\Psi }}_M^T\left(t\right)-L_M{Q}_{M\times M}^{\gamma ,z}{\mathrm{\Psi }}_M^T\left(t\right)+L_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right)& \\ +b\left(F_M^2{\mathrm{\Psi }}_M^T\left(t\right)-K_M{Q}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right)-F_M^3{\mathrm{\Psi }}_M^T\left(t\right)=0& \\ K_M{\mathrm{\Psi }}_M^T\left(t\right)+c\left(F_M^2{\mathrm{\Psi }}_M^T\left(t\right)-K_M{Q}_{M\times M}^{\gamma ,y}{\mathrm{\Psi }}_M^T\left(t\right)+K_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right)& \\ +d\left(F_M^1{\mathrm{\Psi }}_M^T\left(t\right)-{\mathrm{Ł}}_M{Q}_{M\times M}^{\gamma ,z}{\mathrm{\Psi }}_M^T\left(t\right)+{\mathrm{Ł}}_M{G}_{M\times M}^{\gamma }{\mathrm{\Psi }}_M^T\left(t\right)\right)-F_M^4{\mathrm{\Psi }}_M^T\left(t\right)=0.& \end{array}\right.$$

On rearrangement of terms, the above equations give

$$\left\{\begin{array}{ll}{\mathrm{Ł}}_M\left(I_{M\times M}-{\mathit{aQ}}_{M\times M}^{\gamma ,z}+{\mathit{aG}}_{M\times M}^{\gamma }\right)+K_M\left(I_{M\times M}-{\mathit{bQ}}_{M\times M}^{\gamma ,y}+{\mathit{bG}}_{M\times M}^{\gamma }\right)& \\ +{\mathit{aF}}_M^1+{\mathit{bF}}_M^2-F_M^3=0& \\ K_M\left(I_{M\times M}-{\mathit{cQ}}_{M\times M}^{\gamma ,y}+{\mathit{cG}}_{M\times M}^{\gamma }\right)+{\mathrm{Ł}}_M\left(I_{M\times M}-{\mathit{dQ}}_{M\times M}^{\gamma ,z}+{\mathit{dG}}_{M\times M}^{\gamma }\right)& \\ +{\mathit{cF}}_M^2+{\mathit{dF}}_M^1-F_M^4=0.& \end{array}\right.$$

In matrix form, we can write

$$\begin{array}{ll}& \left[{\mathrm{Ł}}_M{K}_M\right]\left[\begin{array}{cc}{I}_{M\times M}-{\mathit{aQ}}_{M\times M}^{\gamma ,z}+{\mathit{aG}}_{M\times M}^{\gamma }& 0\\ 0& I_{M\times M}-{\mathit{cQ}}_{M\times M}^{\gamma ,y}+{\mathit{cG}}_{M\times M}^{\gamma }\end{array}\right]\\ & +\left[\begin{array}{cc}{L}_M& K_M\end{array}\right]\left[\begin{array}{cc}0& I_{M\times M}-{\mathit{dQ}}_{M\times M}^{\gamma ,z}+{\mathit{dG}}_{M\times M}^{\gamma }\\ I_{M\times M}-{\mathit{bQ}}_{M\times M}^{\gamma ,y}+{\mathit{bG}}_{M\times M}^{\gamma }& 0\end{array}\right]\\ & +\left[\begin{array}{c}{\mathit{aF}}_M^1+{\mathit{bF}}_M^2-F_M^3\\ {\mathit{cF}}_M^2+{\mathit{dF}}_M^1-F_M^4\end{array}\right]=0.\end{array}$$

We convert the system to algebraic equation by considering

$$\begin{array}{c}\overline{L}=\left[\begin{array}{cc}{I}_{M\times M}-{\mathit{aQ}}_{M\times M}^{\gamma ,z}+{\mathit{aG}}_{M\times M}^{\gamma }& 0\\ 0& I_{M\times M}-{\mathit{cQ}}_{M\times M}^{\gamma ,y}+{\mathit{cG}}_{M\times M}^{\gamma }\end{array}\right]\\ \overline{M}=\left[\begin{array}{cc}0& I_{M\times M}-{\mathit{dQ}}_{M\times M}^{\gamma ,z}+{\mathit{dG}}_{M\times M}^{\gamma }\\ I_{M\times M}-{\mathit{bQ}}_{M\times M}^{\gamma ,y}+{\mathit{bG}}_{M\times M}^{\gamma }& 0\end{array}\right]\\ \text{and}\overline{N}=\left[\begin{array}{c}{\mathit{aF}}_M^1+{\mathit{bF}}_M^2-F_M^3\\ {\mathit{cF}}_M^2+{\mathit{dF}}_M^1-F_M^4\end{array}\right].\end{array}$$

so that the system is of the form

$$X\overline{L}+X\overline{M}+\overline{N}=0,$$

and solving the given equation for the unknown matrix $X=\left[\begin{array}{cc}{L}_M& K_M\end{array}\right],$ we get the required solution.

Example 4. As an example, we consider the Caputo fractional differential equation for the coupled system with the boundary conditions as

$$\left\{\begin{array}{ll}{}_0^c{D}_t^{\gamma }z\left(t\right)+2z\left(t\right)-2y\left(t\right)-f\left(t\right)=0,& \\ {}_0^c{D}_t^{\gamma }y\left(t\right)-3y\left(t\right)+2z\left(t\right)-g\left(t\right)=0,& \\ z\left(0\right)=4z\left(1\right)=-4,& \\ y\left(0\right)=2,y\left(1\right)=-2.& \end{array}\right.$$

At $\gamma =2,$ the exact solutions are

$$\left\{\begin{array}{ll}z\left(t\right)=t^6+t^5+t^4-t^3+t+1,& \\ y\left(t\right)=t^7-t^6+t^5+t^4+t^3-t^2-t+1.& \end{array}\right.$$

where the source functions are given by

$$\left\{\begin{array}{ll}f\left(t\right)=-2t^7+4t^6+30t^4+16t^3+12t^2-2t+2& \\ g\left(t\right)=-3t^7+12t^6+35t^5-27t^4-19t^3+20t^2+9t-4.& \end{array}\right.$$

We approximate the solution at the considered method by taking scale level $M=5.$ One can see that numerical plot and exact solution plot coincide very well as shown in Figure 5. Similarly the absolute error has been plotted at the given scale $M=5$ in Figure 6, which is very low. The lowest value of absolute error $\parallel z_{\mathit{app}}-z_{\mathit{ex}}\parallel$ and $\parallel y_{\mathit{app}}-y_{\mathit{ex}}\parallel$ indicates efficiency of the proposed method. The table shows the comparison of errors for exact and approximate solutions for fixed scale level $M=5$ and order $\gamma =1.9.$ Further the absolute error has been recorded at different values of space variable in Table 2 which provides the information about efficiency of the proposed method.

Figure 5.

Plots of exact and approximate solution for Case 4, boundary value problem.

Figure 6.

Plots of absolute error for Case 4, boundary value problem.

$t$	Absolute error $\parallel z_{\mathit{app}}-z_{\mathit{ex}}\parallel$	CPU time (s)	Absolute error $\parallel y_{\mathit{app}}-y_{\mathit{ex}}\parallel$	CPU time (s)
0	0.011	49.4	0.010	50.0
0.15	0.0062	50.3	0.0052	52.5
0.35	0.0058	51.2	0.0047	54.6
0.65	0.006	51.5	0.005	55.5
0.85	0.0075	52.6	0.007	56.4
1	0.011	53.8	0.010	56.2

Table 2.

Absolute error at different values of t for Example 4.

5. Conclusion

We have successfully used the class of orthogonal polynomials of Laguerre polynomials to establish a numerical method to compute the numerical solution of FODEs and their coupled systems under some initial and boundary conditions. By using these polynomials, we have obtained some operational matrices corresponding to fractional order derivatives and integration. Also we have computed a new matrix corresponding to boundary conditions for boundary value problems of FODEs. Using the aforementioned matrices, we have converted the considered problem of FODEs to Sylvester-type algebraic equations. To obtain the numerical solution, we easily solved the desired algebraic equations by taking help from MATLAB®. Corresponding to the established procedure, we have provided numbers of examples to demonstrate our results. Also some error analyses have been provided along with graphical representations. By increasing the scale level, the accuracy is increased and vice versa. On the other hand, when the fractional order is approaching to integer value, the solutions tend to the exact solutions of the considered FODE. Therefore in each example, we have compared the exact and approximate solution and found that both the solutions were in closure contact with each other. Hence the established method can be very helpful in solving many classes and systems of FODEs under both initial and boundary conditions. In future the shifted Laguerre polynomials can be used to compute numerical solutions of partial differential equations of fractional order.

Competing interests

We declare that no competing interests exist regarding this manuscript.

Author contribution

All authors equally contributed this paper and approved the final version.