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

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.


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].
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 where 0 < γ 1 , The scalar functions u(t) and v(t) are the solution to be determined. f and g are nonlinear functions of u(t), v(t) 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.

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].
, R) is defined by the following relation.
From the above definition, it is easy to extract The Shifted Legendre Polynomials (LP) These polynomials play a central role in approximation theory. Generally, these polynomials are defined on the domain [0, τ], which is given by where These polynomials enjoys a very important property of orthogonality on the domain [0, τ], which is expressed mathematically as By using Equation (5), any s(t) can be expressed in terms of LP as The above equation has an equivalent vector representation given as where and 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].
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 ∏ M be the space of M terms Legendre polynomials and let u(t) ∈ C m [0, 1], then u m (t) is in space ∏ M ; then, for the m term approximation, C is constant and can be chosen in such a way that u (2m) belongs to ∏ M , where u (m) is defined as where Z is storm livoliel operator of legendre polynomials with u (0) = u(t).
In the next section, we first recall some of our previously designed operational matrices and then develop new operational matrices.

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 Λ τ M (t) be the function vector; the α order integration is generalized as where H τ,α M×M is the OP of integration, defined as .

Proof.
Consider Then, using the previous result, we get Θ a,j,τ ρ τ j (t).

Lemma 3 ([24]
). Let Λ τ M (t) be the function vector as defined in (8); then the derivative of order σ M×M is the operational matrix of derivative of order σ, and G τ,σ M×M is defined as Furthermore, ‫ג‬ (i,k) is similar as defined in (4) and .
with operational matrix of derivative is bounded by the following.
The proof of this corollary is similar as Corollary 1.

Lemma 4 ([24]
). Let u(t) and φ n (t) be smooth functions that are well-defined on [0, τ]. Then where W M is the Legendre coefficients vector of u(t) as defined in (7) and The matrix G τ,σ M×M is the operational matrix of derivative as defined in Lemma 3; the entries of matrix J τ,φ n M×M are defined by the following relation .
with operational matrix of fractional derivative with variable coefficient is bounded by the following.
Then, using the relation above, we get Truncating the sum and writing in modified form we get We can also write it in matrix form as 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(t) be a known function well defined on [0, τ] and φ n c (t) = ct n be a polynomial then, for a function vector Λ τ M (t) as defined in (8), the following result holds where the matrix Q c,n,s(t) M×M is given as . where , d i is the Legendre coefficients of the function s(t) and ‫ג‬ (j,r) is as defined in (4).
Proof. Let s(t) be approximated with Legendre polynomials, as The polynomial τcd i t n 2i+1 can be expressed as a series of Legendre polynomials as where the constant Ω (i,j) is given by Using the definition of ρ τ j (t) and after simplification, we can write Simulating the result for i = 0 : M and j = 0 : M completes the proof of the Lemma.

Lemma 6.
Let φ n c (n) = ct n be a polynomial then for a function vector Λ τ M (t) as defined in (8); The matrix R c,n,τ M×M is defined as Furthermore, the entries are defined as Proof. The general term of the equality can be written as It can be approximated with Legendre polynomials Using the definition of Legendre polynomials we can write which completes the proof of the lemma.

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.

Linear FDEs with Variable Coefficients
We first consider the following coupled system of linear fractional differential equations with variable coefficients where m 1 and m 2 are some real constants. 0 < γ 1 , Assume the solution of (10) in terms of shifted Legendre polynomials, such that the following hold Applying fractional integral of order σ 1 on the first equation of (11) and making use of Lemma 2, we can write 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 . Now using the values of e 0 and e 1 in (12), we can write u(t) as In view of Lemma 5 and Lemma 6, we can write Equation (13) as where u 0 − s 0 t s 1 = F 1 Λ τ M (t). For simplicity of notation, we can write the above equations as where

M×M
Repeating the same procedure for v(t), we can get analogous representation as where Now, in view of (15), (14), (11), and Lemma 4, the equivalent matrix form for system (10) is given as Canceling out the common term and after simplification of notation, we can write where and The above equations can also be written as We see that (18) is a linear system of matrix equations. By solving (17), we will get the required coefficients vector K M and L M , which can be used in (14) and (15) to get approximation to the solution of the main problem.

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 (t) and v 0 (t). • Linearize the nonlinear part at u 0 (t) and v 0 (t). 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 (t) and v 1 (t). • Repeat step 1.
Consider the following nonlinear FDEs.
separating the linear and nonlinear parts, we get First solve the linear part: Its solution is named u 0 (t) and v 0 (t). The next step is to linearize the nonlinear part.
We get a system of FDEs with variable coefficients. The whole process can be expressed as The boundary conditions can be written as It can be easily noted that (23) is fractional differential equation with variable coefficients.

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.

Error Bound for Single Differential Equation
Consider the following fractional differential equation.
subject to the following initial and boundary conditions Our aim is to derive an upper bound for the proposed method. We have to calculate R M defined as The solution of the above problem can be written in terms of shifted Legendre series such that Applying a fractional integral of order σ, using an operational matrix of integration and using Corollary (1), we can write Which can be simplified as We know from (14) and the integral type boundary conditions that we can conclude, Assume thatX = K M E M×M + F 1 . Using Lemma (3) and Corollary (2), we can write Approximating and using (30) in (25) we get Using the bounded property of Legendre polynomail, it follows that In view of Theorem (1), it is evident that if the function u(t) and f (t) are sufficiently smooth functions, then the sequence that defines their coefficient is convergent to zero. Hence, we conclude that as m → ∞ the coefficients u m → 0 and f m → 0. Hence it can be easily observed that the error |R M (t)| → 0. Equation (31) also establishes an upper bound of the error between the exact and approximate solution.

Error Bound for Coupled System of Fractional Differential Equations
Consider the following system of FDEs.
subject to the following initial and boundary conditions Our aim is to derive an upper bound for the proposed method. We have to calculate R u M and R v M , defined as We know from (14) that Corollary (2), we can write Approximating d k ρ τ k (t) and using (30) in (25) we get where R u M (t) and R v M (t) is defined by the relation Using the bounded property of the Legendre polynomial, it follows that The above equation establishes an error bound for the solution u(t) and v(t). It also ensures the convergence of the proposed method for the solution of coupled system of FDEs.

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.
We check the accuracy of the boundary condition by measuring the following error norms.
In the above bounds, the quantity U M (t) represents the m− term approximation to the solution u(t).
where the exact solution u(t) = t 3 + t 2 + t + 4, and the source term
where the exact solution u(t) = sin(t) and v(t) = e t , and the source term

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. One can easily note that at M = 5, the value of is ||E u || 2 = 0.1314 × 10 −1 , ||E u || max = 0.6431 × 10 −1 and ||E u || b = 8.0032 × 10 −17 . While increasing the scale levels, these values start to decrease with great speed. At M = 13, these values become 4.5147 × 10 −10 , 2.4318 × 10 −9 and 2.4362 × 10 −17 , respectively.  The analysis of the second problem is given in Table 2. We fix the orders of equation to σ = 1.8 and γ = 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 ||E u || 2 at M = 5 is 1.1044 × 10 −1 and at M = 13, this value drops to 4.5177 × 10 −10 , which is very high accuracy. At the same scale, the value of ||E v || 2 is 1.3237 × 10 −10 . The value of norms ||E u || max and ||E v || max are 1.4501 × 10 −9 and 2.2449 × 10 −10 , respectively. Similarly the approximate solution satisfy the boundary condition with high accuracy. The errors in the boundary condition are 9.7194 × 10 −17 and 6.2399 × 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 σ 1 and σ 2 are analyzed by calculating the error norm ||E v || 2 . The error norm for the solution u(t) and v(t) 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 Tables 3 and 4. The same conclusion is reached for these two examples. The error bounds are shown in Figures 3 and 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 σ and γ and calculate the the error norms ||E u || 2 at different iterations. We use seven iterations for the approximation of solution. The results are displayed in Figures 5 and 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.

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.

Data Availability Statement:
The data used to support the findings of this study are included within the article.

Conflicts of Interest:
The authors declare no conflict of interest.