Approximate Analytical Solutions for Systems of Fractional Nonlinear Integro-Differential Equations Using the Polynomial Least Squares Method

: We employ the Polynomial Least Squares Method as a relatively new and very straightforward and efﬁcient method to ﬁnd accurate approximate analytical solutions for a class of systems of fractional nonlinear integro-differential equations. A comparison with previous results by means of an extensive list of test-problems illustrate the simplicity and the accuracy of the method.


Introduction
The notion of a fractional derivative as a derivative of any arbitrary real or complex order has a long history, starting with the works of early titans of mathematics such as Leibniz, Abel, Liouville and Heaviside. The theory and applications of the fractional calculus expanded steadily during the 20th century and exploded in the recent decades due to the numerous applications of the fractional differential equations in various field of science and engineering, such as thermal engineering, acoustics, electromagnetism, control, robotics, viscoelasticity, signal processing etc.
The fractional integro-differential equations are in the present the focus of intensive research due to their pivotal role in the modeling of many phenomena and processes from science and engineering such as, for example, electric-circuit analysis, floating structures and viscoelastic material dynamics.
The class of systems of equations studied in this paper is: b a K f k (t, s, x 1 (s), . . ., x n (s))ds, t a K vk (t, s, x 1 (s), . . ., x n (s))ds , together with a set of conditions of the type: Here, for q ∈ N * , D α denotes the Caputo fractional derivative of order α, namely: The kernel functions K f k , K vk and the functions F k are assumed to be suitably defined on the [a, b] interval such that the problem consisting of the Equation (1) together with the conditions (2) admits a solution.
This class of systems of equations is evidently a very general one since it includes both Fredholm and Volterra type equations, linear and nonlinear. Unfortunately, with the exception of a relatively small number of simple cases (such as most of the test problems included as examples), the exact solutions of a nonlinear integro-differential system of equations of the type (1) cannot be found. Thus numerical solutions or (preferably) approximate analytical solutions must be computed.
In recent years many methods and techniques to find approximate solutions for systems of integro-differential equations were proposed, among which we mention: • In 2006, Momani and Qaralleh used the Adomian Decomposition Method [1] to find approximate solutions for a type of system of equations similar to (1) but including only the Volterra term (and not the Fredholm one). The method decomposes the solution of the problem into a rapidly convergent series and replaces the nonlinear term by a series of the Adomian polynomials. The method worked well but unfortunately no error reports were included in the paper. • In 2009, Zurigat et al. employed the well-known Homotopy Analysis Method [2] and in 2020 Akbar et al. employed the Optimal Homotopy Analysis Method [3] to find approximate solutions for systems of fractional integro-differential equations also of Volterra type only. These methods are based on the concept of the homotopy from topology and generate a convergent series solution for nonlinear systems. The experimental results were good with low errors. An interesting feature of the Homotopy Analysis Method is the existence of the auxiliary parameters which can be adjusted to control the convergence region of the solution series. Unfortunately the best choices of these parameters are not always clear. On the other hand, in [2] the authors showed that for certain values of these parameters the Homotopy Analysis Method can also replicate the results obtained in [1] by means of the Adomian Decomposition Method. • In 2010, Saeed and Sdeq employed the widely used Homotopy Perturbation Method [4] to find approximate solutions for systems of linear Fredholm fractional integrodifferential equations. The Homotopy Perturbation Method combines the traditional perturbation method and the concept of homotopy from topology and its best feature is that it does not require the existence of a small parameter regarding the perturbation.
The errors corresponding to the approximations are presented and the convergence of the method is fast: from the experimental data it follows that the order of the error seems to be proportional to the number of terms considered in the sum of the series. • The Chebyshev Pseudo-Spectral Method, employed in 2013 by Khader and Sweilam [5] and in 2017 by Zedan et al. [6] to find approximate for a type of system of equations similar to (1), but including only Volterra integrals, uses the properties of Chebyshev polynomials to reduce the problem to a linear or non-linear system of algebraic equations. From the numerical data the convergence again appears fast, with an order of the error roughly proportional to the numbers of terms in the series sum. • In 2014, Bushnaq et al. presented the Reproducing Kernel Hilbert Space Method [7], which is a kernel-based approximation method, to find approximate solutions for systems of Volterra fractional integro-differential equations. The errors presented are relatively low but, while the absolute convergence of the method is proved, from the experimental data no information about the speed of the convergence can be extracted. • Chebyshev Wavelets Expansion Methods were used in 2014 by Heydary et al. [8] to solve systems of nonlinear singular fractional Volterra integro-differential equations, in 2018 by Zhou and Xu [9] and in 2021 by Bargamadi et al. [10] to solve fractional Volterra-Fredholm integro-differential equations. This category of methods utilizes Chebyshev wavelets as a basis and transforms the problem into a system of algebraic equations. These methods usually yield solutions with very low errors which converge relatively fast. The rest of the paper is structured as follows: In Section 2, we present the Polynomial Least Squares Method (denoted from this point forward as PLSM), in Section 3, we present the results of an extensive testing process involving most of the usual test problems included in similar studies and in Section 4 we present the conclusions of the study.

The Polynomial Least Squares Method
In the following we will denote by the problem (1)+(2) the system of Equation (1) together with the conditions (2). Letx 1 (t), . . .,x n (t) be a set of approximate solutions of the system (1). The overall error obtained by replacing the set of exact solutions x 1 (t), . . ., x n (t) with these approximations can be described by the remainder: We will find a set of approximate polynomial solutions of (1)+(2) on the [a, b] interval, solutions which satisfy the following conditions: We call a set of ε-approximate polynomial solutions of the problem (1)+(2) a set of approximate polynomial solutionsx 1 , . . .,x n satisfying the relations (5) and (6).

Definition 2.
We call a set of weak δ-approximate polynomial solutions of the problem (1)+(2) a set of approximate polynomial solutionsx 1 , . . .,x n satisfying the relation b a R(t,x 1 , . . .,x n )dt ≤ δ together with the initial conditions (6).
The following convergence theorem holds: Theorem 1. The necessary condition for the problem (1)+(2) to admit a set of sequences of polynomials P km (t) convergent to the solutions of this problem is: lim where T km (t) are a set of weak ε-approximate polynomial solutions of the problem (1)+(2).
Proof. We will find a set of weak ε-polynomial solutions of the type: where the constants c k0 , c k1 , . . ., c km , k = 1, . . ., n are calculated using the steps outlined in the following.

Remark 1.
Any set of ε-approximate polynomial solutions of the problem (1)+(2) is also a set of weak ε 2 · (b − a)-approximate polynomial solutions, but the opposite is not always true. It follows that the set of weak approximate solutions of the problem (1)+(2) also contains the set of approximate solutions of the problem.
Taking into account the above remark, in order to find a set of ε-approximate polynomial solutions of the problem (1)+(2) by the Polynomial Least Squares Method, we will first determine weak approximate polynomial solutions, T 1m (t), . . ., T nm (t).

Numerical Examples
In this section we present the results obtained by using PLSM for many (if not most) of the usual test-problems included in similar papers.
In the first application, we present the computations in detail including some remarks about the practical implementation of the method, while in the most of the subsequent applications we only present the results of the computations.

Application 1: System of Fractional Fredholm Integro-Differential Equations
The first application consists of the pair of equations [4,14,20]: The exact solutions of this problem are x 1e (t) = t, x 2e (t) = t 2 .
We will follow the steps of the algorithm described in the Proof from the previous section. First, by choosing the polynomialsx 1 (t) = c 11 · t + c 10 andx 2 (t) = c 22 · t 2 + c 21 · t + c 20 , from the initial conditions it follows that c 10 = 0 and c 20 = 0, hencex 1 The corresponding remainder (8) is: and the corresponding functional (9) is: In order to find the minimum of the functional (9), we solve the system consisting of the equations dJ dc 11 = 0, dJ dc 21 = 0 and dJ dc 22 = 0 and we obtain the solution c 11 = 1, c 21 = 0, c 22 = 1 which is a set of stationary (critical, equilibrium) points for the functional. It is easy to show that this set corresponds indeed to a minimum of J and we obtain the values From the initial conditions, obviously c 0 10 = 0 and c 0 20 = 0 and thusx 1 (t) = t andx 2 (t) = t 2 which means that PLSM is able to find, in a very simple manner, the exact solution of the problem. We remark that the previous methods in [4] (Homotopy Perturbation Method) and [14] (Iterative Refinement Method) were only able to find approximate solutions.
Regarding the practical implementation of the method, we wish to make the following remarks: • Regarding the choice of the degree of the polynomial approximation, in the computations we usually start with the lowest degree (i.e., first degree polynomial) and compute successively higher degree approximation, until the error (see next item) is considered low enough from a practical point of view for the given problem (or, in the case of a test problem, until the error is lower than the error corresponding to the solutions obtained by other methods). Of course, in the case of a test problem when the known solution is a polynomial, one may start directly with the corresponding degree, but this is just a shortcut and by no means necessary when using the method.
• If the exact solution of the problem is not known, as would be the case of a real-life problem, and as a consequence the error can not be computed, then instead of the actual error we can consider as an estimation of the error the value of the remainder R (4) corresponding to the computed approximation, as mentioned in the end of Section 2.

•
If the problem has an (unknown) exact polynomial solution it is easy to see if PLSM finds it since the value of the minimum of the functional in this case is actually zero.
In this situation, if we keep increasing the degree (even though there is no point in that), from the computation we obtain that the coefficients of the higher degrees are actually zero. • Regarding the choice of the optimization method used for the computation of the minimum of the functional (9), if the solution of the problem is a known polynomial (such as in the case of this application and several of the following ones) we usually employ the critical (stationary) points method, because in this way by using PLSM we can easily find the exact solution. Such problems are relatively simple ones; the expression of the functional (9) is also not very complicated and indeed the solutions can usually be computed even by hand (as in the case of this application) and in general no concerns of conditioning or stability arise. However, for a more complicated (real-life) problem, when the solution is not known (or even if the exact solution is known but not polynomial), we would not use the critical points method. In fact, we would not even use an iterative-type method but rather a heuristic algorithm such as Differential Evolution or Simulated Annealing. In our experience with this type of problem even a simple Nelder-Mead type algorithm works well (as was the case for the following Application 7, Application 8 and Application 9).
The corresponding functional (9) is: Again, by computing the stationary (critical, equilibrium) points for the functional, it is easy to show that the minimum corresponds to c 0 11 = 1, c 0 21 = −1, thus finding again the exact solution of the problem.
In the following applications, when the exact solutions are also polynomial, we will omit the details of the computations, since they are very similar.
We also remark again that the previous method in [17] (Block-Pulse Functions Method) was only able to find approximate solutions.

Application 3: System of Fractional Volterra-Fredholm Integro-Differential Equations
The third application consists of the equations [17]: together with the initial conditions x 1 (0) = 0, x 2 (0) = 0. The exact solutions of this problem are x 1e (t) = t 2 , x 2e (t) = −t 2 . Again, using PLSM we are able to find the exact solutions of the problem, while the previous method in [17] (the Block-Pulse Functions Method) was only able to find approximate solutions.

Application 4: System of Fractional Volterra-Fredholm Integro-Differential Equations with a Weakly Singular Kernel
The next application consists of the equations [21]: together with the boundary conditions The exact solutions of this problem are x 1e (t) = t 2 + t, x 2e (t) = t 3 − 1.
Again, using PLSM we are able to find the exact solutions of the problem, while the previous method in [21] (Müntz-Legendre Wavelet Method) was only able to find approximate solutions.
The exact solutions of this problem are x 1e (t) = t, x 2e (t) = −t. Using PLSM we can easily find the exact solutions of the problem, while the previous method in the [10] (Second Chebyshev Wavelets Method) was only able to find approximate solutions.

Application 6: System of Fractional Volterra-Fredholm Integro-Differential Equations
The next application consists of the equations [9]: together with the boundary conditions: The exact solutions of this problem are x 1e (t) = t 2 , x 2e (t) = t 3 .
Once again, using PLSM we can easily find the exact solutions of the problem, while the previous method in [9] (Chebyshev Wavelets Method) was only able to find approximate solutions.
The exact solutions of this problem are only known for α 1 = α 2 = 1, in which case they are x 1e (t) = t + e t , x 2e (t) = t − e t .
This problem is one of the most frequently used test-problems for this type of equation: in 2006 in [1] approximations were computed by means of the Adomian Decomposition Method (but unfortunately no estimations of the error were presented); in 2009 in [2], approximations were computed using the Homotopy Analysis Method (but again no values for the errors were given); in 2013 in [5], approximations were computed by employing the Chebyshev Pseudo-Spectral Method; in 2014 in [7], approximations were computed by using a Reproducing Kernel Hilbert Space Method; in 2015 in [13], approximations were computed by using the Triangular Functions Operational Matrix Method; in 2017 in [6], approximations were computed by employing a Chebyshev Spectral Method; in 2018 in [16], approximations were computed using a Hybrid Bernstein Block-Pulse Functions Method; in 2021 in [22], approximations were computed by employing a Finite Integration Method based on Shifted Chebyshev Polynomials.
The problem was solved in 2020 in [21] where approximations were computed by means of the Müntz-Legendre Wavelets Method.
The exact solutions of this problem are only known for α 1 = α 2 = 1, in which case they are x 1e (t) = sinh(t), x 2e (t) = cosh(t). For this case, in Table 3 we will compare the absolute errors corresponding to the approximations computed in [21] with the absolute errors corresponding to several approximations computed by using PLSM. Since unfortunately in [21] the results were not presented for a given set of values of t, but as plots of the absolute errors, in Table 5 we also present the corresponding maximal absolute errors. Again, the results illustrate the accuracy and the fast convergence of PLSM: Table 5. Comparison of the maximal absolute errors of the approximate solutions x 1 (first row) and x 2 (second row) for problem (18) in the case α 1 = α 2 = 1. We also computed several pairs of approximate solutions corresponding to fractional values for α 1 and α 2 , whose plots are presented in Figures 3 and 4 together with the plots corresponding to the exact solutions in the case α 1 = α 2 = 1.

Application 9: System of Fractional Volterra-Fredholm Integro-Differential Equations
The next application consists of the equations [17]: together with the initial conditions x 1 (0) = 0, x 2 (0) = 0. The problem was studied in 2019 in [17] where approximations were computed by means of a Block-Pulse Functions Method, but unfortunately no errors corresponding to the approximations were presented.
The exact solutions of this problem are only known for α 1 = α 2 = 1, in which case they are x 1e (t) = sin(π · t), x 2e (t) = sinh(t). For this case, in Table 6 we will present the maximal absolute errors corresponding to several approximations computed by using PLSM. Table 6. Comparison of the maximal absolute errors of the approximate solutions x 1 (first row) and x 2 (second row) for problem (19) in the case α 1 = α 2 = 1.
We also computed several pairs of approximate solutions corresponding to fractional values for α 1 and α 2 , whose plots are presented in Figures 5 and 6 together with the plots corresponding to the exact solutions in the case α 1 = α 2 = 1.

Application 10: System of Singular Fractional Volterra Integro-Differential Equations
The last application consists of the equations [8,16]: together with the initial conditions x 1 (0) = 0, x 2 (0) = 0, x 3 (0) = 0. The problem was studied in 2014 in [8], where approximations were computed by means of a Chebyshev Wavelets Expansion Method and in 2018 in [16] by means of a Block-Pulse Function Method.
The exact solutions of this problem are only known for α 1 = α 2 = α 3 = 1, in which case they are x 1e (t) = t 2 − t, x 2e (t) = t 2 , x 3e = t 3 . For this case, using PLSM we are able to compute the exact solution.
We also computed several pairs of approximate solutions corresponding to fractional values for α 1 and α 2 , whose plots are presented in Figures 7-9.

Conclusions
The study presents the Polynomial Least Squares Method as a straightforward, efficient and accurate method to compute approximate analytical solutions for a very general class of systems of fractional nonlinear integro-differential equations.
The main advantages of PLSM are: • The simplicity of the method-the computations involved in the use of PLSM are as straightforward as it gets; in fact in the case of a lower degree polynomial the computations sometimes can be performed even by hand, as illustrated by the first application. • The accuracy of the method-clearly illustrated by the applications presented, since by using PLSM we were able to compute approximations at least as good (if not better) than the approximations computed by other methods. • The simplicity of the approximation-since the approximations are polynomial, they also have the simplest possible form and thus any subsequent computation involving the solution can be performed with ease. While it is true that for some approximation methods which work with polynomial approximations the convergence may be very slow, this is not the case here (see for example Application 7, Application 8 and Application 9, which are representative of the performance of the method).
The study includes an extensive application list, containing most of the usual test problems used for this type of equation and compare our results with previous results obtained by using other well-known methods. For the test problems whose exact solution is a polynomial one, PLSM is able to find the exact solution in a simple manner, while most of the other methods previously used were only able to compute approximate solutions. If the solution is not polynomial, PLSM is able to find approximate solutions, again in a very straightforward way, with errors usually smaller that the ones corresponding to the approximations computed by other methods.
Taking into the account the above considerations, the results of this study recommend PLSM as a very useful tool in the research of systems of fractional nonlinear integrodifferential equations.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

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