Using Shifted Jacobi Polynomials to Handle Boundary Value Problems of Fractional Order

Abstract

This paper is concerned about the study of shifted Jacobi polynomials. By means of these polynomials, we construct some operational matrices of fractional order integration and differentiations. Based on these matrices, we develop a numerical scheme for the boundary value problems of fractional order differential equations. The construction of the procedure is new one for the coupled systems of fractional order boundary value problems. In the proposed scheme, we obtain a simple but highly accurate systems of algebraic equations. These systems are easily soluble by means of Matlab or using Mathematica. We provide some examples to which the procedure is applied to verify the applicability of our proposed method.

1. Introduction

Recently, fractional order differential equations have gained much attention from the researchers of mathematics, physics, computer science and engineers. This is due to the large numbers of applications of fractional order differential equations in various fields of science and applied nature. The applications of fractional order differential equations are found in physics, mechanics, viscoelasticity, photography, biology, chemistry, fluid mechanics, image and signal processing phenomenons, etc.; for more details, see [1, 2, 3, 4]. The researchers are giving more attention to fractional order differential equations, because in most cases the fractional order models are more accurate and reliable than classical order models. On the other hand, fractional order models have more degrees of freedom than integer order models. In most cases, the hereditary phenomena and memory description process are accurately modeled by means of fractional order derivatives and integrations as compare with integer order derivative and integrations. Due to these facts, researchers had given much attention to study the existence and uniqueness of positive solutions as well as multiplicity of solutions for fractional order differential equations and their systems. This area is very well explored, and large numbers of research articles can be found in literature; for details, see [5, 6, 7, 8, 9, 10]. It is very difficult to obtain exact solutions for such type of differential equations in all cases always. However, the area involving numerical solutions is in initial stage. Therefore, various numerical schemes have been developed in last few years by different researchers. The aim of these schemes was to achieve better accuracy. These schemes have their own merits and demerits. For example, in [11, 12, 13], the authors developed Adomain decomposition methods, Homotopy analysis methods, and Variational iterational methods to obtained good approximate solutions for certain fractional order differential equations. In recent years, spectral methods have got much attention as they were applied to solve some real-word problems of various fields of science and engineering. High accuracy was obtained for solving such problems [14]. Spectral method needs operational matrices for the numerical solutions, which have been constructed by using some polynomials, for example, in [15], the authors developed an operational matrix for shifted Legendre polynomials corresponding to fractional order derivative. In [16], the authors introduced operational matrix by using shifted Legendre polynomials corresponding to fractional order integrations. Similarly in [17, 18], authors constructed operational matrices for fractional order derivative by using Chebyshev and Jacobi polynomials. In all these cases, these matrices were applied to solve multiterms fractional order differential equations together with Tau-collocation method. In [19], Singh et al. derived an operational matrices for Berestein polynomials corresponding to fractional order integrations. Haar wavelet operational matrices were also developed, and some problems of fractional order differential equations were solved. In all these methods, collocation method was used together with these methods to obtain numerical solutions for fractional order differential equations. They only solved single problems.

Most of the physical and biological phenomena can be modeled in the form of coupled systems. For example, when two masses and two spring systems are modeled, it led to a coupled system of differential equations.

$$\begin{array}{l}{\mathrm{D}}^{\alpha }X+\frac{\lambda }{m_1}{\mathrm{D}}^{{\nu }_1}X-\frac{\mu }{m_1}{\mathrm{D}}^{{\nu }_2}Y+\frac{{\kappa }_1}{m_1}X-\frac{{\kappa }_2}{m_1}Y=f\left(x\right),1<\alpha \le 2,0<{\nu }_1,{\nu }_2\le 1,\\ {\mathrm{D}}^{\beta }Y-\frac{\lambda }{m_2}{\mathrm{D}}^{{\omega }_1}X+\frac{\lambda +\mu }{m_2}{\mathrm{D}}^{{\omega }_2}Y-\frac{{\kappa }_1}{m_2}X+\frac{{\kappa }_1+{\kappa }_2}{m_2}Y=g\left(x\right),1<\beta \le 2,0<{\omega }_1,{\omega }_2\le 1,\end{array}$$

subject to the boundary conditions

$$\begin{array}{ll}X\left(0\right)=0,& X\left(1\right)=0\\ Y\left(0\right)=1,& Y\left(1\right)=-1.\end{array}$$

Many coupled systems can be found in physics and in biological model, chemistry, fluid dynamics, etc.; see [16, 20, 21, 22, 23, 24, 25, 26, 27]. In fractional calculus, we generalize these systems by replacing the integer order by any fractional values that lie in certain range. All the above numerical methods were applied to solve singles problems, and to the best of our knowledge, very few articles are devoted to the numerical solutions of coupled systems or models. In all these cases, the authors have solved initial value problems of coupled systems of fractional order differential and partial differential equations; for details, see [15, 28] and the references therein. For example, in [15], authors have solved some initial value problems of coupled systems of fractional differential equations (FDEs) by using shifted Jacobi polynomials operational matrix method. However, the scheme is not properly applied yet for solving coupled systems of boundary value problems of fractional order differential equations. In this paper, we discuss the shifted Jacobi polynomial operational matrices methods to solve boundary value problems for a coupled systems of fractional order differential equations. In this scheme, we introduced a new matrix corresponding to boundary conditions, which is required for the approximate solutions. First, we solve single boundary value problems given by

$${\mathrm{D}}^{\alpha }x\left(t\right)+A{\mathrm{D}}^{\nu }x\left(t\right)+B{\mathrm{D}}^{\omega }x\left(t\right)+\mathit{Cx}\left(t\right)=f\left(t\right),t\in \left[0,T\right],A,B,C\in \mathbb{R},$$E1

subject the boundary conditions

$$x\left(0\right)=x_0,x\left(T\right)=x_1,$$

where $f\left(t\right)\in C\left[0,T\right]$ is a source term and $1<\alpha \le 2,0<\nu ,\omega <1.$ Then, we extend our technique to solve coupled system of fractional differential equations (FDEs) given as

$$\begin{array}{l}{\mathrm{D}}^{\alpha }x\left(t\right)+A_1{\mathrm{D}}^{{\nu }_1}x\left(t\right)+B_1{\mathrm{D}}^{{\nu }_2}y\left(t\right)+C_{1}x\left(t\right)+D_{1}y\left(t\right)=f\left(t\right),t\in \left[0,T\right],\\ {\mathrm{D}}^{\beta }y\left(t\right)+A_2{\mathrm{D}}^{{\omega }_1}x\left(t\right)+B_2{\mathrm{D}}^{{\omega }_2}y\left(t\right)+C_{2}x\left(t\right)+D_{2}y\left(t\right)=g\left(t\right),t\in \left[0,T\right],\end{array}$$E2

subject to the boundary conditions given by

$$\begin{array}{l}x\left(0\right)=x_0,x\left(T\right)=x_1,\\ y\left(0\right)=y_0,y\left(T\right)=y_1,\end{array}$$E3

where $1<\alpha ,\beta \le 2,0<{\nu }_i,{\omega }_i<1\left(i=1,2\right).$ Further, $f\left(t\right),g\left(t\right)\in C\left[0,T\right]$ are the source term of the system (42), (43) and $x_0,y_0,x_1,y_1$ are real constants, $A_i,B_i,C_i,D_i\left(i=1,2\right)$ are any constants. We convert the differential equations/systems of differential equations to a coupled systems of simple algebraic equations of Selvester or Laypunove type, which are easily soluble for unknown matrix required for approximations. In our proposed method, we do no need collation method, which is required for abovementioned methods. We applied our technique to some practical problems to verify the applicability of the proposed method.

Our article is organized as follows: Some fundamental concepts and results of fractional calculus and Jacobi polynomials needed are provided in sections 2. Section 3 contains some operational matrices for shifted Jacobi polynomials corresponding to fractional order derivative and integrations. Also in this section, we provide new operational matrix corresponding to the boundary conditions. In Section 4, we derive main procedure for the general class of boundary value problems of fractional order differential equations as well as for coupled system of FDEs. In Section 5, we apply the proposed method to some examples. Also we plot the images of these examples. In Section 6, we provide a short conclusion of the papers.

2. Preliminaries

In this section, we recall some basic definitions and results of fractional calculus, also we will give fundamental properties, definitions related to Jacobi polynomials as in [1, 4, 6, 7, 8, 14, 17, 18, 28, 29, 30, 31].

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

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

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

$${\mathbf{I}}^{\alpha }{z}^{\delta }=\frac{\mathrm{\Gamma }\left(\delta +1\right)}{\mathrm{\Gamma }\left(\delta +\alpha +1\right)}{z}^{\alpha +\delta }.$$E5

Definition 2.2. The Caputo fractional derivative of order $\alpha >0$ of a function $z\in C^n\left[0,1\right]$ is defined by

$$\begin{array}{l}{\mathrm{D}}_{0+}^{\alpha }z\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(n-\alpha \right)}{\int }_0^t{\left(t-s\right)}^{n-\alpha -1}\frac{d^n}{{\mathit{dt}}^n}z\left(s\right)\mathit{ds},n-1<\alpha \le n,t>0,\\ \mathit{\text{where}}n=⌈\alpha ⌉+1,\end{array}$$E6

provided that the right side is pointwise defined on $\left(0,\infty \right)$. Also one of the important properties of fractional order derivative is given by

$${\mathrm{D}}^{\alpha }{t}^k=\frac{\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(k-\alpha +1\right)}{t}^{k-\alpha }.\mathit{\text{Also for}}\mathit{any}\mathit{\text{cons}}\mathit{tan}tC,\mathit{\text{we have}}{\mathrm{D}}^{\alpha }C=0.$$E7

The following results are needed in the sequel.

Lemma 2.2.1. [6], Let $\alpha >0$ then

$${\mathrm{I}}^{\alpha }{\mathrm{D}}^{\alpha }y\left(t\right)=y\left(t\right)+C_0+C_{1}t+\dots +C_{n-1}{t}^{\alpha -n},$$E8

for arbitrary

$$C_i\in \mathbb{R},i=\mathrm{0,1,2},\dots ,n-1,n=\left[\alpha \right]+1.$$

2.1 The shifted Jacobi polynomials and its fundamental properties

In this section, we provide basic properties of shifted Jacobi polynomials. The famous Jacobi polynomials $P_i^{\left(a,b\right)}\left(z\right)$ are defined over the interval $\left[-1,1\right]$ as given by

$$\begin{array}{l}{P}_i^{\left(a,b\right)}\left(z\right)=\frac{\left(a+b+2i-1\right)\left[a^2-b^2+z\left(a+b+2i-2\right)\left(a+b+2i-2\right)\right]}{2i\left(a+b+i\right)\left(a+b+2i-2\right)}{P}_{i-1}^{\left(a,b\right)}\left(z\right)\\ -\frac{\left(a+i-1\right)\left(b+i-1\right)\left(a+b+2i\right)}{i\left(a+b+i\right)\left(a+b+2i-2\right)}{P}_{i-2}^{\left(a,b\right)}\left(z\right),i=2,3,\dots ,\\ \text{where}{P}_0^{\left(a,b\right)}\left(z\right)=1,P_1^{\left(a,b\right)}\left(z\right)=\frac{a+b+2}{2}z+\frac{a-b}{2}.\end{array}$$E9

By means of the substitution $z=\frac{2t}{\xi }-1$, we get a new version of the above polynomials defined in (9), which is called the shifted Jacobi polynomials over the interval $\left[0,\xi \right]$. The analytical form of this shifted Jacobi polynomials is given by the relations

$$\begin{array}{l}{P}_{\xi ,i}^{\left(a,b\right)}\left(t\right)=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\mathrm{\Gamma }\left(i+b+1\right)\mathrm{\Gamma }\left(i+k+a+b+1\right)}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+1\right){\xi }^k}{t}^k,\\ \text{where}{P}_{\xi ,i}^{\left(a,b\right)}\left(0\right)={\left(-1\right)}^i\frac{\mathrm{\Gamma }\left(i+b+1\right)}{\mathrm{\Gamma }\left(b+1\right)\mathrm{\Gamma }\left(i+1\right)}.\end{array}$$E10

The orthogonality conditions of the shifted Jacobi polynomials are given by

$$\begin{array}{l}{\int }_0^{\xi }{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right)P_{\xi ,i}^{\left(a,b\right)}\left(t\right)W_{\xi }^{\left(a,b\right)}\left(t\right)\mathit{dt}={\mathbf{\Omega }}_{\xi ,j}^{\left(a,b\right)}{\delta }_{\mathit{ji}},\\ \text{where}{\delta }_{\mathit{ji}}=1,\text{if}i=j,\text{other wise}{\delta }_{\mathit{ji}}=0,\\ \text{and the weight function is given}\mathrm{by}{W}_{\xi }^{\left(a,b\right)}\left(t\right)={\left(\xi -t\right)}^a{t}^b,\end{array}$$E11

and

$${\mathbf{\Omega }}_{\xi ,j}^{\left(a,b\right)}\left(t\right)=\frac{{\xi }^{a+b+1}\mathrm{\Gamma }\left(j+a+1\right)\mathrm{\Gamma }\left(j+b+1\right)}{\left(2j+a+b+1\right)\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(j+a+b+1\right)}.$$E12

Any function $f\left(t\right)$ such that $f^2\left(t\right)$ is integrable over $\left[0,\xi \right],$ can be approximated in terms of shifted Jacobi polynomials as given by

$$f\left(t\right)\simeq \sum _{k=0}^m{B}_j{P}_{\xi ,k}^{\left(a,b\right)}\left(t\right)={\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right),$$E13

where the shifted Jacobi coefficient vector is denoted by ${\mathrm{K}}_M$ and ${\mathrm{\Phi }}_M\left(t\right)$ is $M$ terms function vector. Also $M=m+1,$ when $m\to \infty$ the approximation converges to the exact value of the function. The coefficient $b_k$ can be calculated by using (11)–(13) as

$$B_j=\frac{1}{{\mathrm{\Omega }}_{\xi ,j}^{\left(a,b\right)}}{\int }_0^{\xi }{W}_{\xi }^{\left(a,b\right)}\left(t\right)f\left(t\right)P_{\xi ,i}^{\left(a,b\right)}\left(t\right)\mathit{dt},i=0,1,.\dots$$E14

3. Operational matrices

In this section, we provide some operational matrices for fractional order differentiation as well as for fractional order integration. Further, in this section, we give construction of new matrix for the boundary conditions.

Theorem 3.1. Let ${\mathrm{\Phi }}_M\left(t\right)$ be the function vector, then the operational matrix of fractional order derivative $\alpha$ is given by

$${\mathrm{D}}^{\alpha }\left[{\mathrm{\Phi }}_M\left(t\right)\right]\simeq {\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right),$$E15

where ${\mathrm{H}}_{M\times M}^{\left(\alpha \right)}$ is the operational matrix of fractional order derivative $\alpha$ given by

$${\mathbf{H}}_{M\times M}^{\left(\alpha \right)}=\left[\begin{array}{cccccc}{\mathrm{\Upsilon }}_{0,0,k}& {\mathrm{\Upsilon }}_{0,1,k}& \cdots & {\mathrm{\Upsilon }}_{0,j,k}& \cdots & {\mathrm{\Upsilon }}_{0,m,k}\\ {\mathrm{\Upsilon }}_{2,1,k}& {\mathrm{\Upsilon }}_{2,2,k}& \cdots & {\mathrm{\Upsilon }}_{2,r,k}& \cdots & {\mathrm{\Upsilon }}_{1,m,k}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathrm{\Upsilon }}_{i,0,k}& {\mathrm{\Upsilon }}_{i,1,k}& \cdots & {\mathrm{\Upsilon }}_{i,j,k}& \cdots & {\mathrm{\Upsilon }}_{i,m,k}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathrm{\Upsilon }}_{m,0,k}& {\mathrm{\Upsilon }}_{m,1,k}& \cdots & {\mathrm{\Upsilon }}_{m,j,k}& \cdots & {\mathrm{\Upsilon }}_{m,m,k}\end{array}\right],$$E16

where

$$\begin{array}{l}{\mathrm{\Upsilon }}_{i,j,k}=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\mathrm{\Gamma }\left(i+b+1\right)\mathrm{\Gamma }\left(i+k+a+b+1\right)\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+1\right)\mathrm{\Gamma }\left(k+1-\alpha \right){\xi }^k}\\ \times \sum _{l=0}^j\frac{{\left(-1\right)}^{j-l}\left(2j+a+b+1\right)\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(j+l+a+b+1\right)\mathrm{\Gamma }\left(k-\mathit{\text{alpha}}+l+b+1\right)\mathrm{\Gamma }\left(a+1\right){\xi }^{\alpha }}{\mathrm{\Gamma }\left(j+a+1\right)\mathrm{\Gamma }\left(l+b+1\right)\mathrm{\Gamma }\left(j-l+1\right)\mathrm{\Gamma }\left(l+1\right)\mathrm{\Gamma }\left(k-\alpha +l+b+a+2\right)}.\\ \mathit{\text{If}}{\mathrm{\Upsilon }}_{i,j,k}=0,\mathit{\text{for}}i<\alpha .\end{array}$$E17

Proof. The proof of this theorem is same as given in [14, 17] in lemma $\left(3.2\right)$ and $\left(3.4\right)$, respectively.□

Theorem 3.2. Let ${\mathbf{\Phi }}_M\left(t\right)$ be the function vector corresponding to the shifted Jacobi polynomials and $\alpha >0$_, then the operational matrix corresponding to the fractional order integration is given by

$${\mathrm{I}}^{\alpha }\left[{\mathrm{\Phi }}_M\left(t\right)\right]\simeq {\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right),$$E18

where ${\mathrm{H}}_{M\times M}^{\left(\alpha \right)}$ is the operational matrix of fractional order integration $\alpha$ and is given by

$${\mathbf{H}}_{M\times M}^{\left(\alpha \right)}=\left[\begin{array}{cccccc}{\mathrm{\Upsilon }}_{0,0,k}& {\mathrm{\Upsilon }}_{0,1,k}& \cdots & {\mathrm{\Upsilon }}_{0,j,k}& \cdots & {\mathrm{\Upsilon }}_{0,m,k}\\ {\mathrm{\Upsilon }}_{2,1,k}& {\mathrm{\Upsilon }}_{2,2,k}& \cdots & {\mathrm{\Upsilon }}_{2,r,k}& \cdots & {\mathrm{\Upsilon }}_{1,m,k}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathrm{\Upsilon }}_{i,0,k}& {\mathrm{\Upsilon }}_{i,1,k}& \cdots & {\mathrm{\Upsilon }}_{i,j,k}& \cdots & {\mathrm{\Upsilon }}_{i,m,k}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathrm{\Upsilon }}_{m,0,k}& {\mathrm{\Upsilon }}_{m,1,k}& \cdots & {\mathrm{\Upsilon }}_{m,j,k}& \cdots & {\mathrm{\Upsilon }}_{m,m,k}\end{array}\right],$$E19

where

$$\begin{array}{l}{\mathrm{\Upsilon }}_{i,j,k}=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\mathrm{\Gamma }\left(i+b+1\right)\mathrm{\Gamma }\left(i+k+a+b+1\right)\mathrm{\Gamma }\left(k+1\right)}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+1\right)\mathrm{\Gamma }\left(k+1+\alpha \right){\xi }^k}\\ \times \sum _{l=0}^j\frac{{\left(-1\right)}^{j-l}\left(2j+a+b+1\right)\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(j+l+a+b+1\right)\mathrm{\Gamma }\left(k+\alpha +l+b+1\right)\mathrm{\Gamma }\left(a+1\right){\xi }^{\alpha }}{\mathrm{\Gamma }\left(j+a+1\right)\mathrm{\Gamma }\left(l+b+1\right)\mathrm{\Gamma }\left(j-l+1\right)\mathrm{\Gamma }\left(l+1\right)\mathrm{\Gamma }\left(k+\alpha +l+b+a+2\right)}.\end{array}$$E20

Proof. The proof of this theorem is available in [14]. Therefore, we omit the proof.

Theorem 3.3 Let ${\mathrm{\Phi }}_M\left(t\right)$ be a function vector, and let $\varphi \left(t\right)$ be any function defined as $\varphi \left(t\right)={\mathit{ct}}^n,n=\mathrm{0,1,2},\dots ,c\in \mathbb{R}$ and $u\left(t\right)={\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right)$ . Then

$$\varphi {\left(t\right)}_0{\mathrm{I}}_{\xi }^{\alpha }u\left(t\right)={\mathrm{K}}_M^T{\mathrm{Q}}_{M\times M}^{\left(c,\varphi ,\alpha \right)}{\mathrm{\Phi }}_M\left(t\right),$$E21

where ${\mathrm{Q}}_{M\times M}^{\left(c,\varphi ,\alpha \right)}$ is an operational matrix corresponding to some boundary value and is given by

$${\mathbf{Q}}_{M\times M}^{\left(c,\varphi ,\alpha \right)}=\left[\begin{array}{cccccc}{\mathrm{\Upsilon }}_{0,0,k}& {\mathrm{\Upsilon }}_{0,1,k}& \cdots & {\mathrm{\Upsilon }}_{0,j,k}& \cdots & {\mathrm{\Upsilon }}_{0,m,k}\\ {\mathrm{\Upsilon }}_{2,1,k}& {\mathrm{\Upsilon }}_{2,2,k}& \cdots & {\mathrm{\Upsilon }}_{2,r,k}& \cdots & {\mathrm{\Upsilon }}_{1,m,k}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathrm{\Upsilon }}_{i,0,k}& {\mathrm{\Upsilon }}_{i,1,k}& \cdots & {\mathrm{\Upsilon }}_{i,j,k}& \cdots & {\mathrm{\Upsilon }}_{i,m,k}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathrm{\Upsilon }}_{m,0,k}& {\mathrm{\Upsilon }}_{m,1,k}& \cdots & {\mathrm{\Upsilon }}_{m,j,k}& \cdots & {\mathrm{\Upsilon }}_{m,m,k}\end{array}\right],$$E22

where

$${\mathrm{\Upsilon }}_{i,j,k}=\sum _{j=0}^M{\mathrm{\Delta }}_{i,k,\alpha }{B}_j{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right),$$E23

where

$${\mathrm{\Delta }}_{i,k,\alpha }=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)c{\xi }^{\alpha }}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\alpha \right)}$$E24

and

$$B_j=\sum _{l=0}^j\frac{{\left(-1\right)}^{j-l}\left(2j+a+b+1\right)\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(j+l+a+b+1\right)\mathrm{\Gamma }\left(n+b+1\right)\mathrm{\Gamma }\left(a+1\right){\xi }^{n-l}}{\mathrm{\Gamma }\left(j+a+1\right)\mathrm{\Gamma }\left(l+b+1\right)\mathrm{\Gamma }\left(j-l+1\right)\mathrm{\Gamma }\left(l+1\right)\mathrm{\Gamma }\left(n+a+b+1\right)}.$$E25

Proof. Consider the general term of ${\mathrm{\Phi }}_M\left(t\right)$ as

$$\begin{array}{l}{}_0{\mathrm{I}}_{\xi }^{\alpha }{P}_{\xi ,i}^{\left(a,b\right)}\left(t\right)=\frac{1}{\mathrm{\Gamma }\alpha }{\int }_0^{\xi }{\left(\xi -s\right)}^{\alpha -1}{P}_{\xi ,i}^{\left(a,b\right)}\left(s\right)\mathit{ds}\\ =\frac{1}{\mathrm{\Gamma }\alpha }{\int }_0^{\xi }{\left(\xi -s\right)}^{\alpha -1}\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+k+a+b+1\right)s^k}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+1\right){\xi }^k}\mathit{ds}\\ =\frac{1}{\mathrm{\Gamma }\alpha }\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+k+a+b+1\right)}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+1\right){\xi }^k}{\int }_0^{\xi }{\left(\xi -s\right)}^{\alpha -1}{s}^k\mathit{ds}.\end{array}$$E26

Now by applying Laplace transform, we get

$$\begin{array}{l}\mathcal{L}\left({\int }_0^{\xi }{\left(\xi -s\right)}^{\alpha -1}{s}^k\mathit{ds}\right)=\frac{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(k+1\right)}{s^{k+\alpha +1}}\\ \Rightarrow {\int }_0^{\xi }{\left(\xi -s\right)}^{\alpha -1}{s}^k\mathit{ds}=\frac{\mathrm{\Gamma }\left(k+1\right)\mathrm{\Gamma }\alpha {\xi }^{k+\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}.\end{array}$$E27

Now putting (27) in (26), we get

$${}_0{\mathrm{I}}_{\xi }^{\alpha }{P}_{\xi ,i}^{\left(a,b\right)}\left(t\right)=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+a+b+k+1\right){\xi }^{\alpha }}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\alpha \right)}.$$E28

Now if $\varphi \left(t\right)=t^n$, then

$$\mathit{c\varphi }{\left(t\right)}_0{\mathrm{I}}_{\xi }^{\alpha }{P}_{\xi ,i}^{\left(a,b\right)}\left(t\right)=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+a+b+k+1\right){\xi }^{\alpha }{\mathit{ct}}^n}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\alpha \right)}.$$E29

Representing $t^n$ in terms of shifted Jacobi polynomials as

$$t^n=\sum _{j=0}^M{B}_j{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right).$$E30

By the use of orthogonality relation (11), we have

$$\begin{array}{l}{B}_j=\frac{1}{{\mathrm{\Omega }}_{\xi ,j}^{\left(a,b\right)}}{\int }_0^{\xi }{t}^n{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right)W_{\xi }^{\left(a,b\right)}\left(\xi \right)\mathit{dt}\\ =\frac{1}{{\mathrm{\Omega }}_{\xi ,j}^{\left(a,b\right)}}\sum _{l=0}^j\frac{{\left(-1\right)}^{\left(j-l\right)}\mathrm{\Gamma }\left(j+b+1\right)\mathrm{\Gamma }\left(j+l+a+b+1\right)}{\mathrm{\Gamma }\left(l+b+1\right)\mathrm{\Gamma }\left(j+a+b+1\right)\mathrm{\Gamma }\left(j-l+1\right)\mathrm{\Gamma }\left(l+1\right){\xi }^l}{\int }_0^{\xi }{t}^{n+b}{\left(\xi -t\right)}^a\mathit{dt}.\end{array}$$E31

Now by means of Laplace transform, we have

$$\begin{array}{l}\mathrm{L}\left({\int }_0^{\xi }{t}^{n+b}{\left(\xi -t\right)}^a\mathit{dt}\right)=\frac{\mathrm{\Gamma }\left(n+b+1\right)\mathrm{\Gamma }\left(a+1\right)}{{\xi }^{n+b+a+2}}\\ \Rightarrow {\int }_0^{\xi }{t}^{n+b}{\left(\xi -t\right)}^a\mathit{dt}=\frac{\mathrm{\Gamma }\left(n+b+1\right)\mathrm{\Gamma }\left(a+1\right)}{\mathrm{\Gamma }\left(n+a+b+1\right)}{\xi }^{n+a+b+1}.\end{array}$$E32

Thus putting (32) in (31) and using (12), we get

$$B_j=\sum _{l=0}^j\frac{{\left(-1\right)}^{j-l}\left(2j+a+b+1\right)\mathrm{\Gamma }\left(j+1\right)\mathrm{\Gamma }\left(j+l+a+b+1\right)\mathrm{\Gamma }\left(n+b+1\right)\mathrm{\Gamma }\left(a+1\right){\xi }^{n-l}}{\mathrm{\Gamma }\left(j+a+1\right)\mathrm{\Gamma }\left(l+b+1\right)\mathrm{\Gamma }\left(j-l+1\right)\mathrm{\Gamma }\left(l+1\right)\mathrm{\Gamma }\left(n+a+b+1\right)}.$$E33

Now (29) implies that

$$\begin{array}{l}\varphi {\left(t\right)}_0{\mathrm{I}}_{\xi }^{\alpha }{P}_{\xi ,i}^{\left(a,b\right)}\left(t\right)=\sum _{j=0}^M\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+a+b+k+1\right){\xi }^{\alpha }c}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i-k+1\right)\mathrm{\Gamma }\left(k+\alpha \right)}{B}_j{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right)\\ =\sum _{j=0}^M{\mathrm{\Delta }}_{i,k,\alpha }{B}_j{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right)=\sum _{j=0}^M{\mathrm{\Upsilon }}_{i,j,k}{P}_{\xi ,j}^{\left(a,b\right)}\left(t\right),\end{array}$$E34

where

$$\begin{array}{l}{\mathrm{\Upsilon }}_{i,j,k}=\sum _{j=0}^M{\mathrm{\Delta }}_{i,k,\alpha }{B}_j,\\ {\mathrm{\Delta }}_{i,k,\alpha }=\sum _{k=0}^i\frac{{\left(-1\right)}^{i-k}\left(i+b+1\right)\mathrm{\Gamma }\left(i+a+b+k+1\right){\xi }^{\alpha }c}{\mathrm{\Gamma }\left(k+b+1\right)\mathrm{\Gamma }\left(i+a+b+1\right)\mathrm{\Gamma }\left(i+k+1\right)\mathrm{\Gamma }\left(k+\alpha \right)}.\end{array}$$

Evaluating the result (34) for different values of $j$ we get the required result.

4. Applications of operational matrices

In this section, we show fundamental importance of operational matrices of fractional order derivative and integration. We apply them to solve some multiterms fractional order boundary value problems of fractional differential equations. Consider the following general fractional differential equation with constant coefficient and given boundary conditions

$${\mathrm{D}}^{\alpha }x\left(t\right)+A{\mathrm{D}}^{\nu }x\left(t\right)+B{\mathrm{D}}^{\omega }x\left(t\right)+\mathit{Cx}\left(t\right)=f\left(t\right),t\in \left[0,T\right],A,B,C\in \mathbb{R},$$E35

subject the boundary conditions

$$x\left(0\right)=x_0,x\left(T\right)=x_1,$$

where $f\left(t\right)\in C\left[0,T\right]$ is a source terms and $1<\alpha \le 2,0<\nu ,\omega <1.$.

To obtain the solutions of the (35) in terms of shifted Jacobi polynomials, we proceed as

$$D^{\alpha }x\left(t\right)={\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right).$$E36

Applying ${\mathrm{I}}^{\alpha }$, definition (7), and lemma (2.2.1), then we get

$$x\left(t\right)={\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)+c_0+c_{1}t.$$E37

By means of boundary conditions we get

$${\left.x\left(t\right)={\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)+x_0+t\frac{\left(x_1-x_0\right)}{T}-\frac{t}{T}{\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)\right|}_{t=T}$$E38

using the approximation $x_0+t\frac{\left(x_1-x_0\right)}{T}\simeq {\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right),\frac{t}{T}=\varphi \left(t\right),$ also using ${\left.\mathit{c\varphi }\left(t\right){\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)\right|}_{t=T}={\mathrm{K}}_M^T{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}{\mathrm{\Phi }}_M\left(t\right),$then (38) implies that

$$\begin{array}{c}{\left.x\left(t\right)={\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1\mathrm{\Phi }\left(t\right)-\mathit{c\varphi }\left(t\right){\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)\right|}_{t=T}\\ ={\mathrm{K}}_M^T{\mathrm{H}}_{M\times M}^{\left(\alpha \right)}{\mathrm{\Phi }}_M\left(t\right)-{\mathrm{K}}_M^T{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right)\\ ={\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\alpha }-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right).\end{array}$$E39

Now from (39) we have

$$\begin{array}{l}{\mathrm{D}}^{\nu }x\left(t\right)={\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left(\nu \right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left(\nu \right)}{\mathrm{\Phi }}_M\left(t\right)\\ {\mathrm{D}}^{\omega }x\left(t\right)={\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left(\omega \right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left(\omega \right)}{\mathrm{\Phi }}_M\left(t\right)\\ \text{and approximating the source term}\mathrm{as}f\left(t\right)\simeq {\mathrm{F}}_M^2{\mathrm{\Phi }}_M\left(t\right).\end{array}$$E40

Putting (36)–(40) in (35), which yields

$$\begin{array}{l}{\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right)+A\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left(\nu \right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left(\nu \right)}{\mathrm{\Phi }}_M\left(t\right)\right]\\ +B\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left(\omega \right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left(\omega \right)}{\mathrm{\Phi }}_M\left(t\right)\right]\\ +C\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right)\right]-{\mathrm{F}}_M^2{\mathrm{\Phi }}_M\left(t\right)=0\\ {\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right)\left[A{\mathrm{G}}_{M\times M}^{\left(\nu \right)}+B{\mathrm{G}}_{M\times M}^{\left(\omega \right)}+\mathit{CI}\right]{\mathrm{\Phi }}_M\left(t\right)\\ +\left[{\mathrm{F}}_M^1\left(A{\mathrm{G}}_{M\times M}^{\left(\nu \right)}+B{\mathrm{G}}_{M\times M}^{\left(\omega \right)}\right)+C{\mathrm{F}}_M^1-{\mathrm{F}}_M^2\right]{\mathrm{\Phi }}_M\left(t\right)=0\\ {\mathrm{K}}_M^T+{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\alpha }-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right)\left[A{\mathrm{G}}_{M\times M}^{\left(\nu \right)}+B{\mathrm{G}}_{M\times M}^{\left(\omega \right)}+\mathit{CI}\right]\\ +{\mathrm{F}}_M^1\left(A{\mathrm{G}}_{M\times M}^{\left(\nu \right)}+B{\mathrm{G}}_{M\times M}^{\left(\omega \right)}+\mathit{CI}\right)-{\mathrm{F}}_M^2=0.\end{array}$$E41

Which is a simple algebraic equation. Solving this equation for ${\mathbf{K}}_M$ and putting it in (39) we get the required approximate solution of the corresponding boundary value problem.

4.1 Coupled system of boundary value problems for fractional order differential equations

In this subsection, we use operational matrices to derive procedure for the numerical solutions of coupled system. We consider the following general coupled system of FDEs as

$$\begin{array}{l}{\mathrm{D}}^{\alpha }x\left(t\right)+A_1{\mathrm{D}}^{{\nu }_1}x\left(t\right)+B_1{\mathrm{D}}^{{\nu }_2}y\left(t\right)+C_{1}x\left(t\right)+D_{1}y\left(t\right)=f\left(t\right),t\in \left[0,T\right],\\ {\mathrm{D}}^{\beta }y\left(t\right)+A_2{\mathrm{D}}^{{\omega }_1}x\left(t\right)+B_2{\mathrm{D}}^{{\omega }_2}y\left(t\right)+C_{2}x\left(t\right)+D_{2}y\left(t\right)=g\left(t\right),t\in \left[0,T\right],\end{array}$$E42

subject to the boundary conditions given by

$$\begin{array}{l}x\left(0\right)=x_0,x\left(T\right)=x_1,\\ y\left(0\right)=y_0,y\left(T\right)=y_1,\end{array}$$E43

where $1<\alpha ,\beta \le 2,0<{\nu }_i,{\omega }_i<1\left(i=1,2\right).$ Further, $f\left(t\right),g\left(t\right)\in C\left[0,T\right]$ are the source term of the system (42), (43) and $x_0,y_0,x_1,y_1$ are real constants, $A_i,B_i,C_i,D_i\left(i=1,2\right)$ are any constants. To approximate the above system in terms of shifted Jacobi polynomials, we take

$${\mathrm{D}}^{\alpha }x\left(t\right)={\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right)$$E44

and

$${\mathrm{D}}^{\beta }y\left(t\right)={\mathrm{L}}_M^T{\mathrm{\Phi }}_M\left(t\right)$$E45

Applying ${\mathrm{I}}^{\alpha },{\mathrm{I}}^{\beta }$ and lemma (1) and boundary conditions, the above Eqs. (44), (45) imply that

$$x\left(t\right)={\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right)$$E46

$$y\left(t\right)=\mathrm{L}{T}_M\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{\Phi }}_M\left(t\right).$$E47

Now taking fractional order derivative of (46), (47), we have

$$\begin{array}{l}{\mathrm{D}}^{{\nu }_1}x\left(t\right)=K_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\nu }_1\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left({\nu }_1\right)}{\mathrm{\Phi }}_M\left(t\right)\\ {\mathrm{D}}^{{\omega }_1}x\left(t\right)={\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\omega }_1\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left({\omega }_1\right)}{\mathrm{\Phi }}_M\left(t\right)\end{array}$$E48

and

$$\begin{array}{l}{\mathrm{D}}^{{\nu }_2}y\left(t\right)={\mathrm{L}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\nu }_2\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{G}}_{M\times M}^{\left({\nu }_2\right)}{\mathrm{\Phi }}_M\left(t\right)\\ {\mathrm{D}}^{\left({\omega }_2\right)}y\left(t\right)={\mathrm{L}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\omega }_2\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{G}}_{M\times M}^{\left({\omega }_2\right)}{\mathrm{\Phi }}_M\left(t\right).\end{array}$$E49

Putting (45)–(49) together with the use of approximation $f\left(t\right)\simeq {\mathrm{F}}_M^3{\mathrm{\Phi }}_M\left(t\right),g\left(t\right)\simeq {\mathrm{F}}_M^4{\mathrm{\Phi }}_M\left(t\right)$ in (42), we obtain

$$\begin{array}{l}{\mathrm{K}}_M^T{\mathrm{\Phi }}_M\left(t\right)+A_1\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\nu }_1\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left({\nu }_1\right)}{\mathrm{\Phi }}_M\left(t\right)\right]\\ +B_1\left[{\mathrm{L}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\nu }_2\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{G}}_{M\times M}^{\left({\nu }_2\right)}{\mathrm{\Phi }}_M\left(t\right)\right]\\ +C_1\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right)\right]\\ +D_1\left[{\mathrm{L}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{\Phi }}_M\left(t\right)\right]\\ -{\mathrm{F}}_M^3{\mathrm{\Phi }}_M\left(t\right)=0\end{array}$$E50

$$\begin{array}{l}{\mathrm{L}}_M^T{\mathrm{\Phi }}_M\left(t\right)+A_2\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\omega }_1\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{G}}_{M\times M}^{\left({\omega }_1\right)}{\mathrm{\Phi }}_M\left(t\right)\right]\\ +B_2\left[{\mathrm{L}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{G}}_{M\times M}^{\left({\omega }_2\right)}{\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{G}}_{M\times M}^{\left({\omega }_2\right)}{\mathrm{\Phi }}_M\left(t\right)\right]\\ +C_2\left[{\mathrm{K}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\alpha \right)}-{\mathrm{Q}}_{M\times M}^{\left(\alpha ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^1{\mathrm{\Phi }}_M\left(t\right)\right]\\ +D_2\left[{\mathrm{L}}_M^T\left({\mathrm{H}}_{M\times M}^{\left(\beta \right)}-{\mathrm{Q}}_{M\times M}^{\left(\beta ,\varphi \right)}\right){\mathrm{\Phi }}_M\left(t\right)+{\mathrm{F}}_M^2{\mathrm{\Phi }}_M\left(t\right)\right]\\ -{\mathrm{F}}_M^4{\mathrm{\Phi }}_M\left(t\right)=0.\end{array}$$

Rearranging the terms in above system (50) and using the notation for simplicity

$$\begin{array}{l}\mathrm{S}={\mathrm{F}}_M^1\left(A_1{\mathrm{G}}^{\left({\nu }_1\right)}+C_{1}I\right)+{\mathrm{F}}_M^2\left(B_1{\mathrm{G}}^{\left({\nu }_2\right)}+D_{1}I\right)-{\mathrm{F}}_M^3\\ \mathrm{R}={\mathrm{F}}_M^1\left(A_2{\mathrm{G}}^{\left({\omega }_1\right)}+C_{2}I\right)+{\mathrm{F}}_M^2\left(B_2{\mathrm{G}}^{\left({\omega }_2\right)}+D_{2}I\right)-{\mathrm{F}}_M^4\\ {\mathrm{H}}^{\left(\alpha \right)}-{\mathrm{Q}}^{\left(\alpha ,\varphi \right)}=\overbrace{P},{\mathrm{H}}^{\left(\beta \right)}-{\mathrm{Q}}^{\left(\beta ,\varphi \right)}=\overbrace{Q}\end{array}$$E50a

then witting in matrix form, we have

$$\begin{array}{l}\left[{\mathrm{K}}_M^T{\mathrm{L}}_M^T\right]\left[\begin{array}{cc}{\mathrm{\Phi }}_M\left(t\right)& 0\\ 0& {\mathrm{\Phi }}_M\left(t\right)\end{array}\right]\\ +\left[\begin{array}{cc}{\mathrm{K}}_M^T& {\mathrm{L}}_M^T\end{array}\right]\left[\begin{array}{cc}\overbrace{\mathrm{P}}\left(A_1{\mathrm{G}}^{{\nu }_1}+C_{1}I\right)& \overbrace{\mathrm{Q}}\left(B_1{\mathrm{G}}^{{\nu }_2}+D_{1}I\right)\\ \overbrace{\mathrm{P}}\left(A_2{\mathrm{G}}^{{\omega }_1}+C_{2}I\right)& \overbrace{\mathrm{Q}}\left(B_2{\mathrm{G}}^{{\omega }_2}+C_{2}I\right)\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Phi }}_M\left(t\right)& 0\\ 0& {\mathrm{\Phi }}_M\left(t\right)\end{array}\right]\\ +\left[\begin{array}{cc}\mathrm{S}& \mathrm{R}\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Phi }}_M\left(t\right)& 0\\ 0& {\mathrm{\Phi }}_M\left(t\right)\end{array}\right]=0\\ \Rightarrow \left[\begin{array}{cc}{\mathrm{K}}_M^T& {\mathrm{L}}_M^T\end{array}\right]+\left[\begin{array}{cc}{\mathrm{K}}_M^T& {\mathrm{L}}_M^T\end{array}\right]\left[\begin{array}{cc}\overbrace{\mathrm{P}}\left(A_1{\mathrm{G}}^{{\nu }_1}+C_{1}I\right)& \overbrace{\mathrm{Q}}\left(B_1{\mathrm{G}}^{{\nu }_2}+D_{1}I\right)\\ \overbrace{\mathrm{P}}\left(A_2{\mathrm{G}}^{{\omega }_1}+C_{2}I\right)& \overbrace{\mathrm{Q}}\left(B_2{\mathrm{G}}^{{\omega }_2}+C_{2}I\right)\end{array}\right]+\left[\begin{array}{cc}\mathrm{S}& \mathrm{R}\end{array}\right]=0\end{array}$$E51

which is an algebraic equation and can easily be solved for unknown matrices ${\mathrm{K}}_M,{\mathrm{L}}_M$, then putting it in (38) and (39), we get required approximations for $x\left(t\right),y\left(t\right)$.

5. Numerical examples

In this section, we will apply our proposed scheme to some practical problems.

Example 1. Consider the problems of Buckling of a thin vertical column such that when a constant compressive force $F$ is applied to the said thin column of length $L$ hanged at both ends whose general differential equations in the absence of source term are given by

$$\begin{array}{l}\mathit{EL}\frac{d^{\alpha }y}{{\mathit{dx}}^{\alpha }}+\mathit{Fy}=f\left(t\right),1<\alpha \le 2,\\ \mathit{\text{subject to the boundary conditions}}y\left(0\right)=y_1,y\left(L\right)=y_2,\end{array}$$E52

where $E,I$ are Young modulus and moment of inertia, respectively. Where $E=2.6\times {10}^7\mathit{lb}/\mathit{\text{in}},I=0.25{\pi }^2{R}^4,$ where $R$ is the radius of the column. Taking $L=1\mathit{\text{in}},R=1\mathit{\text{in}}$ the exact deflection at $\alpha =2$ is $y\left(x\right)=sin\left(\mathit{\pi x}\right)$. We apply our technique introduced in section $\left(4\right)$ for various choices of $a,b,M$. The observations are given in Table 1 for the absolute error $\mid y_{\mathit{app}}-y_{\mathit{\text{exact}}}\mid$.

$M$	$a$	$b$	$\mid y_{\mathit{app}}-y_{\mathit{\text{exact}}}\mid$	$a$	$b$	$\mid y_{\mathit{app}}-y_{\mathit{\text{exact}}}\mid$
$2$	$0$	$0$	$8\times {10}^{-2}$	$1$	$1$	$4\times {10}^{-2}$
$4$	$0$	$0$	$6\times {10}^{-3}$	$1$	$1$	$3\times {10}^{-3}$
$6$	$2$	$2$	$1.2\times {10}^{-3}$	$1$	1	$2\times {10}^{-3}$
$8$	$0.5$	$0.5$	$8\times {10}^{-4}$	$1.5$	$1.5$	$2\times {10}^{-4}$
$8$	$0$	$0$	$1\times {10}^{-5}$	$-0.5$	$-0.5$	$2\times {10}^{-5}$
$10$	$2$	$0$	$2.3\times {10}^{-5}$	$0$	$2$	$4\times {10}^{-5}$

Table 1.

Maximum absolute error for $T=1,\alpha =1.8,{\nu }_1={\nu }_2=0$ for example, (1).

In Table 2, we have shown the maximum absolute error for various choices of $a,b,M$. Also in Figures 1 and 2, we have shown the images of their comparison between exact and approximate values. Further, we observed that as the $\alpha \to 2$, the error is reducing and one can check that at $\alpha =2$, the absolute error is below $1.8\times {10}^{-13}$ as given in the next table. From above tables, we observe that as order $\alpha \to 2$ the approximation converges to the value at $\alpha =2$. This phenomenon shows the accuracy of the approximate solutions.

$M$	$\alpha$	${\nu }_1$	${\nu }_2$	$\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$
$4$	$2$	$1$	$1$	$1\times {10}^{-3}$
$6$	$2$	$1$	$1$	$7\times {10}^{-17}$
	$1.8$	$0.9$	$0.8$	$2\times {10}^{-4}$
$8$	$1.7$	$1$	$0.8$	$2\times {10}^{-4}$
$10$	$1.9$	$0.999$	$0.999$	$3\times {10}^{-5}$

Table 2.

Maximum absolute error at $a=b=0,$ and for different $\alpha ,{\nu }_1,{\nu }_2$ for example, (2).

Figure 1.

Subplot (a) represents comparison of approximate and exact $y\left(x\right)$ of example 1. Stars and dots represent exact solution while the red curve represents approximate solution. Setting $\alpha =1.8,M=4,T=1,a=b=1.$ the subplot represents their absolute error at $M=8$.

Figure 2.

Subplot (a) represents comparison of approximate and exact $y\left(x\right)$ of example 1.Red dots represent exact solution while the blue curve represents approximate solution. Setting $\alpha =1.8,a=b=0,M=8,T=1.$ the subplot (b) represents their absolute error at $M=8$.

Example 2. Consider the following general problem

$$\begin{array}{l}{\mathrm{D}}^{\alpha }Y+A{\mathrm{D}}^{{\nu }_1}Y+B{\mathrm{D}}^{{\nu }_2}Y+\mathit{CY}=f\left(x\right),1<\alpha \le 2,0<{\nu }_1,{\nu }_2\le 1,\\ \mathit{\text{subject to the boundary conditions}}Y\left(0\right)=1,Y\left(1\right)=5.\end{array}$$E53

Let the exact solution $Y\left(x\right)=x^4+x^3+x^2+x+1,$ and $A=1,B=2,C=3,$ the source terms

$$\begin{array}{l}f\left(x\right)=3x^4+3x^3+15x^2+9x+5+\frac{24}{\mathrm{\Gamma }\frac{21}{5}}{x}^{\frac{16}{5}}+\frac{6}{\mathrm{\Gamma }\frac{16}{5}}{x}^{\frac{11}{5}}+\frac{2}{\mathrm{\Gamma }\frac{11}{5}}{x}^{\frac{6}{5}}+\frac{1}{\mathrm{\Gamma }\frac{6}{5}}{x}^{\frac{1}{5}}\\ +\frac{48}{\mathrm{\Gamma }\frac{23}{5}}{x}^{\frac{18}{5}}+\frac{12}{\mathrm{\Gamma }\frac{18}{5}}{x}^{\frac{13}{5}}+\frac{4}{\mathrm{\Gamma }\frac{13}{5}}{x}^{\frac{8}{5}}+\frac{2}{\mathrm{\Gamma }\frac{8}{5}}{x}^{\frac{3}{5}}.\end{array}$$E54

We approximate the solutions of this problem with proposed method by setting $\alpha =2,{\nu }_1=0.8,{\nu }_2=0.4,a=b=1,M=8,T=1.$ and observed that the scheme provides much more accurate approximations for very small-scale level even at scale level $M=8$ maximum absolute error is $9\times {10}^{-17}$ as shown in Figure 3.

Figure 3.

Subplot (a) represents comparison of approximate and exact $Y\left(x\right)$ of example 2. Blue curve represents exact solution while the red square dots represent approximate solution. The subplot (b) represents their absolute error at $M=8$.

The absolute error for different choices of $a,b$ and $M$ and for $\alpha =2,{\nu }_1=0.8,{\nu }_2=0.4$ is given in Table 3.

$M$	$a$	$b$	$\alpha$	$\mid y_{\mathit{app}}-y_{\mathit{\text{exact}}}\mid$
$4$	$0$	$0$	$1.7$	$4\times {10}^{-3}$
$6$			$1.8$	$4\times {10}^{-5}$
	$1$	$1$	$1.9$	$1.5\times {10}^{-4}$
			$2$	$3.5\times {10}^{-13}$
$6$	$0$	$0$	$2$	$1.8\times {10}^{-13}$

Table 3.

Maximum absolute error at various values of $\alpha ,a,b$ and $M$ for example, (1).

From Table 4, we observed that as the orders of derivatives tend to their corresponding integer values, the approximate solutions tend to its exact value, which demonstrate the accuracy of numerical solutions obtained in our proposed method.

$M$	$a$	$b$	$\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$
$4$	$0$	$0$	$1\times {10}^{-3}$
$6$	$0.5$	$0.5$	$1.6\times {10}^{-15}$
$8$	$1$	$1$	$9\times {10}^{-17}$
$8$	$0.5$	$-0.5$	$2.5\times {10}^{-15}$
$10$	$0$	$2$	$4\times {10}^{-16}$

Table 4.

Maximum absolute error at $\alpha =2,{\nu }_1=0.8,{\nu }_2=0.4$ for example, (2).

In the following examples, we solve some coupled systems under the given conditions by our proposed method.

Example 3. Consider the coupled system given by

$$\begin{array}{l}{\mathrm{D}}^{\alpha }X+A_1{\mathrm{D}}^{{\nu }_1}X+B_1{\mathrm{D}}^{{\nu }_2}Y+C_{1}X+D_{1}Y=f\left(x\right),1<\alpha \le 2,0<{\nu }_1,{\nu }_2\le 1,\\ {\mathrm{D}}^{\beta }Y+A_2{\mathrm{D}}^{{\omega }_1}X+B_2{\mathrm{D}}^{{\omega }_2}Y+C_{2}X+D_{2}Y=g\left(x\right),1<\beta \le 2,0<{\omega }_1,{\omega }_2\le 1,\end{array}$$E55

subject to the boundary conditions

$$\begin{array}{l}X\left(0\right)=0.5,X\left(1\right)=0.5\\ Y\left(0\right)=0.6,Y\left(1\right)=0.6.\end{array}$$E56

We solve this problem under the given parameters

$$\begin{array}{l}{A}_1=2,A_2=2,B_1=2,B_2=2,C_1=3,C_2=4,D_1=5,D_2=5,\\ \alpha =\beta =2,{\nu }_1={\omega }_2=1,{\nu }_1={\omega }_1=1,a=b=0.\end{array}$$

Let the exact solutions at $\alpha =\beta =2,{\nu }_i={\omega }_i=1\left(i=1,2\right)$ are given by

$$\begin{array}{l}X\left(x\right)=x^2{\left(T-x\right)}^2+0.5,Y\left(x\right)=x{\left(T-x\right)}^3+0.6,\end{array}$$E57

the source terms are given by

$$\begin{array}{l}f\left(x\right)=4x\left(2x-2\right)-6x{\left(x-1\right)}^2-5x{\left(x-1\right)}^3+2{\left(x-1\right)}^2-2{\left(x-1\right)}^3+3x^2{\left(x-1\right)}^2\\ +\frac{\left(3172393274221175x^{\left(\frac{7}{5}\right)}\left(50x^2-85x+34\right)\right)}{33495522228568064}+2x^2+\frac{9}{2}\end{array}$$E58

$$\begin{array}{c}g\left(x\right)=2x^2\left(2x-2\right)-2x{\left(x-1\right)}^2-5x{\left(x-1\right)}^3-6{\left(x-1\right)}^2-2{\left(x-1\right)}^3-3x\left(2x-2\right)\\ +3x^2{\left(x-1\right)}^2+\frac{9}{2}\end{array}$$E59

We approximate the coupled systems at scale level $M=6$ for the given parameters. We see from Figures 4 and 5 that our method works very well, and the absolute error is about $3\times {10}^{-16}$ for $X\left(x\right)$ and for $Y\left(x\right)$ this amount is $7\times {10}^{-16}.$ Which is very small numbers, which prove the applicability of our methods. Absolute error for various $a,b,M$ is given in next Table 5.

Figure 4.

Subplot (a) represents comparison of approximate and exact $X\left(x\right)$ of example 3. Blue curve represents exact solution while the red circles represent approximate solution. The subplot (b) represents their absolute error at $M=6$.

Figure 5.

Subplot (c) represents comparison of approximate and exact $Y\left(x\right)$ of example 3. Blue dots represent exact solution while the red curve represents approximate solution. The subplot (d) represents their absolute error at $M=6$.

$M$	$a$	$b$	$\mid X_{\mathit{app}}-X_{\mathit{\text{exact}}}\mid$	$\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$
$4$	$0$	$0$	$1\times {10}^{-3}$	$1\times {10}^{-3}$
$6$			$3\times {10}^{-16}$	$7\times {10}^{-16}$
$6$	$1$	$1$	$4.5\times {10}^{-16}$	$7\times {10}^{-16}$
$8$	$0.5$	$1$	$4.5\times {10}^{-16}$	$6\times {10}^{-16}$
$8$	$1$	$0.5$	$3.5\times {10}^{-16}$	$1.2\times {10}^{-15}$

Table 5.

Maximum absolute error at various values of $a,b$ and $M$ for example, (3).

Absolute error for specific $a,b$ and for various $\alpha ,\beta ,{\nu }_1,{\nu }_2,{\omega }_1,{\omega }_2$ is given in Table 6.

$M$	$\alpha$	$\beta$	${\nu }_1$	${\nu }_2$	${\omega }_1$	${\omega }_2$	$\mid X_{\mathit{app}}-X_{\mathit{\text{exact}}}\mid$	$\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$
$4$	$1.5$	$1.5$	$0.5$	$0.5$	$0.5$	$0.5$	$2.5\times {10}^{-3}$	$1\times {10}^{-3}$
$6$	$1.6$	$1.6$	$0.6$	$0.6$	$0.6$	$0.6$	$3\times {10}^{-3}$	$2\times {10}^{-3}$
$6$	$1.7$	$1.7$	$0.7$	$0.7$	$0.7$	$0.7$	$6\times {10}^{-4}$	$1.2\times {10}^{-3}$
$8$	$1.8$	$1.8$	$0.8$	$0.8$	$0.8$	$0.8$	$4.5\times {10}^{-4}$	$1\times {10}^{-4}$
$8$	$1.9$	$1.9$	$0.9$	$0.9$	$0.9$	$0.9$	$9\times {10}^{-5}$	$2\times {10}^{-5}$
$10$	$2$	$2$	$1$	$1$	$1$	$1$	$4.5\times {10}^{-16}$	$6\times {10}^{-16}$

Table 6.

Maximum absolute error at specific values of $a=b=0$ for example, (3).

From Table 5, we see the effect of $a,b$ at various scale level $M$. While Table 6 indicates the effect of various choices of $\alpha ,\beta ,{\nu }_1,{\nu }_2,{\omega }_1,{\omega }_2$ at different scale level. As the scale level increases, values of $\alpha ,\beta ,{\nu }_1,{\nu }_2,{\omega }_1,{\omega }_2$ tend to their integer values. The absolute error is reducing and the solutions are tending to their exact value, which demonstrate the applicability of the proposed method.

Example 4. Consider the mathematical models of fractionally damped coupled system of spring masses system whose model is given by

$$\begin{array}{l}{\mathrm{D}}^{\alpha }X+\frac{\lambda }{m_1}{\mathrm{D}}^{{\nu }_1}X-\frac{\mu }{m_1}{\mathrm{D}}^{{\nu }_2}Y+\frac{{\kappa }_1}{m_1}X-\frac{{\kappa }_2}{m_1}Y=f\left(x\right),1<\alpha \le 2,0<{\nu }_1,{\nu }_2\le 1,\\ {\mathrm{D}}^{\beta }Y-\frac{\lambda }{m_2}{\mathrm{D}}^{{\omega }_1}X+\frac{\lambda +\mu }{m_2}{\mathrm{D}}^{{\omega }_2}Y-\frac{{\kappa }_1}{m_2}X+\frac{{\kappa }_1+{\kappa }_2}{m_2}Y=g\left(x\right),1<\beta \le 2,0<{\omega }_1,{\omega }_2\le 1,\end{array}$$E60

subject to the boundary conditions

$$\begin{array}{l}X\left(0\right)=0,X\left(1\right)=0\\ Y\left(0\right)=1,Y\left(1\right)=-1.\end{array}$$E61

Where the source terms are given by

$$\begin{array}{l}f\left(x\right)=cos\left(\mathit{\pi x}\right)-sin\left(\mathit{\pi x}\right)-\pi cos\left(\mathit{\pi x}\right)-\pi sin\left(\mathit{\pi x}\right)-{\pi }^{2}sin\left(\mathit{\pi x}\right)\\ g\left(x\right)=sin\left(\mathit{\pi x}\right)-2cos\left(\mathit{\pi x}\right)+\pi cos\left(\mathit{\pi x}\right)+2\pi sin\left(\mathit{\pi x}\right)-{\pi }^{2}cos\left(\mathit{\pi x}\right).\end{array}$$E62

Where $\lambda ,\mu$ are damping parameters and ${\kappa }_1,{\kappa }_2$ are spring constants and $m_1,m_2$ are masses of objects and springs are hanged from both ends. We solve this problem under the following values:

$$\begin{array}{l}\lambda =2,\mu =2,{\kappa }_1=2,{\kappa }_2=2,m_1=m_2=2,\\ \mathit{\text{then calculating the coefficient of the problems}}\mathit{as}\\ A_1=1,A_2=-1,B_1=-1,B_2=2,C_1=1,C_2=-1,D_1=-1,D_2=1.\end{array}$$

Let the exact solutions at $\alpha =\beta =2,{\nu }_i={\omega }_i=1\left(i=1,2\right)$ are given by

$$X\left(x\right)=sin\left(\mathit{\pi x}\right),Y\left(x\right)=cos\left(\mathit{\pi x}\right)$$E63

The above model is approximated for the solutions by our proposed methods (Figures 6 and 7). We observed that our method provides best approximate solutions to the problems for small-scale level $M=5.$ We also find numerical solutions for fractional values of $\alpha ,\beta ,{\nu }_i,{\omega }_i\left(i=1,2\right)$. We observe that as these orders tend to their integer values, the solutions are tending to the exact, which prove the applicability of our method in the following tables.

Figure 6.

Subplot (a) represents comparison of approximate and exact $X\left(x\right)$ at $M=5$. The subplot (b) represents comparison of approximate and exact $Y\left(x\right)$ at $M=5$ for example, (4). Blue curve represents exact solution while the red curve represents approximate solution for $a=b=0$.

Figure 7.

Subplot (c) represents absolute error $\mid X_{\mathit{app}}-X_{\mathit{\text{exact}}}\mid$_, and subplot (d) represents their absolute error $\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$ at $M=6$ for example, (4) for $a=b=0$.

Absolute error for various $a,b,M$ and $\alpha =\beta =2,{\nu }_i={\omega }_i=1\left(i=1,2\right)$ is given in Table 7.

$M$	$a$	$b$	$\mid X_{\mathit{app}}-X_{\mathit{\text{exact}}}\mid$	$\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$
$4$	$0$	$0$	$1\times {10}^{-3}$	$1\times {10}^{-3}$
$5$			$10\times {10}^{-4}$	$6\times {10}^{-3}$
$6$			$9\times {10}^{-7}$	$9\times {10}^{-6}$
$6$	$1$	$1$	$4\times {10}^{-3}$	$7\times {10}^{-3}$
$8$	$0.5$	$.5$	$2\times {10}^{-3}$	$1\times {10}^{-3}$
$8$	$1$	$0.5$	$3.5\times {10}^{-3}$	$1.2\times {10}^{-3}$

Table 7.

Maximum absolute error at various values of $a,b$ and $M$ for example, (4).

We observe from Table 7 that values of $a,b$ play important role in the approximate solution. By giving integral value to $a,b,$ methods give best approximate solution for the given problem. Similarly for $a=b=0$ and for various choices of $\alpha ,\beta ,{\nu }_1,{\nu }_2,{\omega }_1,{\omega }_2$_, the absolute errors are given in Table 8.

$M$	$\alpha$	$\beta$	${\nu }_1$	${\nu }_2$	${\omega }_1$	${\omega }_2$	$\mid X_{\mathit{app}}-X_{\mathit{\text{exact}}}\mid$	$\mid Y_{\mathit{app}}-Y_{\mathit{\text{exact}}}\mid$
$4$	$1.5$	$1.5$	$0.5$	$0.5$	$0.5$	$0.5$	$2.5\times {10}^{-2}$	$1\times {10}^{-2}$
$6$	$1.6$	$1.6$	$0.6$	$0.6$	$0.6$	$0.6$	$4\times {10}^{-3}$	$2\times {10}^{-3}$
$6$	$1.7$	$1.7$	$0.7$	$0.7$	$0.7$	$0.7$	$1\times {10}^{-3}$	$1\times {10}^{-3}$
$8$	$1.8$	$1.8$	$0.8$	$0.8$	$0.8$	$0.8$	$4.5\times {10}^{-4}$	$1\times {10}^{-4}$
$8$	$1.9$	$1.9$	$0.9$	$0.9$	$0.9$	$0.9$	$1.2\times {10}^{-5}$	$2\times {10}^{-5}$
$10$	$2$	$2$	$1$	$1$	$1$	$1$	$9\times {10}^{-7}$	$8\times {10}^{-6}$

Table 8.

Maximum absolute error at $a=b=0$ and for various $\alpha ,\beta ,{\nu }_1,{\nu }_2,{\omega }_1,{\omega }_2$ for example, (4).

From Table 8, it is obvious that when the orders of the derivatives approach to their integral values, the error is reducing and the approximate solutions are converging to the exact solutions. This phenomenon indicates that the proposed method is highly accurate.

6. Conclusion and discussion

In this article, we have studied shifted Jacobi polynomials. Based on these polynomials, we recalled some already existing matrices of fractional order derivative and integrations from the literature as well as we constructed a new operational matrix corresponding to boundary conditions. By means of these operational matrices, we converted the system of fractional differential equations to simple and easily soluble system of algebraic equations. There is no need of Tau-collocation method. The simple algebraic equations were easily solved, and the result were plotted. From the plot, we observe that our proposed method is highly accurate and can be applied to a variety of problems of fractional order ordinary as well as partial differential equations. We also compared our results to the exact solutions and observed that our method gave satisfactory results. The proposed method can easily and accurately can be applied to a variety of problems of applied mathematics, engineering, and physics, etc.

Acknowledgments

The authors Aziz Khan, Kamal Shah and Thabet Abdeljawad would like to thank Prince Sultan university for the support through TAS research lab.

Mathematics Subject Classification:

26A33; 34A08; 34K37.