Extension of Operational Matrix Technique for the Solution of Nonlinear System of Caputo Fractional Differential Equations Subjected to Integral Type Boundary Constrains

Abstract

We extend the operational matrices technique to design a spectral solution of nonlinear fractional differential equations (FDEs). The derivative is considered in the Caputo sense. The coupled system of two FDEs is considered, subjected to more generalized integral type conditions. The basis of our approach is the most simple orthogonal polynomials. Several new matrices are derived that have strong applications in the development of computational scheme. The scheme presented in this article is able to convert nonlinear coupled system of FDEs to an equivalent S-lvester type algebraic equation. The solution of the algebraic structure is constructed by converting the system into a complex Schur form. After conversion, the solution of the resultant triangular system is obtained and transformed back to construct the solution of algebraic structure. The solution of the matrix equation is used to construct the solution of the related nonlinear system of FDEs. The convergence of the proposed method is investigated analytically and verified experimentally through a wide variety of test problems.

1. Introduction

Fractional calculus has a long and rich history. Its wide range of applications made this field an active area of mathematical research. Frequently investigated properties of FDEs include existence and uniqueness problems, multiplicity of solutions, stability of solutions and analytical study of the analytical properties of solution. Parallel to these area of research, one of the most explored and interested areas of research in this field is the design of new numerical methods for finding the approximate solutions to problems of this category. Many scientists and mathematicians are trying to design efficient and reliable techniques to find possible estimates to solutions of FDEs or their coupled systems.

There are many analytical, semi-analytical, numerical, and spectral approximations of solution to FDEs and their coupled systems. Among others, one of the easiest method is the differential transformation method (DTM). In [1], DTM is successfully applied to solve simple nonlinear FDEs with simple initial conditions. In [2], DTM is designed to solve the fractional-order counterpart of Korteweg De Vries (KDV) equation. The method is further improved for the solution of fractional-order boundary value problems in [3]. Solutions to fractional-order boundary value problems are also attempted with analytical methods such as the homotopy perturbation method; see for example [4,5,6]. The Adomian decomposition method is also a powerful analytical method [7,8,9]. Spectral methods have gained the attention of scholars in recent decades. Compared to other methods, spectral methods are easy to design and implement [10,11,12,13,14,15,16,17,18,19,20]. The availability of a wide range of orthogonal polynomials makes this method more interesting. They have the ability to solve fractional order problems, whose solutions are difficult or sometimes impossible to obtain with other traditional methods. For new readers, we strongly recommend studying the results obtained in [21,22,23,24,25,26] for a clear understanding and developing a good base. However, to the best of our knowledge, the spectral method becomes difficult and sometimes fails to handle the situation when boundary conditions are given in more complicated forms such as local conditions, nonlocal m-point terminal conditions, integral type terminal conditions, and radiation boundary conditions. Such conditions have solid application in various problems of traveling waves, heat conduction and electromagnetism. One can find a good example of application of integral type boundary condition to heat conduction phenomena in a rod of fixed length in a recent article [27].

Nonlocal FDEs arise in mathematical modeling of various problems in physics, engineering, ecology, and biological sciences [28,29,30]. Some of the numerical investigations regarding FDEs with nonlocal constrains are discussed in [31,32,33,34,35]. Numerical approaches such as finite difference and radial base function also remain a focus of interest. Application of these methods to one-dimensional heat-like equations has been studied in [32,36,37,38]. Two-dimensional diffusion problems [33,39,40] and Laplace equations with integral constraints are explored in [31].

Keeping in view the increasing interest of mathematicians in fractional calculus, we are strongly motivated to design a new spectral approximation technique for complicated problems such as

$$\begin{array}{cc}\hfill & {}^c{D}^{{\sigma }_1}u\left(t\right)=f(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \\ \hfill & {}^c{D}^{{\sigma }_2}v\left(t\right)=g(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \\ \hfill & u\left(0\right)=u_0,u\left(\tau \right)=m_1{\int }_0^{{\omega }_1}s\left(t\right)u\left(t\right)dt,0<{\omega }_1\le \tau ,\hfill \\ \hfill & v\left(0\right)=v_0,v\left(\tau \right)=m_2{\int }_0^{{\omega }_2}r\left(t\right)v\left(t\right)dt,0<{\omega }_2\le \tau .\hfill \end{array}$$

(1)

where $0<{\gamma }_1,{\gamma }_2\le 1,$ $1<{\sigma }_1,{\sigma }_2\le 2,$ $t\in [0,\tau ].$ The scalar functions $u\left(t\right)$ and $v\left(t\right)$ are the solution to be determined. f and g are nonlinear functions of $u\left(t\right),$ $v\left(t\right)$ and its fractional order derivatives and is assumed to be in such a form that the solution of the problem exists.

We start our discussion by introducing some definitions and preliminary results.

2. Preliminaries

The following definitions and notations are important for our further analysis. More details and theoretical understanding of these results, see [41,42,43,44,45].

Definition 1.

The Riemann–Liouville $\alpha -$order integral of $\varphi \in (L^1[a,b],R)$ is defined by the following relation.

$$I^{\alpha }\varphi \left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_a^t{(t-s)}^{\alpha -1}\varphi \left(s\right)ds,$$

Definition 2.

For $\varphi \left(t\right)\in C^n[a,b]$, the α order Caputo derivative is defined as

$$\begin{array}{cc}\hfill {}^c{D}^{\alpha }\varphi \left(t\right)=& \frac{1}{\Gamma (n-\alpha )}{\int }_a^t\frac{{\varphi }^{\left(n\right)}\left(s\right)}{{(t-s)}^{\alpha +1-n}}ds,\hfill \\ \hfill & n-1\le \alpha <n,n\in N,\hfill \end{array}$$

where $n=\left[\alpha \right]+1$.

From the above definition, it is easy to extract $I^{\alpha }{t}^b=\frac{\Gamma (1+b)}{\Gamma (1+b+\alpha )}{t}^{b+\alpha }\mathrm{for}\alpha >0,k\ge 0$, ${}^c{D}^{\alpha }C=0,$ and

$${}^c{D}^{\alpha }{t}^b=\frac{\Gamma (1+b)}{\Gamma (1+b-\alpha )}{t}^{b-\alpha },\mathrm{for}b\ge \left[\alpha \right].$$

(2)

The Shifted Legendre Polynomials (LP)

These polynomials play a central role in approximation theory. Generally, these polynomials are defined on the domain $[0,\tau ]$, which is given by

$${\rho }_l^{\tau }\left(t\right)=\sum _{b=0}^l{\gimel }_{(l,b)}{t}^b,$$

(3)

where

$${\gimel }_{(l,b)}=\frac{{(-1)}^{l+b}(l+b)!}{(l-b)!{\tau }^b{(b!)}^2}.$$

(4)

These polynomials enjoys a very important property of orthogonality on the domain $[0,\tau ],$ which is expressed mathematically as

$${\int }_0^{\tau }{\rho }_i^{\tau }\left(t\right){\rho }_j^{\tau }\left(t\right)dt=\left\{\begin{array}{cc}\frac{\tau }{2j+1},\hfill & \mathrm{if}\hspace{0.17em}i=j,\hfill \\ 0,\hfill & \mathrm{if}\hspace{0.17em}i\ne j.\hfill \end{array}\right.$$

(5)

By using Equation (5), any $s\left(t\right)$ can be expressed in terms of LP as

$$s\left(t\right)\approx \sum _{l=0}^m{c}_l{\rho }_l^{\tau }\left(t\right),\mathrm{where}{c}_l=\frac{(2l+1)}{\tau }{\int }_0^{\tau }s\left(t\right){\rho }_l^{\tau }\left(t\right)dt.$$

(6)

The above equation has an equivalent vector representation given as

$$s\left(t\right)\approx C_M{\Lambda }_M^{\tau }\left(t\right),$$

(7)

where

$${\Lambda }_M^{\tau }\left(t\right)={\left[\begin{array}{cccccc}{\rho }_0^{\tau }\left(t\right)& {\rho }_1^{\tau }\left(t\right)& \cdots & {\rho }_i^{\tau }\left(t\right)& \cdots & {\rho }_m^{\tau }\left(t\right)\end{array}\right]}^T$$

(8)

and

$$C_M\hspace{0.17em}=\hspace{0.17em}\left[\begin{array}{cccccc}{c}_0& c_1& \cdots & c_i& \cdots & c_m\end{array}\right].$$

(9)

The following useful constant is important in the derivation of the operational matrices. We only recall the definition of the constant. The detailed derivation of which can be found in [25].

Lemma 1.

The integral of product of any three LP is a constant, represented by, ${\wp }^{(i,j,k)},$ defined as

$${\int }_0^{\tau }{\rho }_i^{\tau }\left(t\right){\rho }_j^{\tau }\left(t\right){\rho }_k^{\tau }\left(t\right)dt={\wp }^{(i,j,k)},$$

where

$${\wp }^{(i,j,k)}=\sum _{a=0}^i\sum _{p=0}^j\sum _{r=0}^k{\gimel }_{(i,a)}{\gimel }_{(i,p)}{\gimel }_{(i,r)}{\mathrm{Y}}_{(a,p,r)}.$$

${\gimel }_{(.,.)}$ are as defined in (4) and

$${\mathrm{Y}}_{(a,p,r)}=\frac{{\tau }^{(a+p+r+1)}}{(a+p+r+1)}.$$

The constant defined above is important in the solution of FDEs. We recall one more important result of the Legendre polynomials, which is their application in the study of convergence and developing of error bounds.

Theorem 1

([21]). Let ${\prod }_M$ be the space of M terms Legendre polynomials and let $u\left(t\right)\in C^m[0,1],$ then $u_m\left(t\right)$ is in space ${\prod }_M$; then, for the m term approximation,

$$u\left(t\right)=\sum _{i=0}^m{c}_i{\rho }_i^{\tau }\left(t\right),$$

there exists a constant C such that

$$c_k\simeq \frac{C}{{\lambda }_k}\left|u^{\left(m\right)}\right|.$$

and

$$|u\left(t\right)-\sum _{i=0}^m{c}_i{\rho }_i^{\tau }{\left(t\right)|}^2\le \sum _{k=m+1}^{\infty }{\lambda }_k{c}_k^2,$$

where

$$c_k=\frac{{\lambda }_k+1}{\tau }{\int }_0^{\tau }u\left(t\right){\rho }_k^{\tau }\left(t\right)dt,{\lambda }_k=k(k+1).$$

C is constant and can be chosen in such a way that $u^{\left(2m\right)}$ belongs to ${\prod }_M$, where $u^{\left(m\right)}$ is defined as

$$u^{\left(m\right)}=\mathbf{Z}\left(u^{(m-1)}\right)={\mathbf{Z}}^m\left(u^{\left(0\right)}\right)$$

where Z is storm livoliel operator of legendre polynomials with $u^{\left(0\right)}=u\left(t\right).$

In the next section, we first recall some of our previously designed operational matrices and then develop new operational matrices.

3. Operational Matrices (OP)

OP acts as a basic block in developing approximation techniques. The purpose of operational matrices is to replace a given derivative term with its matrix notation. The following matrices are important in our further investigation.

Lemma 2

([24]). Let ${\Lambda }_M^{\tau }\left(t\right)$ be the function vector; the α order integration is generalized as

$$I^{\alpha }\left({\Lambda }_M^{\tau }\left(t\right)\right)\simeq {\mathbb{H}}_{M\times M}^{\tau ,\alpha }{\Lambda }_M^{\tau }\left(t\right),$$

where ${\mathbb{H}}_{M\times M}^{\tau ,\alpha }$ is the OP of integration, defined as

$${\mathbb{H}}_{M\times M}^{\tau ,\alpha }\hspace{0.17em}=\hspace{0.17em}\left[\begin{array}{c}{\Theta }_{i,j,\tau }\end{array}\right],$$

where

$${\Theta }_{i,j,\tau }=\sum _{a=0}^i{s}_{(a,j)}{\gimel }_{(i,a)}\frac{\Gamma (a+1)}{\Gamma (a+\alpha +1)}.$$

where

$$s_{(a,j)}=\frac{(2j+1)}{\tau }\sum _{l=0}^j\frac{{(-1)}^{j+l}(j+l)!{\left(\tau \right)}^{a+l+\alpha +1}}{\left({\tau }^l\right)(j-l){(l!)}^2(a+l+\alpha +1)}.$$

Corollary 1.

The error $|{\mathbf{E}}_M|=|I^{\alpha }u\left(t\right)-C_M{\mathbb{H}}_{M\times M}^{\tau ,\alpha }{\Lambda }_M^{\tau }\left(t\right)|$ is bounded by the following relation

$$|{\mathbf{E}}_M|\le |\sum _{a=m+1}^{\infty }{c}_k\{\sum _{i=0}^m{\Theta }_{a,j,\tau }\}|.$$

Proof. 

Consider

$$u\left(t\right)=\sum _{a=0}^{\infty }{c}_a{\rho }_a^{\tau }\left(t\right).$$

Then, using the previous result, we get

$$I^{\alpha }u\left(t\right)=\sum _{a=0}^{\infty }{c}_a\sum _{j=0}^m{\Theta }_{a,j,\tau }{\rho }_j^{\tau }\left(t\right).$$

Simplifying the above estimate

$$I^{\alpha }u\left(t\right)-\sum _{a=0}^m{c}_a\sum _{j=0}^m{\Theta }_{a,j,\tau }{\rho }_j^{\tau }\left(t\right)=\sum _{a=m+1}^{\infty }{c}_a\sum _{j=0}^m{\Theta }_{a,j,\tau }{\rho }_j^{\tau }\left(t\right).$$

In matrix notation

$$I^{\alpha }u\left(t\right)-C_M{\mathbb{H}}_{M\times M}^{\tau ,\alpha }{\Lambda }_M^{\tau }\left(t\right)=\sum _{a=m+1}^{\infty }{c}_a\sum _{j=0}^m{\Theta }_{a,j,\tau }{\rho }_j^{\tau }\left(t\right).$$

Using the fact $|{\rho }_j^{\tau }\left(t\right)|\le 1$ for $t\in [0,\tau ]$, therefore, we can write

$$|I^{\alpha }u\left(t\right)-C_M{\mathbb{H}}_{M\times M}^{\tau ,\alpha }{\Lambda }_M^{\tau }\left(t\right)|\le |\sum _{a=m+1}^{\infty }{c}_a\sum _{j=0}^m{\Theta }_{a,j,\tau }|.$$

 □

Lemma 3

([24]). Let ${\Lambda }_M^{\tau }\left(t\right)$ be the function vector as defined in (8); then the derivative of order σ of ${\Lambda }_M^{\tau }\left(t\right)$ is generalized as

$${}^c{D}^{\sigma }\left({\Lambda }_M^{\tau }\left(t\right)\right)\simeq {\mathbb{G}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right),$$

where ${\mathbb{G}}_{M\times M}^{\tau ,\sigma }$ is the operational matrix of derivative of order σ, and ${\mathbb{G}}_{M\times M}^{\tau ,\sigma }$ is defined as

$${\mathbb{G}}_{M\times M}^{\tau ,\sigma }\hspace{0.17em}=\hspace{0.17em}\left[\begin{array}{c}{\Phi }_{i,j,\tau }\end{array}\right]$$

where

$${\Phi }_{i,j,\tau }=\sum _{k=\lceil \sigma \rceil }^i{s}_{(k,j)}{\gimel }_{(i,k)}\frac{\Gamma (k+1)}{\Gamma (k-\sigma +1)}$$

with ${\Phi }_{i,j,\tau }=0$ if $i<\lceil \sigma \rceil$.

Furthermore, ${\gimel }_{(i,k)}$ is similar as defined in (4) and

$$s_{(k,j)}=\frac{(2j+1)}{\tau }\sum _{l=0}^j\frac{{(-1)}^{j+l}(j+l)!{\left(\tau \right)}^{k+l-\sigma +1}}{\left({\tau }^l\right)(j-l){(l!)}^2(k+l-\sigma +1)}.$$

Corollary 2.

The error $|{\mathbf{E}}_M{|\hspace{0.17em}=\hspace{0.17em}|}^c{D}^{\sigma }u\left(t\right)-C_M{\mathbb{G}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right)|$ in approximating $D^{\sigma }u\left(t\right)$ with operational matrix of derivative is bounded by the following.

$$|{\mathbf{E}}_M|\le |\sum _{k=m+1}^{\infty }{u}_k\{\sum _{i=\lceil \sigma \rceil }^m{\Phi }_{i,j,\tau }\}|$$

Proof. 

The proof of this corollary is similar as Corollary 1. □

Lemma 4

([24]). Let $u\left(t\right)$ and ${\varphi }_n\left(t\right)$ be smooth functions that are well-defined on $[0,\tau ]$. Then

$${\varphi }_n\left(t\right){}^c{D}^{\sigma }u\left(t\right)={\mathbf{W}}_M{\mathbf{B}}_{{\varphi }_n}^{\sigma }{\Lambda }_M^{\tau }\left(t\right).$$

where ${\mathbf{W}}_M$ is the Legendre coefficients vector of $u\left(t\right)$ as defined in (7) and

$${\mathbf{B}}_{{\varphi }_n}^{\sigma }={\mathbb{G}}_{M\times M}^{\tau ,\sigma }{J}_{M\times M}^{\tau ,{\varphi }_n}.$$

$${\mathbf{B}}_{{\varphi }_n}^{\sigma }\hspace{0.17em}=\hspace{0.17em}\left[\begin{array}{c}\overbrace{{\Theta }_{(r,j)}}\end{array}\right]$$

where

$$\overbrace{{\Theta }_{(i,j)}}=\sum _{l=0}^m{\Phi }_{(i,l)}{\Omega }_{(l,j)}$$

The matrix ${\mathbb{G}}_{M\times M}^{\tau ,\sigma }$ is the operational matrix of derivative as defined in Lemma 3; the entries of matrix $J_{M\times M}^{\tau ,{\varphi }_n}$ are defined by the following relation

$${\Omega }_{(l,j)}=\frac{2j+1}{\tau }\sum _{i=0}^m{c}_i{\wp }_{}^{(i,l,j)}.$$

and $c_i={\int }_0^{\tau }{\varphi }_n\left(t\right){\rho }_i\left(t\right)dt$.

Corollary 3.

The error $|{\mathbf{E}}_M|\hspace{0.17em}=\hspace{0.17em}|{\varphi }_n{\left(t\right)}^c{D}^{\sigma }u\left(t\right)-C_M{\mathbf{B}}_{{\varphi }_n}^{\sigma }{\Lambda }_M^{\tau }\left(t\right)|$ in approximating ${\varphi }_n{\left(t\right)}^c{D}^{\sigma }u\left(t\right)$ with operational matrix of fractional derivative with variable coefficient is bounded by the following.

$$|{\mathbf{E}}_M|\le |\sum _{k=m+1}^{\infty }{c}_k\sum _{j=0}^m\overbrace{{\Theta }_{(k,j)}}|.$$

Proof. 

Consider

$$u\left(t\right)=\sum _{k=0}^{\infty }{c}_k{\rho }_k^{\tau }\left(t\right).$$

Then, using the relation above, we get

$${\varphi }_n{\left(t\right)}^c{D}^{\sigma }u\left(t\right)=\sum _{k=0}^{\infty }{c}_k\sum _{j=0}^m\overbrace{{\Theta }_{(k,j)}}{\rho }_j^{\tau }\left(t\right)$$

Truncating the sum and writing in modified form we get

$${\varphi }_n{\left(t\right)}^c{D}^{\sigma }u\left(t\right)-\sum _{k=0}^m{c}_k\sum _{j=0}^m\overbrace{{\Theta }_{(k,j)}}{\rho }_j^{\tau }\left(t\right)=\sum _{k=m+1}^{\infty }{c}_k\sum _{j=0}^m\overbrace{{\Theta }_{(k,j)}}{\rho }_j^{\tau }\left(t\right).$$

We can also write it in matrix form as

$${\varphi }_n{\left(t\right)}^c{D}^{\sigma }u\left(t\right)-C_M{\mathbf{B}}_{{\varphi }_n}^{\sigma }{\Lambda }_M^{\tau }\left(t\right)=\sum _{k=m+1}^{\infty }{c}_k\sum _{j=0}^m\overbrace{{\Theta }_{(k,j)}}{\rho }_j^{\tau }\left(t\right)$$

Using the fact ${\rho }_j^{\tau }\left(t\right)\le 1$ for $t\in [0,\tau ]$, therefore, we can write

$$|{\varphi }_n{\left(t\right)}^c{D}^{\sigma }u\left(t\right)-C_M{\mathbf{B}}_{{\varphi }_n}^{\sigma }{\Lambda }_M^{\tau }\left(t\right)|\hspace{0.17em}\le \hspace{0.17em}|\sum _{k=m+1}^{\infty }{c}_k\sum _{j=0}^m\overbrace{{\Theta }_{(k,j)}}|.$$

Furthermore, hence the proof is complete. □

The above matrices have been successfully applied to solve fractional-order differential equations (FDEs) under the effect of initial conditions. However these matrices do not have the ability to solve FDEs with integral types of boundary conditions. Therefore, the invention of new matrices that can easily handle integral types of boundary conditions are of basic importance. In the forthcoming discussion, we will design two new operational matrices having the ability to deal with integral type boundary conditions.

Lemma 5.

Let $s\left(t\right)$ be a known function well defined on $[0,\tau ]$ and ${\varphi }_c^n\left(t\right)=c{t}^n$ be a polynomial then, for a function vector ${\Lambda }_M^{\tau }\left(t\right)$ as defined in (8), the following result holds

$${\varphi }_c^n\left(t\right){\int }_0^{\tau }s\left(t\right){\Lambda }_M^{\tau }\left(t\right)dt={\mathbb{Q}}_{M\times M}^{c,n,s\left(t\right)}{\Lambda }_M^{\tau }\left(t\right),$$

where the matrix ${\mathbb{Q}}_{M\times M}^{c,n,s\left(t\right)}$ is given as

$${\mathbb{Q}}_{M\times M}^{c,n,s\left(t\right)}=\left(\begin{array}{cccc}{\Omega }_{(0,0)}& {\Omega }_{(0,1)}& \cdots & {\Omega }_{(0,m)}\\ {\Omega }_{(1,0)}& {\Omega }_{(1,1)}& \cdots & {\Omega }_{(1,m)}\\ ⋮& ⋮& \ddots & ⋮\\ {\Omega }_{(m,0)}& {\Omega }_{(m,1)}& \cdots & {\Omega }_{(m,m)}\end{array}\right).$$

where

$${\Omega }_{(i,j)}=\sum _{r=0}^j\frac{\tau c{d}_i{\tau }^{r+n+1}{\gimel }_{(j,r)}}{(2i+1)(r+n+1)},$$

$d_i$ is the Legendre coefficients of the function $s\left(t\right)$ and ${\gimel }_{(j,r)}$ is as defined in (4).

Proof. 

Let $s\left(t\right)$ be approximated with Legendre polynomials, as

$$s\left(t\right)=\sum _{l=0}^m{d}_l{\varphi }_c^n\left(t\right).$$

We can write the $i^{th}$ term of ${\varphi }_c^n\left(t\right){\int }_0^{\tau }s\left(t\right){\Lambda }_M^{\tau }\left(t\right)dt$ as

$$\begin{array}{cc}\hfill c{t}^n{\int }_0^{\tau }\sum _{l=0}^m{d}_l{\rho }_l^{\tau }\left(t\right){\rho }_i^{\tau }\left(t\right)dt& =c{t}^n\sum _{l=0}^m{d}_l{\int }_0^{\tau }{\rho }_l^{\tau }\left(t\right){\rho }_i^{\tau }\left(t\right)dt,\hfill \\ \hfill & =\frac{\tau c{d}_i{t}^n}{2i+1}.\hfill \end{array}$$

The polynomial $\frac{\tau c{d}_i{t}^n}{2i+1}$ can be expressed as a series of Legendre polynomials as

$$\frac{\tau c{d}_i{t}^n}{2i+1}=\sum _{j=0}^m{\Omega }_{(i,j)}{\rho }_j\left(t\right).$$

where the constant ${\Omega }_{(i,j)}$ is given by

$${\Omega }_{(i,j)}=\frac{\tau c{d}_i}{2i+1}{\int }_0^{\tau }{\rho }_j^{\tau }\left(t\right)t^{n}dt.$$

Using the definition of ${\rho }_j^{\tau }\left(t\right)$ and after simplification, we can write

$$\begin{array}{cc}\hfill {\Omega }_{(i,j)}& =\frac{\tau c{d}_i}{2i+1}\sum _{r=0}^j{\gimel }_{(j,r)}{\int }_0^{\tau }{t}^{r+n}dt,\hfill \\ \hfill & =\sum _{r=0}^j\frac{\tau c{d}_i{\tau }^{r+n+1}{\gimel }_{(j,r)}}{(2i+1)(r+n+1)}.\hfill \end{array}$$

Simulating the result for $i=0:M$ and $j=0:M$ completes the proof of the Lemma. □

Lemma 6.

Let ${\varphi }_c^n\left(n\right)=c{t}^n$ be a polynomial then for a function vector ${\Lambda }_M^{\tau }\left(t\right)$ as defined in (8); the following holds

$${\varphi }_c^n\left(t\right){\Lambda }_M^{\tau }\left(\tau \right)={\mathbb{R}}_{M\times M}^{c,n,\tau }{\Lambda }_M^{\tau }\left(t\right)$$

The matrix ${\mathbb{R}}_{M\times M}^{c,n,\tau }$ is defined as

$${\mathbb{R}}_{M\times M}^{c,n,\tau }=\left(\begin{array}{cccc}{\Im }_{(0,0)}& {\Im }_{(0,1)}& \cdots & {\Im }_{(0,m)}\\ {\Im }_{(1,0)}& {\Im }_{(1,1)}& \cdots & {\Im }_{(1,m)}\\ ⋮& ⋮& \ddots & ⋮\\ {\Im }_{(m,0)}& {\Im }_{(m,1)}& \cdots & {\Im }_{(m,m)}\end{array}\right).$$

Furthermore, the entries are defined as

$${\Im }_{(i,j)}=\sum _{k=0}^i\sum _{l=0}^j\frac{c{\tau }^{k+n+l+1}{\gimel }_{(i,k)}{\gimel }_{(j,l)}}{n+l+1}.$$

Proof. 

The general term of the equality can be written as

$${\varphi }_c^n\left(t\right){\rho }_i^{\tau }\left(\tau \right)=\sum _{k=0}^i{\gimel }_{(i,k)}{\tau }^{k}c{t}^n.$$

It can be approximated with Legendre polynomials

$$\sum _{k=0}^i{\gimel }_{(i,k)}{\tau }^{k}c{t}^n=\sum _{j=0}^m{\Im }_{(i,j)}{\rho }_j^{\tau }\left(t\right),$$

where

$${\Im }_{(i,j)}={\int }_0^{\tau }\sum _{k=0}^i{\gimel }_{(i,k)}{\tau }^{k}c{t}^n{\rho }_j^{\tau }\left(t\right)dt.$$

Using the definition of Legendre polynomials we can write

$${\Im }_{(i,j)}=\sum _{k=0}^i\sum _{l=0}^j\frac{c{\tau }^{k+n+l+1}{\gimel }_{(i,k)}{\gimel }_{(j,l)}}{n+l+1}.$$

which completes the proof of the lemma. □

4. Application of Operational Matrices

The operational matrix method for the solution of fractional differential equations is, in fact, a spectral method. The main aim of the spectral method is to convert a typical differential equation to system of easily solvable algebraic equations, which can be solved to obtain the solution in the series form of some orthogonal polynomials. The application of these methods to nonlinear differential equations results in a nonlinear system of algebraic equations, which can be solved using some iterative algorithms (the Newton–Raphson method is a frequently used method), see for example [46,47,48,49,50,51,52,53].

In this paper, we implement a different approach. We first use the Taylor series method to linearize the nonlinear part f and g to convert the nonlinear fractional differential equation into a recurrence relation of linear fractional differential equations with variable coefficients.

4.1. Linear FDEs with Variable Coefficients

We first consider the following coupled system of linear fractional differential equations with variable coefficients

$$\begin{array}{cc}\hfill & {}^c{D}^{{\sigma }_1}u\left(t\right)=c_0\left(t\right)u+c_1\left(t\right)u^{\left({\gamma }_1\right)}+c_2\left(t\right)v+c_3\left(t\right)v^{\left({\gamma }_1\right)}+h\left(t\right),\hfill \\ \hfill & {}^c{D}^{{\sigma }_2}v\left(t\right)=d_0\left(t\right)u+d_1\left(t\right)u^{\left({\gamma }_1\right)}+d_2\left(t\right)v+d_3\left(t\right)v^{\left({\gamma }_1\right)}+k\left(t\right),\hfill \\ \hfill & u\left(0\right)=u_0,u\left(\tau \right)=m_1{\int }_0^{\tau }s\left(t\right)u\left(t\right)dt,\hfill \\ \hfill & v\left(0\right)=v_0,v\left(\tau \right)=m_2{\int }_0^{\tau }r\left(t\right)v\left(t\right)dt,\hfill \end{array}$$

(10)

where $m_1$ and $m_2$ are some real constants. $0<{\gamma }_1,{\gamma }_2\le 1,$ $1<{\sigma }_1,{\sigma }_2\le 2,$ $t\in [0,\tau ]$ and $c_i\left(t\right),d_i\left(t\right),s\left(t\right),r\left(t\right),h\left(t\right)$ and $k\left(t\right)$ are bounded, continuous, and well-defined functions on the domain $[0,\tau ].$

Assume the solution of (10) in terms of shifted Legendre polynomials, such that the following hold

$${}^c{D}^{{\sigma }_1}u\left(t\right)={\mathbf{K}}_M{\Lambda }_M^{\tau }\left(t\right),{}^c{D}^{{\sigma }_2}v\left(t\right)={\mathbf{L}}_M{\Lambda }_M^{\tau }\left(t\right).$$

(11)

Applying fractional integral of order ${\sigma }_1$ on the first equation of (11) and making use of Lemma 2, we can write

$$u\left(t\right)={\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\Lambda }_M^{\tau }\left(t\right)+e_0+e_{1}t.$$

(12)

By the application of initial conditions, we get $e_0=u_0,$ and to get the value of $e_1$, we use the integral-type boundary conditions, and after simplification it follows that

$$\begin{array}{cc}\hfill {\mathbf{K}}_M& {\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\Lambda }_M^{\tau }\left(\tau \right)+u_0+e_1\tau =\hfill \\ & m_1{\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\int }_0^{\tau }s\left(t\right){\Lambda }_M^{\tau }\left(t\right)dt\hfill \\ & +m_1{\int }_0^{\tau }s\left(t\right)u_{0}dt+m_1{e}_1{\int }_0^{\tau }s\left(t\right)tdt,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & e_1=\frac{1}{s_1}({\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\Lambda }_M^{\tau }\left(\tau \right)\hfill \\ & -m_1{\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\int }_0^{\tau }s\left(t\right){\Lambda }_M^{\tau }\left(t\right)dt-s_0),\hfill \end{array}$$

where $s_1=(m_1{\int }_0^{\tau }s\left(t\right)tdt-\tau )$ and $s_0=m_1{\int }_0^{\tau }s\left(t\right)u_{0}dt-u_0.$

Now using the values of $e_0$ and $e_1$ in (12), we can write $u\left(t\right)$ as

$$\begin{array}{cc}\hfill & u\left(t\right)={\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\Lambda }_M^{\tau }\left(t\right)+u_0+\hfill \\ & \frac{t}{s_1}\left({\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\Lambda }_M^{\tau }\left(\tau \right)-m_1{\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\int }_0^{\tau }s\left(t\right){\Lambda }_M^{\tau }\left(t\right)dt-s_0\right),\hfill \end{array}$$

(13)

In view of Lemma 5 and Lemma 6, we can write Equation (13) as

$$\begin{array}{cc}\hfill u\left(t\right)=& {\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\Lambda }_M^{\tau }\left(t\right)+{\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\mathbb{R}}_{M\times M}^{\frac{1}{s_1},1,\tau }{\Lambda }_M^{\tau }\left(t\right)\hfill \\ & -{\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}{\mathbb{Q}}_{M\times M}^{\frac{m_1}{s_1},1,s\left(t\right)}{\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{F}}_{\mathbf{1}}}_M{\Lambda }_M^{\tau }\left(t\right).\hfill \end{array}$$

where $u_0-\frac{s_{0}t}{s_1}={\mathbf{F}}_{\mathbf{1}}{\Lambda }_M^{\tau }\left(t\right).$ For simplicity of notation, we can write the above equations as

$$u\left(t\right)={\mathbf{K}}_M{{\mathbf{E}}^{\left(\mathbf{1}\right)}}_{M\times M}{\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{F}}_{\mathbf{1}}}_M{\Lambda }_M^{\tau }\left(t\right).$$

(14)

where

$${{\mathbf{E}}^{\left(\mathbf{1}\right)}}_{M\times M}={\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_1}\left({\mathbf{I}}_{M\times M}+{\mathbb{R}}_{M\times M}^{\frac{1}{s_1},1,\tau }-{\mathbb{Q}}_{M\times M}^{\frac{m_1}{s_1},1,s\left(t\right)}\right)$$

Repeating the same procedure for $v\left(t\right)$, we can get analogous representation as

$$v\left(t\right)={\mathbf{L}}_M{{\mathbf{E}}^{\left(\mathbf{2}\right)}}_{M\times M}{\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{F}}_{\mathbf{2}}}_M{\Lambda }_M^{\tau }\left(t\right).$$

(15)

where

$${\mathbf{E}}_{M\times M}^{\left(\mathbf{2}\right)}={\mathbb{H}}_{M\times M}^{\tau ,{\sigma }_2}\left({\mathbf{I}}_{M\times M}+{\mathbb{R}}_{M\times M}^{\frac{1}{r_1},1,\tau }-{\mathbb{Q}}_{M\times M}^{\frac{m_2}{r_1},1,r\left(t\right)}\right)$$

Now, in view of (15), (14), (11), and Lemma 4, the equivalent matrix form for system (10) is given as

$$\begin{array}{cc}\hfill {\mathbf{K}}_M{\Lambda }_M^{\tau }\left(t\right)=& {\mathbf{K}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{1}\right)}\left({\mathbf{B}}_{M\times M}^{c_0,0}+{\mathbf{B}}_{M\times M}^{c_1,{\gamma }_1}\right){\Lambda }_M^{\tau }\left(t\right)+\hfill \\ & {{\mathbf{F}}_{\mathbf{1}}}_M\left({\mathbf{B}}_{M\times M}^{c_0,0}+{\mathbf{B}}_{M\times M}^{c_1,{\gamma }_1}\right){\Lambda }_M^{\tau }\left(t\right)+\hfill \\ \hfill & {\mathbf{L}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{2}\right)}\left({\mathbf{B}}_{M\times M}^{c_2,0}+{\mathbf{B}}_{M\times M}^{c_3,{\gamma }_2}\right){\Lambda }_M^{\tau }\left(t\right)+\hfill \\ & {{\mathbf{F}}_{\mathbf{2}}}_M\left({\mathbf{B}}_{M\times M}^{c_2,0}+{\mathbf{B}}_{M\times M}^{c_3,{\gamma }_2}\right){\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{Z}}_{\mathbf{1}}}_M{\Lambda }_M^{\tau }\left(t\right),\hfill \\ \hfill {\mathbf{L}}_M{\Lambda }_M^{\tau }\left(t\right)=& {\mathbf{K}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{1}\right)}\left({\mathbf{B}}_{M\times M}^{d_0,0}+{\mathbf{B}}_{M\times M}^{d_1,{\gamma }_1}\right){\Lambda }_M^{\tau }\left(t\right)+\hfill \\ & {{\mathbf{F}}_{\mathbf{1}}}_M\left({\mathbf{B}}_{M\times M}^{d_0,0}+{\mathbf{B}}_{M\times M}^{d_1,{\gamma }_1}\right){\Lambda }_M^{\tau }\left(t\right)+\hfill \\ \hfill & {\mathbf{L}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{2}\right)}\left({\mathbf{B}}_{M\times M}^{d_2,0}+{\mathbf{B}}_{M\times M}^{d_3,{\gamma }_2}\right){\Lambda }_M^{\tau }\left(t\right)+\hfill \\ & {{\mathbf{F}}_{\mathbf{2}}}_M\left({\mathbf{B}}_{M\times M}^{d_2,0}+{\mathbf{B}}_{M\times M}^{d_3,{\gamma }_2}\right){\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{Z}}_{\mathbf{2}}}_M{\Lambda }_M^{\tau }\left(t\right)\hfill \end{array}$$

(16)

Canceling out the common term and after simplification of notation, we can write

$$\begin{array}{cc}\hfill {\mathbf{K}}_M=& {\mathbf{K}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{1}\right)}\left({\mathbf{B}}_{M\times M}^{c_0,0}+{\mathbf{B}}_{M\times M}^{c_1,{\gamma }_1}\right)+\hfill \\ & {\mathbf{L}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{2}\right)}\left({\mathbf{B}}_{M\times M}^{c_2,0}+{\mathbf{B}}_{M\times M}^{c_3,{\gamma }_2}\right)+{{\mathbf{Y}}_{\mathbf{1}}}_M\hfill \\ \hfill {\mathbf{L}}_M=& {\mathbf{K}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{1}\right)}\left({\mathbf{B}}_{M\times M}^{d_0,0}+{\mathbf{B}}_{M\times M}^{d_1,{\gamma }_1}\right)+\hfill \\ & {\mathbf{L}}_M{\mathbf{E}}_{M\times M}^{\left(\mathbf{2}\right)}\left({\mathbf{B}}_{M\times M}^{d_2,0}+{\mathbf{B}}_{M\times M}^{d_3,{\gamma }_2}\right)+{{\mathbf{Y}}_{\mathbf{2}}}_M,\hfill \end{array}$$

(17)

where

$${{\mathbf{Y}}_{\mathbf{1}}}_M={{\mathbf{F}}_{\mathbf{2}}}_M\left({\mathbf{B}}_{M\times M}^{c_2,0}+{\mathbf{B}}_{M\times M}^{c_3,{\gamma }_2}\right)+$$

$${{\mathbf{F}}_{\mathbf{1}}}_M\left({\mathbf{B}}_{M\times M}^{c_0,0}+{\mathbf{B}}_{M\times M}^{c_1,{\gamma }_1}\right)+{{\mathbf{Z}}_{\mathbf{1}}}_M,$$

and

$${{\mathbf{Y}}_{\mathbf{2}}}_M={{\mathbf{F}}_{\mathbf{2}}}_M\left({\mathbf{B}}_{M\times M}^{d_2,0}+{\mathbf{B}}_{M\times M}^{d_3,{\gamma }_2}\right)+$$

$${{\mathbf{F}}_{\mathbf{1}}}_M\left({\mathbf{B}}_{M\times M}^{d_0,0}+{\mathbf{B}}_{M\times M}^{d_1,{\gamma }_1}\right)+{{\mathbf{Z}}_{\mathbf{2}}}_M.$$

The above equations can also be written as

$$\begin{array}{cc}\hfill & \left[{\mathbf{K}}_M\hspace{1em}{\mathbf{L}}_M\right]\hspace{0.17em}=\hspace{0.17em}\left[{\mathbf{K}}_M\hspace{1em}{\mathbf{L}}_M\right]\hfill \\ & \left[\begin{array}{cc}{\mathbf{E}}^{\left(\mathbf{1}\right)}\left(B^{c_0,0}+{\mathbf{B}}^{c_1,{\gamma }_1}\right)& {\mathbf{E}}^{\left(\mathbf{1}\right)}\left({\mathbf{B}}^{d_0,0}+{\mathbf{B}}^{d_1,{\gamma }_1}\right)\\ {\mathbf{E}}^{\left(\mathbf{2}\right)}\left({\mathbf{B}}^{c_2,0}+{\mathbf{B}}^{c_3,{\gamma }_2}\right)& {\mathbf{E}}^{\left(\mathbf{2}\right)}\left({\mathbf{B}}^{d_2,0}+{\mathbf{B}}^{d_3,{\gamma }_2}\right)\end{array}\right]\hfill \\ & +\left[\begin{array}{cc}{{\mathbf{Y}}_{\mathbf{1}}}_M& {{\mathbf{Y}}_{\mathbf{2}}}_M\end{array}\right]\hfill \end{array}$$

(18)

We see that (18) is a linear system of matrix equations. By solving (17), we will get the required coefficients vector ${\mathbf{K}}_M$ and ${\mathbf{L}}_M$, which can be used in (14) and (15) to get approximation to the solution of the main problem.

4.2. Nonlinear FDEs

Nonlinear FDEs cannot be directly solved using the OP method; however, combining it with the quasilinearization method makes it easy to recursively solve nonlinear FDEs. The procedure of this technique is as given as follows.

• Approximate the initial solution, the solution of the linear part, by the method presented in previous section and name it $u_0\left(t\right)$ and $v_0\left(t\right).$
• Linearize the nonlinear part at $u_0\left(t\right)$ and $v_0\left(t\right)$. This will convert the system of nonlinear FDEs into a system of linear FDEs that is easily solvable with the method devolved. Solve it and name the solution as $u_1\left(t\right)$ and $v_1\left(t\right)$.
• Repeat step 1.

Consider the following nonlinear FDEs.

$$\begin{array}{cc}\hfill & {}^c{D}^{{\sigma }_1}u\left(t\right)=f(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \\ \hfill & {}^c{D}^{{\sigma }_2}v\left(t\right)=g(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \\ \hfill & u\left(0\right)=a_0,u\left(\tau \right)={\int }_0^{{\omega }_1}s\left(t\right)u\left(t\right)dt,0<{\omega }_1\le \tau ,\hfill \\ \hfill & v\left(0\right)=b_0,v\left(\tau \right)={\int }_0^{{\omega }_2}r\left(t\right)v\left(t\right)dt,0<{\omega }_2\le \tau .\hfill \end{array}$$

(19)

separating the linear and nonlinear parts, we get

$$\begin{array}{cc}\hfill & {}^c{D}^{{\sigma }_1}u\left(t\right)=L_1(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)})+N_1(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \\ \hfill & {}^c{D}^{{\sigma }_2}v\left(t\right)=L_2(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)})+N_2(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \end{array}$$

(20)

First solve the linear part:

$$\begin{array}{cc}\hfill & {}^c{D}^{{\sigma }_1}u\left(t\right)=L_1(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \\ \hfill & {}^c{D}^{{\sigma }_2}v\left(t\right)=L_2(u,v,u^{\left({\gamma }_1\right)},v^{\left({\gamma }_2\right)}),\hfill \end{array}$$

(21)

Its solution is named $u_0\left(t\right)$ and $v_0\left(t\right)$. The next step is to linearize the nonlinear part.

$$\begin{array}{cc}\hfill {}^c{D}^{{\sigma }_1}{u}_1\left(t\right)& =L_1(u_1,v_1,u_1^{\left({\gamma }_1\right)},v_1^{\left({\gamma }_2\right)})+N_1(u_0,v_0,u_0^{\left({\gamma }_1\right)},v_0^{\left({\gamma }_2\right)})\hfill \\ & +(u_1-u_0)\frac{\partial N_1}{\partial u_0}+(v_1-v_0)\frac{\partial N_1}{\partial v_0}+\hfill \\ \hfill & (u_1^{\left({\gamma }_1\right)}-u_0^{\left({\gamma }_1\right)})\frac{\partial N_1}{\partial u_0^{\left({\gamma }_1\right)}}+(v_1^{\left({\gamma }_2\right)}-v_0^{\left({\gamma }_2\right)})\frac{\partial N_1}{\partial v_0^{\left({\gamma }_2\right)}},\hfill \\ \hfill {}^c{D}^{{\sigma }_2}{v}_1\left(t\right)& =L_2(u_1,v_1,u_1^{\left({\gamma }_1\right)},v_1^{\left({\gamma }_2\right)})+N_2(u_0,v_0,u_0^{\left({\gamma }_1\right)},v_0^{\left({\gamma }_2\right)})\hfill \\ & +(u_1-u_0)\frac{\partial N_2}{\partial u_0}+(v_1-v_0)\frac{\partial N_2}{\partial v_0}+\hfill \\ \hfill & (u_1^{\left({\gamma }_1\right)}-u_0^{\left({\gamma }_1\right)})\frac{\partial N_2}{\partial u_0^{\left({\gamma }_1\right)}}+(v_1^{\left({\gamma }_2\right)}-v_0^{\left({\gamma }_2\right)})\frac{\partial N_2}{\partial v_0^{\left({\gamma }_2\right)}}.\hfill \end{array}$$

(22)

We get a system of FDEs with variable coefficients. The whole process can be expressed as

$$\begin{array}{cc}\hfill {}^c{D}^{{\sigma }_1}{u}_{r+1}\left(t\right)& =L_1(u_{r+1},v_{r+1},u_{r+1}^{\left({\gamma }_1\right)},v_{r+1}^{\left({\gamma }_2\right)})+\hfill \\ & N_1(u_r,v_r,u_r^{\left({\gamma }_1\right)},v_r^{\left({\gamma }_2\right)})\hfill \\ & +(u_{r+1}-u_r)\frac{\partial N_1}{\partial u_r}+(v_{r+1}-v_r)\frac{\partial N_1}{\partial v_r}+\hfill \\ \hfill & (u_{r+1}^{\left({\gamma }_1\right)}-u_r^{\left({\gamma }_1\right)})\frac{\partial N_1}{\partial u_r^{\left({\gamma }_1\right)}}+(v_{r+1}^{\left({\gamma }_2\right)}-v_r^{\left({\gamma }_2\right)})\frac{\partial N_1}{\partial v_r^{\left({\gamma }_2\right)}}\hfill \\ \hfill {}^c{D}^{{\sigma }_2}{v}_{r+1}\left(t\right)& =L_2(u_{r+1},v_{r+1},u_{r+1}^{\left({\gamma }_1\right)},v_{r+1}^{\left({\gamma }_2\right)})+\hfill \\ & N_2(u_r,v_r,u_r^{\left({\gamma }_1\right)},v_r^{\left({\gamma }_2\right)})+(u_{r+1}-u_r)\frac{\partial N_2}{\partial u_r}\hfill \\ & +(v_{r+1}-v_r)\frac{\partial N_2}{\partial v_r}+(u_{r+1}^{\left({\gamma }_1\right)}-u_r^{\left({\gamma }_1\right)})\frac{\partial N_2}{\partial u_r^{\left({\gamma }_1\right)}}\hfill \\ & +(v_{r+1}^{\left({\gamma }_2\right)}-v_r^{\left({\gamma }_2\right)})\frac{\partial N_2}{\partial v_r^{\left({\gamma }_2\right)}}\hfill \end{array}$$

(23)

The boundary conditions can be written as

$$\begin{array}{cc}\hfill & u_{r+1}\left(0\right)=a_0,u_{r+1}\left(\tau \right)={\int }_0^{{\omega }_1}s\left(t\right)u_{r+1}\left(t\right)dt,0<{\omega }_1\le \tau ,\hfill \\ \hfill & v_{r+1}\left(0\right)=b_0,v_{r+1}\left(\tau \right)={\int }_0^{{\omega }_2}r\left(t\right)v_{r+1}\left(t\right)dt,0<{\omega }_2\le \tau .\hfill \end{array}$$

(24)

It can be easily noted that (23) is fractional differential equation with variable coefficients.

5. Error Bound of the Approximate Solution and Convergence

In this section, we calculate a upper bound for error of approximation of solution with the proposed method.

5.1. Error Bound for Single Differential Equation

Consider the following fractional differential equation.

$${}^c{D}^{\sigma }u\left(t\right)=c_0\left(t\right)u\left(t\right)+c_1\left(t\right)u{\left(t\right)}^{\left(\gamma \right)}+h\left(t\right),$$

(25)

subject to the following initial and boundary conditions

$$u\left(0\right)=u_0,u\left(\tau \right)=m_1{\int }_0^{\tau }s\left(t\right)u\left(t\right)dt.$$

Our aim is to derive an upper bound for the proposed method. We have to calculate $R_M$ defined as

$$R_M={|}^c{D}^{\sigma }u\left(t\right)-{\mathbf{K}}_M{\mathbb{G}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right)|.$$

(26)

The solution of the above problem can be written in terms of shifted Legendre series such that

$${}^c{D}^{\sigma }u\left(t\right)=\sum _{k=0}^{\infty }{u}_k{\rho }_k^{\tau }\left(t\right)={\mathbf{K}}_M{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k{\rho }_k^{\tau }\left(t\right).$$

(27)

Applying a fractional integral of order $\sigma$, using an operational matrix of integration and using Corollary (1), we can write

$$u\left(t\right)-c_0-c_{1}t={\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{c}_k\sum _{i=0}^m{\Theta }_{i,k,\tau }{\rho }_i^{\tau }\left(t\right).$$

(28)

Which can be simplified as

$$u\left(t\right)={\mathbf{K}}_M{\mathbb{H}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right)+c_0+c_{1}t+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m{\Theta }_{i,k,\tau }{\rho }_i^{\tau }\left(t\right),$$

(29)

We know from (14) and the integral type boundary conditions that we can conclude,

$$u\left(t\right)={\mathbf{K}}_M{\mathbf{E}}_{M\times M}{\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{F}}_{\mathbf{1}}}_M{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m{\Theta }_{i,k,\tau }{\rho }_i^{\tau }\left(t\right),$$

Assume that $\widehat{X}={\mathbf{K}}_M{\mathbf{E}}_{M\times M}+{\mathbf{F}}_{\mathbf{1}}.$ Using Lemma (3) and Corollary (2), we can write

$$c_l\left(t\right)u{\left(t\right)}^{\left(\gamma \right)}={\widehat{X}}_M^T{B}_{c_l}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(l\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right).$$

(30)

Approximating $h\left(t\right)={\mathbf{F}}_{\mathbf{2}}{\Lambda }_M^{\tau }\left(t\right)+{\sum }_{k^{\prime }=m+1}^{\infty }{f}_{k^{\prime }}{\rho }_{k^{\prime }}^{\tau }\left(t\right)$ and using (30) in (25) we get

$$\begin{array}{cc}\hfill & {\mathbf{K}}_M{\Lambda }_M^{\tau }\left(t\right)-{\widehat{X}}_M^T{B}_{c_0}^0{\Lambda }_M^{\tau }\left(t\right)-{\widehat{X}}_M^T{B}_{c_1}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)-{\mathbf{F}}_{\mathbf{2}}{\Lambda }_M^{\tau }\left(t\right)=R_M\left(t\right).\hfill \end{array}$$

where $R_M\left(t\right)$ is defined by relation

$$\begin{array}{cc}\hfill R_M\left(t\right)=& \sum _{k=m+1}^{\infty }{u}_k{\rho }_k^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\hfill \\ & \sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k^{\prime }=m+1}^{\infty }{f}_{k^{\prime }}{\rho }_{k^{\prime }}^{\tau }\left(t\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill R_M\left(t\right)=& \sum _{k=m+1}^{\infty }{u}_k[{\rho }_k^{\tau }\left(t\right)+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\hfill \\ & \sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)]+\sum _{k^{\prime }=m+1}^{\infty }{f}_{k^{\prime }}{\rho }_{k^{\prime }}^{\tau }\left(t\right).\hfill \end{array}$$

Using the bounded property of Legendre polynomail, it follows that

$$\begin{array}{cc}\hfill |R_M\left(t\right)|\le & \sum _{k=m+1}^{\infty }|u_k\left|\right|[1+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}]|+\sum _{k^{\prime }=m+1}^{\infty }\left|f_{k^{\prime }}\right|.\hfill \end{array}$$

(31)

In view of Theorem (1), it is evident that if the function $u\left(t\right)$ and $f\left(t\right)$ are sufficiently smooth functions, then the sequence that defines their coefficient is convergent to zero. Hence, we conclude that as $m\to \infty$ the coefficients $u_m\to 0$ and $f_m\to 0$. Hence it can be easily observed that the error $|R_M\left(t\right)|\to 0.$ Equation (31) also establishes an upper bound of the error between the exact and approximate solution.

5.2. Error Bound for Coupled System of Fractional Differential Equations

Consider the following system of FDEs.

$$\begin{array}{cc}\hfill & {}^c{D}^{\sigma }u\left(t\right)=c_0\left(t\right)u\left(t\right)+c_1\left(t\right)v\left(t\right)+c_2\left(t\right)u{\left(t\right)}^{\left(\gamma \right)}+c_3\left(t\right)v{\left(t\right)}^{\left(\gamma \right)}+h\left(t\right),\hfill \\ \hfill & {}^c{D}^{\sigma }v\left(t\right)=d_0\left(t\right)u\left(t\right)+d_1\left(t\right)v\left(t\right)+d_2\left(t\right)u{\left(t\right)}^{\left(\gamma \right)}+d_3\left(t\right)v{\left(t\right)}^{\left(\gamma \right)}+g\left(t\right),\hfill \end{array}$$

(32)

subject to the following initial and boundary conditions

$$\begin{array}{cc}\hfill & u\left(0\right)=u_0,u\left(\tau \right)=m_1{\int }_0^{\tau }s\left(t\right)u\left(t\right)dt.\hfill \\ \hfill & v\left(0\right)=v_0,v\left(\tau \right)=m_2{\int }_0^{\tau }r\left(t\right)v\left(t\right)dt.\hfill \end{array}$$

Our aim is to derive an upper bound for the proposed method. We have to calculate $R_M^u$ and $R_M^v$, defined as

$$\begin{array}{cc}\hfill & R_M^u={|}^c{D}^{\sigma }u\left(t\right)-{\mathbf{K}}_M{\mathbb{G}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right)|,\hfill \\ \hfill & R_M^v={|}^c{D}^{\sigma }v\left(t\right)-{\mathbf{L}}_M{\mathbb{G}}_{M\times M}^{\tau ,\sigma }{\Lambda }_M^{\tau }\left(t\right)|.\hfill \end{array}$$

(33)

We know from (14) that

$$\begin{array}{cc}\hfill & u\left(t\right)={\mathbf{K}}_M{\mathbf{E}}_{M\times M}{\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{F}}_{\mathbf{1}}}_M{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m{\Theta }_{i,k,\tau }{\rho }_i^{\tau }\left(t\right),\hfill \\ \hfill & v\left(t\right)={\mathbf{L}}_M{\mathbf{E}}_{M\times M}{\Lambda }_M^{\tau }\left(t\right)+{{\mathbf{F}}_{\mathbf{2}}}_M{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{v}_k\sum _{i=0}^m{\Theta }_{i,k,\tau }{\rho }_i^{\tau }\left(t\right),\hfill \end{array}$$

Assume, $\widehat{X}={\mathbf{K}}_M{\mathbf{E}}_{M\times M}+{\mathbf{F}}_{\mathbf{1}}$ and $\widehat{Y}={\mathbf{L}}_M{\mathbf{E}}_{M\times M}+{\mathbf{F}}_{\mathbf{2}}.$ Using Lemma (3) and Corollary (2), we can write

$$\begin{array}{cc}\hfill & c_l\left(t\right)u{\left(t\right)}^{\left(\gamma \right)}={\widehat{X}}_M^T{B}_{c_l}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(l\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right),\hfill \\ \hfill & d_l\left(t\right)v{\left(t\right)}^{\left(\gamma \right)}={\widehat{Y}}_M^T{B}_{d_l}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{v}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(l\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right).\hfill \end{array}$$

(34)

Approximating $h\left(t\right)=D_1{\Lambda }_M^{\tau }\left(t\right)+{\sum }_{k=m+1}^{\infty }{d}_k{\rho }_k^{\tau }\left(t\right)$ and $g\left(t\right)=D_2{\Lambda }_M^{\tau }\left(t\right)+{\sum }_{k=m+1}^{\infty }$ $d_k^{\prime }{\rho }_k^{\tau }\left(t\right)$ and using (30) in (25) we get

$$\begin{array}{cc}\hfill & {\mathbf{K}}_M{\Lambda }_M^{\tau }\left(t\right)-{\widehat{X}}_M^T{B}_{c_0}^0{\Lambda }_M^{\tau }\left(t\right)-{\widehat{Y}}_M^T{B}_{c_1}^0{\Lambda }_M^{\tau }\left(t\right)-{\widehat{X}}_M^T{B}_{c_2}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)-{\widehat{Y}}_M^T{B}_{c_3}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)-D_1{\Lambda }_M^{\tau }\left(t\right)=R_M^u\left(t\right),\hfill \\ \hfill & {\mathbf{L}}_M{\Lambda }_M^{\tau }\left(t\right)-{\widehat{X}}_M^T{B}_{d_0}^0{\Lambda }_M^{\tau }\left(t\right)-{\widehat{Y}}_M^T{B}_{d_1}^0{\Lambda }_M^{\tau }\left(t\right)-{\widehat{X}}_M^T{B}_{d_2}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)-{\widehat{Y}}_M^T{B}_{d_3}^{\gamma }{\Lambda }_M^{\tau }\left(t\right)-D_2{\Lambda }_M^{\tau }\left(t\right)=R_M^v\left(t\right).\hfill \end{array}$$

where $R_M^u\left(t\right)$ and $R_M^v\left(t\right)$ is defined by the relation

$$\begin{array}{cc}\hfill R_M^u\left(t\right)=& \sum _{k=m+1}^{\infty }{u}_k{\rho }_k^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{v}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\hfill \\ & \sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{v}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{d}_k{\rho }_k^{\tau }\left(t\right),\hfill \\ \hfill R_M^v\left(t\right)=& \sum _{k=m+1}^{\infty }{v}_k{\rho }_k^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{v}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\hfill \\ & \sum _{k=m+1}^{\infty }{u}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{v}_k\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}{\rho }_{i^{\prime }}^{\tau }\left(t\right)+\sum _{k=m+1}^{\infty }{d}_k^{\prime }{\rho }_k^{\tau }\left(t\right),\hfill \end{array}$$

Using the bounded property of the Legendre polynomial, it follows that

$$\begin{array}{cc}\hfill |R_M^u\left(t\right)|\le & \sum _{k=m+1}^{\infty }|u_k\left|\right|[1+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}]|+\hfill \\ & \sum _{k=m+1}^{\infty }|v_k\left|\right|[\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}]|+\sum _{k=m+1}^{\infty }\left|d_k\right|,\hfill \\ \hfill |R_M^v\left(t\right)|\le & \sum _{k=m+1}^{\infty }|v_k\left|\right|[1+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}]|+\hfill \\ & \sum _{k=m+1}^{\infty }|u_k\left|\right|[\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(0\right)}}+\sum _{i=0}^m\sum _{i^{\prime }=0}^m{\Theta }_{i,k,\tau }\overbrace{{\Theta }_{k,i^{\prime },\tau }^{\left(1\right)}}]|+\sum _{k=m+1}^{\infty }\left|d_k^{\prime }\right|.\hfill \end{array}$$

(35)

The above equation establishes an error bound for the solution $u\left(t\right)$ and $v\left(t\right).$ It also ensures the convergence of the proposed method for the solution of coupled system of FDEs.

6. Test Problems

We solve one single equation, three systems of linear FDEs with variable equations, and two systems of nonlinear problems, and analyze the convergence of the approximate solution by measuring the following error norms.

$$\left|\right|E_u{\left|\right|}_2=\frac{1}{\tau }{\int }_0^{\tau }{U}_M\left(t\right)-u\left(t\right)dt$$

and

$$\left|\right|E_u{\left|\right|}_{max}=ma{x}_{x\in [0,\tau ]}|U_M\left(t\right)-u\left(t\right)|.$$

We check the accuracy of the boundary condition by measuring the following error norms.

$$\left|\right|E_u{\left|\right|}_b=|U_M\left(\tau \right)-m_1{\int }_0^{\tau }{U}_M\left(t\right)|.$$

In the above bounds, the quantity $U_M\left(t\right)$ represents the $m-$ term approximation to the solution $u\left(t\right)$.

Test Problem 1.

$$\begin{array}{cc}\hfill {}^c{D}^{1.2}u\left(t\right)=(t^2+sin\left(t\right))u\left(t\right)& +(2t-t^3)u^{0.7}\left(t\right)+f\left(t\right),\hfill \\ \hfill u\left(0\right)=4,& u\left(1\right)=1.1216{\int }_0^{1}cos\left(t\right)u\left(t\right)dt,\hfill \end{array}$$

where the exact solution $u\left(t\right)=t^3+t^2+t+4$, and the source term $f\left(t\right)=\frac{38683084397149375t^{\frac{4}{5}}\left(50t^2-35t+21\right)}{378302368699121664}-\left(sin\left(t\right)+t^2\right)\left(t^4-t^3+t^2+4\right)-\frac{150543064388819875t^{\frac{13}{10}}\left(2t-t^3\right)\left(400t^2-330t+253\right)}{22218508761632342016}.$

Test Problem 2.

$$\begin{array}{cc}\hfill {}^c{D}^{\sigma }u\left(t\right)=(t+1)u\left(t\right)& +(1-t)u^{{\sigma }_1}\left(t\right)+\left(2t\right)v\left(t\right)\hfill \\ & +\left(t^2\right)v^{\prime }\left(t\right)+f\left(t\right),\hfill \\ \hfill {}^c{D}^{\sigma }v\left(t\right)=(t^2+1)u\left(t\right)& +(1-t^2)u^{{\sigma }_1}\left(t\right)+\left(3t\right)v\left(t\right)\hfill \\ & +\left(t^3\right)v^{\prime }\left(t\right)+g\left(t\right),\hfill \\ \hfill u\left(0\right)=1,& u\left(1\right)=2.1270{\int }_0^{1}sin\left(t\right)u\left(t\right)dt,\hfill \\ \hfill v\left(0\right)=-1,& v\left(1\right)=-1.8925{\int }_0^{1}cos\left(t\right)v\left(t\right)dt.\hfill \end{array}$$

where the exact solution $u\left(t\right)=t^2+t^3+e^t$ and $v\left(t\right)=t^2-t^3-e^t$, and the source term $f\left(t\right)=6t+{\mathrm{e}}^t+2t\left({\mathrm{e}}^t-t^2+t^3\right)-\left(t+1\right)\left({\mathrm{e}}^t+t^2+t^3\right)+\left(t-1\right)\left(2t+{\mathrm{e}}^t+3t^2\right)+t^2\left({\mathrm{e}}^t-2t+3t^2\right)+2$ and $g\left(t\right)=3t\left({\mathrm{e}}^t-t^2+t^3\right)-{\mathrm{e}}^t-6t+t^3\left({\mathrm{e}}^t-2t+3t^2\right)-\left(t^2+1\right)\left({\mathrm{e}}^t+t^2+t^3\right)+\left(t^2-1\right)\left(2t+{\mathrm{e}}^t+3t^2\right)+2$.

Test Problem 3.

$$\begin{array}{cc}\hfill {}^c{D}^{\sigma }u\left(t\right)=& sin\left(t\right)u\left(t\right)+cos\left(t\right)u^{{\sigma }_1}\left(t\right)+\hfill \\ & (sin\left(t\right)+cos\left(t\right))v\left(t\right)+sin\left(2t\right)v^{\prime }\left(t\right)+f\left(t\right),\hfill \\ \hfill {}^c{D}^{\sigma }v\left(t\right)=& \left(cos\left(2t\right)\right))u\left(t\right)+\left(tsin\left(t\right)\right)u^{{\sigma }_1}\left(t\right)\hfill \\ & +\left(tcos\left(t\right)\right)v\left(t\right)+\left(t^{2}sin\left(t\right)\right)v^{\prime }\left(t\right)+g\left(t\right),\hfill \\ \hfill u\left(0\right)=0,& u\left(1\right)=0.8{\int }_0^1{e}^{2t}sin\left(t\right)u\left(t\right)dt,\hfill \\ \hfill v\left(0\right)=1,& v\left(1\right)=0.71{\int }_0^1{e}^{-2t}cos\left(t\right)v\left(t\right)dt.\hfill \end{array}$$

where the exact solution $u\left(t\right)=e^{t}sin\left(t\right)$ and $v\left(t\right)=e^{-t}cos\left(t\right)$, and the source term $f\left(t\right)=2{\mathrm{e}}^{t}cos\left(t\right)-{\mathrm{e}}^t{sin\left(t\right)}^2-cos\left(t\right)\left({\mathrm{e}}^{t}cos\left(t\right)+{\mathrm{e}}^{t}sin\left(t\right)\right)+sin\left(2t\right)\left({\mathrm{e}}^{-t}cos\left(t\right)+{\mathrm{e}}^{-t}sin\left(t\right)\right)-{\mathrm{e}}^{-t}cos\left(t\right)\left(cos\left(t\right)+sin\left(t\right)\right)$ and $g\left(t\right)=2{\mathrm{e}}^{-t}sin\left(t\right)-t{\mathrm{e}}^{-t}{cos\left(t\right)}^2-tsin\left(t\right)\left({\mathrm{e}}^{t}cos\left(t\right)+{\mathrm{e}}^{t}sin\left(t\right)\right)-cos\left(2t\right){\mathrm{e}}^{t}sin\left(t\right)+t^{2}sin\left(t\right)\left({\mathrm{e}}^{-t}cos\left(t\right)+{\mathrm{e}}^{-t}sin\left(t\right)\right).$

Test Problem 4.

$$\begin{array}{cc}\hfill {}^c{D}^{\sigma }u\left(t\right)=e^{\left(t\right)}u\left(t\right)& +e^{(-t)}{u}^{{\sigma }_1}\left(t\right)+\hfill \\ & (e^{\left(t\right)}+e^{(-t)})v\left(t\right)+e^{\left(2t\right)}{v}^{\prime }\left(t\right)+f\left(t\right),\hfill \\ \hfill {}^c{D}^{\sigma }v\left(t\right)=(e^{(-2t)})& u\left(t\right)+(t{e}^{\left(t\right)})u^{{\sigma }_1}\left(t\right)\hfill \\ & +(t{e}^{\left(t\right)})v\left(t\right)+(t^2{e}^{\left(t\right)})v^{\prime }\left(t\right)+g\left(t\right),\hfill \\ \hfill u\left(0\right)=1,& u\left(1\right)=0.5{\int }_0^1(1-t)e^{2t}u\left(t\right)dt,\hfill \\ \hfill v\left(0\right)=1,& v\left(1\right)=10{\int }_0^1(1-t^2)e^{-2t}v\left(t\right)dt.\hfill \end{array}$$

where the exact solution $u\left(t\right)=t^4{(2-t)}^3$ and $v\left(t\right)=\left(t^3\right)(3-t)$, and the source term $f\left(t\right)={\mathrm{e}}^{-t}\left(4t^3{\left(t-2\right)}^3+3t^4{\left(t-2\right)}^2\right)-12t^2{\left(t-2\right)}^3-24t^3{\left(t-2\right)}^2-{\mathrm{e}}^{2t}\left(t^3\left(2t-6\right)+3t^2{\left(t-3\right)}^2\right)-3t^4\left(2t-4\right)+t^4{\mathrm{e}}^t{\left(t-2\right)}^3-t^3\left({\mathrm{e}}^{-t}+{\mathrm{e}}^t\right){\left(t-3\right)}^2$ and $g\left(t\right)=6t{\left(t-3\right)}^2+6t^2\left(2t-6\right)+2t^3-t^4{\mathrm{e}}^t{\left(t-3\right)}^2-t^{2}sin\left(t\right)\left(t^3\left(2t-6\right)+3t^2{\left(t-3\right)}^2\right)+t^4{\mathrm{e}}^{-2t}{\left(t-2\right)}^3+t{\mathrm{e}}^t(4t^3{\left(t-2\right)}^3+3t^4{\left(t-2\right)}^2).$

Test Problem 5.

$$\begin{array}{cc}\hfill & {}^c{D}^{\sigma }u\left(t\right)=u\left(t\right)+v\left(t\right)+u^2\left(t\right)+v^{\prime }\left(t\right)u^{\prime }\left(t\right)+f\left(t\right),\hfill \\ \hfill & {}^c{D}^{\sigma }v\left(t\right)=u^{\prime }\left(t\right)+v^{\prime }\left(t\right)+v{\left(t\right)}^2+u\left(t\right)u^{\prime }\left(t\right)+g\left(t\right),\hfill \\ \hfill & u\left(0\right)=0,u\left(1\right)=5.4069{\int }_0^{1}sin\left(t\right)u\left(t\right)dt,\hfill \\ \hfill & v\left(0\right)=1,v\left(1\right)=0.\hfill \end{array}$$

where the exact solution $u\left(t\right)=(t^2+1)t^2$ and $v\left(t\right)=t^3(1-t)$, and the source term $f\left(t\right)=t^3\left(t-1\right)+\left(3t^2\left(t-1\right)+t^3\right)\left(2t\left(t^2+1\right)+2t^3\right)-t^2\left(t^2+1\right)-t^4{\left(t^2+1\right)}^2+12t^2+2$ and $g\left(t\right)=3t^2\left(t-1\right)-2t\left(t^2+1\right)-t^6{\left(t-1\right)}^2-6t\left(t-1\right)-6t^2-t^3-t^2\left(t^2+1\right)\left(2t\left(t^2+1\right)+2t^3\right).$

Test Problem 6.

$$\begin{array}{cc}\hfill {}^c{D}^{\sigma }u\left(t\right)=& (t+1)u\left(t\right)+(1-t)u^{\prime }\left(t\right)+2tv\left(t\right)\hfill \\ & +t^2{v}^{\prime }\left(t\right)+u{\left(t\right)}^2+v^{\prime }\left(t\right)u^{\prime }\left(t\right)-v^3\left(t\right)+f\left(t\right),\hfill \\ \hfill {}^c{D}^{\sigma }v\left(t\right)=& (t^2+1)u\left(t\right)+(1-t^2)u^{\prime }\left(t\right)+3tv\left(t\right)\hfill \\ & +t^3{v}^{\prime }\left(t\right)+v{\left(t\right)}^2+u\left(t\right)u^{\prime }\left(t\right)-u^4\left(t\right)+g\left(t\right),\hfill \\ \hfill u\left(0\right)=0,& u\left(1\right)=0.6173{\int }_0^1(2t+1)u\left(t\right)dt,\hfill \\ \hfill v\left(0\right)=1,& v\left(1\right)=0.7311{\int }_0^1(3t+1)v\left(t\right)dt.\hfill \end{array}$$

where the exact solution $u\left(t\right)=sin\left(t\right)$ and $v\left(t\right)=e^t$, and the source term $f\left(t\right)={\mathrm{e}}^{3t}-sin\left(t\right)+cos\left(t\right)\left(t-1\right)-t^2{\mathrm{e}}^t-{\mathrm{e}}^{t}cos\left(t\right)-sin\left(t\right)\left(t+1\right)-{sin\left(t\right)}^2-2t{\mathrm{e}}^t$ and $g\left(t\right)={\mathrm{e}}^t-{\mathrm{e}}^{2t}-t^3{\mathrm{e}}^t-cos\left(t\right)sin\left(t\right)+{sin\left(t\right)}^4+cos\left(t\right)\left(t^2-1\right)-sin\left(t\right)\left(t^2+1\right)-3t{\mathrm{e}}^t.$

7. Results and Discussion

The first test problem is solved using the proposed method. The problem is linear and is relatively easy to solve. We compare the approximate solution with the exact solution of the problem. We observe that by increasing the scale level of approximation, the approximate solution draws closer to the exact solution as expected; see for example Figure 1a. At $M=8$ the approximate solution (black line) coincides with the exact solution (red dots). In the second part of the same figure, we plot the absolute difference of the exact and approximate solution considering different scales. It is observed that at $M=11$, the absolute error is less than ${10}^{-9}$. This means accuracy up to the ninth decimal place is achieved. We also calculate all the three error norms at different scales. The results are displayed in

Table 1. Error norms of Test Problem 1 at scale level $M=5:13$.

M	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{max}}$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{b}}$
5	$0.1314\times {10}^{-1}$	$0.6431\times {10}^{-1}$	$8.0032\times {10}^{-17}$
6	$1.1355\times {10}^{-2}$	$9.3415\times {10}^{-1}$	$0.4367\times {10}^{-16}$
7	$2.3631\times {10}^{-3}$	$8.3640\times {10}^{-3}$	$3.0669\times {10}^{-16}$
8	$3.9474\times {10}^{-5}$	$3.6641\times {10}^{-4}$	$9.2291\times {10}^{-17}$
9	$1.5963\times {10}^{-5}$	$0.5933\times {10}^{-5}$	$8.3432\times {10}^{-17}$
10	$1.6057\times {10}^{-6}$	$0.5364\times {10}^{-6}$	$7.7554\times {10}^{-17}$
11	$1.5823\times {10}^{-7}$	$1.0592\times {10}^{-7}$	$6.9857\times {10}^{-17}$
12	$7.4855\times {10}^{-9}$	$9.0092\times {10}^{-8}$	$5.4743\times {10}^{-17}$
13	$4.5147\times {10}^{-10}$	$2.4318\times {10}^{-9}$	$2.4362\times {10}^{-17}$

. One can easily note that at $M=5$, the value of is $\left|\right|E_u{\left|\right|}_2=0.1314\times {10}^{-1}$, $\left|\right|E_u{\left|\right|}_{max}=0.6431\times {10}^{-1}$ and $\left|\right|E_u{\left|\right|}_b=8.0032\times {10}^{-17}$. While increasing the scale levels, these values start to decrease with great speed. At $M=13$, these values become $4.5147\times {10}^{-10},$ $2.4318\times {10}^{-9}$ and $2.4362\times {10}^{-17}$, respectively.

The analysis of the second problem is given in

Table 2. Error norms of Test Problem 2 at scale level $M=5:13$.

M	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{max}}$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_{\mathit{max}}$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{b}}$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_{\mathit{b}}$
5	$1.1044\times {10}^{-1}$	$4.3598\times {10}^{-2}$	$3.6499\times {10}^{-1}$	$6.5235\times {10}^{-2}$	$7.9467\times {10}^{-17}$	$6.2402\times {10}^{-16}$
6	$1.9005\times {10}^{-2}$	$9.3636\times {10}^{-3}$	$1.6085\times {10}^{-1}$	$1.6796\times {10}^{-2}$	$9.7023\times {10}^{-16}$	$6.2370\times {10}^{-16}$
7	$2.3621\times {10}^{-3}$	$7.6605\times {10}^{-4}$	$7.5980\times {10}^{-3}$	$1.6380\times {10}^{-3}$	$9.7237\times {10}^{-17}$	$6.2396\times {10}^{-16}$
8	$3.9834\times {10}^{-5}$	$4.9422\times {10}^{-5}$	$1.0331\times {10}^{-4}$	$1.0113\times {10}^{-4}$	$9.7191\times {10}^{-17}$	$6.2400\times {10}^{-16}$
9	$1.4293\times {10}^{-5}$	$5.3436\times {10}^{-6}$	$4.5571\times {10}^{-5}$	$8.6999\times {10}^{-6}$	$9.7195\times {10}^{-17}$	$6.2400\times {10}^{-16}$
10	$1.9227\times {10}^{-6}$	$5.0884\times {10}^{-7}$	$6.1295\times {10}^{-6}$	$7.3992\times {10}^{-7}$	$9.7194\times {10}^{-17}$	$6.2399\times {10}^{-16}$
11	$1.2515\times {10}^{-7}$	$3.1891\times {10}^{-8}$	$4.0112\times {10}^{-7}$	$6.2640\times {10}^{-8}$	$9.7194\times {10}^{-17}$	$6.2399\times {10}^{-16}$
12	$7.1330\times {10}^{-9}$	$1.9583\times {10}^{-9}$	$2.2596\times {10}^{-8}$	$3.6441\times {10}^{-9}$	$9.7194\times {10}^{-17}$	$6.2399\times {10}^{-16}$
13	$4.5177\times {10}^{-10}$	$1.3237\times {10}^{-10}$	$1.4501\times {10}^{-9}$	$2.2449\times {10}^{-10}$	$9.7194\times {10}^{-17}$	$6.2399\times {10}^{-16}$

. We fix the orders of equation to $\sigma =1.8$ and $\gamma =0.8$ and solve the problem using scale level ranges from $M=5$ to $M=13.$ We observe that the error norm decreases rapidly with the increase in scale level. For example, the value $\left|\right|E_u{\left|\right|}_2$ at $M=5$ is $1.1044\times {10}^{-1}$ and at $M=13$, this value drops to $4.5177\times {10}^{-10},$ which is very high accuracy. At the same scale, the value of $\left|\right|E_v{\left|\right|}_2$ is $1.3237\times {10}^{-10}$. The value of norms $\left|\right|E_u{\left|\right|}_{max}$ and $\left|\right|E_v{\left|\right|}_{max}$ are $1.4501\times {10}^{-9}$ and $2.2449\times {10}^{-10}$, respectively. Similarly the approximate solution satisfy the boundary condition with high accuracy. The errors in the boundary condition are $9.7194\times {10}^{-17}$ and $6.2399\times {10}^{-16}$. The error in the boundary condition is observed to be constant, that is, we have the same accuracy at all scale level. The accuracy of the proposed method for all possible values of ${\sigma }_1$ and ${\sigma }_2$ are analyzed by calculating the error norm $\left|\right|E_v{\left|\right|}_2$. The error norm for the solution $u\left(t\right)$ and $v\left(t\right)$ are displayed in Figure 2. We observe that the method produces excellent approximation to the solution at almost all values of parameters.

The same analysis is performed for Test Problem 3 and Test Problem 4. The results are shown in

Table 3. Error norms of Test Problem 3 at scale level $M=5:13$.

N	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{max}}$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_{\mathit{max}}$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{b}}$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_{\mathit{b}}$
5	$6.5505\times {10}^{-3}$	$5.1204\times {10}^{-4}$	$1.5519\times {10}^{-2}$	$8.6002\times {10}^{-4}$	$7.0670\times {10}^{-17}$	$1.8084\times {10}^{-18}$
6	$3.2583\times {10}^{-4}$	$2.8746\times {10}^{-5}$	$7.7139\times {10}^{-4}$	$4.7281\times {10}^{-5}$	$7.0336\times {10}^{-17}$	$1.6297\times {10}^{-18}$
7	$4.6407\times {10}^{-6}$	$4.3808\times {10}^{-7}$	$1.0949\times {10}^{-5}$	$7.1483\times {10}^{-7}$	$7.0320\times {10}^{-17}$	$1.6203\times {10}^{-18}$
8	$4.6870\times {10}^{-7}$	$4.6599\times {10}^{-8}$	$1.1065\times {10}^{-6}$	$7.8519\times {10}^{-8}$	$7.0320\times {10}^{-17}$	$1.6201\times {10}^{-18}$
9	$4.3761\times {10}^{-8}$	$3.4302\times {10}^{-9}$	$1.0365\times {10}^{-7}$	$5.7780\times {10}^{-9}$	$7.0320\times {10}^{-17}$	$1.6201\times {10}^{-18}$
10	$1.3030\times {10}^{-9}$	$1.1448\times {10}^{-10}$	$3.0748\times {10}^{-9}$	$1.8688\times {10}^{-10}$	$7.0320\times {10}^{-17}$	$1.6201\times {10}^{-18}$
11	$1.2252\times {10}^{-11}$	$1.1526\times {10}^{-12}$	$2.8927\times {10}^{-11}$	$1.8785\times {10}^{-12}$	$7.0320\times {10}^{-17}$	$1.6201\times {10}^{-18}$
12	$7.2561\times {10}^{-13}$	$7.2734\times {10}^{-14}$	$1.7118\times {10}^{-12}$	$1.1828\times {10}^{-13}$	$7.0320\times {10}^{-17}$	$1.6201\times {10}^{-18}$
13	$4.6262\times {10}^{-14}$	$3.6748\times {10}^{-15}$	$1.0959\times {10}^{-13}$	$6.1876\times {10}^{-15}$	$7.0320\times {10}^{-17}$	$1.6201\times {10}^{-18}$

and

Table 4. Error norms of Test Problem 4 at scale level $M=5:13$.

N	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_2$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathbf{max}}$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_{\mathbf{max}}$	$\left|\right|{\mathit{E}}_{\mathit{u}}{\left|\right|}_{\mathit{b}}$	$\left|\right|{\mathit{E}}_{\mathit{v}}{\left|\right|}_{\mathit{b}}$
5	$0.0108$	$0.0036249$	$0.030961$	$0.0095518$	$-3.5187\times {10}^{-17}$	$3.7849\times {10}^{-15}$
6	$0.0048582$	$0.0020882$	$0.01299$	$0.0075678$	$3.4932\times {10}^{-17}$	$3.78\times {10}^{-15}$
7	$0.00055459$	$0.00024661$	$0.001525$	$0.00089155$	$3.4845\times {10}^{-17}$	$3.7773\times {10}^{-15}$
8	9.6413 × ${10}^{-16}$	$7.5388\times {10}^{-16}$	$1.8345\times {10}^{-15}$	$1.2652\times {10}^{-15}$	$3.486\times {10}^{-17}$	$3.7777\times {10}^{-15}$
9	$6.8784\times {10}^{-16}$	$1.4702\times {10}^{-15}$	$1.3601\times {10}^{-15}$	$2.8364\times {10}^{-15}$	$3.4861\times {10}^{-17}$	$3.7777\times {10}^{-15}$
10	$4.5115\times {10}^{-15}$	$5.6863\times {10}^{-16}$	$1.302\times {10}^{-14}$	$8.6167\times {10}^{-16}$	$3.4861\times {10}^{-17}$	$3.7777\times {10}^{-15}$
11	$8.4757\times {10}^{-17}$	$2.7414\times {10}^{-16}$	$1.6027\times {10}^{-16}$	$4.7286\times {10}^{-16}$	$3.4861\times {10}^{-17}$	$3.7777\times {10}^{-15}$
12	$1.5906\times {10}^{-14}$	$3.1092\times {10}^{-15}$	$4.6821\times {10}^{-14}$	$6.3783\times {10}^{-15}$	$3.4861\times {10}^{-17}$	$3.7777\times {10}^{-15}$
13	$4.926\times {10}^{-15}$	$1.0866\times 10-15$	$1.4378\times {10}^{-14}$	$2.0953\times {10}^{-15}$	$3.4861\times {10}^{-17}$	$3.7777\times {10}^{-15}$

. The same conclusion is reached for these two examples. The error bounds are shown in Figure 3 and Figure 4. We observe that the method yields an almost accurate solution for all values of these parameters. In Figure 2, we can easily see that the error norm is less than ${10}^{-3}$. Note that for these problems we have set $M=5$. It is always possible to get a more accurate solution by selecting the highest choices of scale level.

The nonlinear Test Problem 5 is solved with the proposed iterative scheme. We use three different choices of the parameters $\sigma$ and $\gamma$ and calculate the the error norms $\left|\right|E_u{\left|\right|}_2$ at different iterations. We use seven iterations for the approximation of solution. The results are displayed in Figure 5 and Figure 6. We observe that the error decreases with respect to the number of iterations and is highly convergent. Note that for this example, we fix the scale level $M=10$. It is clear form the figure that the proposed method is highly efficient, especially in the solution of nonlinear equations. The same phenomenon is observed for Test Problem 6. The results are displayed in Figure 7.

8. Conclusions and Future Work

This article presents a new algorithm for the solution of fractional order differential equations. The Newton Raphson method is combined with the operational matrix method for the solution of these problems. The convergence of the proposed method is checked analytically and is confirmed by solving several test problems. It is found that the approximate solution is highly accurate, and one can get high accuracy by using high scale levels. The mathematical proof of convergence and error analysis is our future plan of research.