An Enhanced Adaptive Bernstein Collocation Method for Solving Systems of ODEs

In this paper, we introduce two new methods to solve systems of ordinary differential equations. The first method is constituted of the generalized Bernstein functions, which are obtained by Bernstein polynomials, and operational matrix of differentiation with collocation method. The second method depends on tau method, the generalized Bernstein functions and operational matrix of differentiation. These methods produce a series which is obtained by non-polynomial functions set. We give the standard Bernstein polynomials to explain the generalizations for both methods. By applying the residual correction procedure to the methods, one can estimate the absolute errors for both methods and may obtain more accurate results. We apply the methods to some test examples including linear system, non-homogeneous linear system, nonlinear stiff systems, non-homogeneous nonlinear system and chaotic Genesio system. The numerical shows that the methods are efficient and work well. Increasing m yields a decrease on the errors for all methods. One can estimate the errors by using the residual correction procedure.


Introduction
Many real life phenomena can be modeled by systems of ordinary differential equations (ODEs). For instance, the mathematical models of circuits and mechanical systems involving several springs connected in series can be given by a system of differential equations. Generally, such systems are frequently encountered in chemical, ecological, biological and engineering applications [1]. Various phenomena in chemical kinetics and engineering are modeled with the stiff systems [2]. Explicit numerical methods may solve these problems with some limitations on the step size which yields computational complexity [3]. In control theory, ODE systems also have chaotic behaviors [4,5]. A chaotic system is a structure that exhibits a sensitive dependence on initial conditions and is a nonlinear deterministic system with complex and unpredictable behavior. The Genesio system is an example of such a system [6]. It is one of the chaos paradigms as it has many properties of chaotic systems.
A system of ODEs can be expressed in the form u j = f j (x, u 1 , . . . , u r ), u j (x 0 ) = u 0,j , j = 1, 2, . . . , r, where f j are real-valued functions, x 0 and u 0,j are real numbers. By applying the variable transformation x → x + x 0 , the systems (1) can be defined around the origin, and so we consider solving the equations of the form u j = f j (x, u 1 , . . . , u r ), u j (0) = α j , j = 1, 2, . . . , r.
Different numerical integration algorithms such as Runge-Kutta method for approximating solutions of the systems (2) have been proposed in the literature. However, these algorithms calculate the values of approximate solutions on the nodes instead of giving a solution over the interval. Approximate analytical solutions of certain classes of systems of ODEs based on the homotopy analysis method and homotopy perturbation method have been given in [3,7] respectively. A relatively new analytical method based on the Bernstein polynomials has been shown to be a promising method for solving linear and non-linear equations. Isik et al. [8] presented an approximate method based on the Bernstein polynomials for solving high order linear differential equations. Approximate Bernstein series solutions of fractional heat-and wave-like equations were given by Rostamy and Karimi [9]. Yuzbasi [10] and Baleanu et al. [11] presented approximate analytical methods constituted of the Bernstein polynomials for solving fractional Riccati type differential equations. Bernstein series solutions of Lane-Emden type equations were given by Pandey and Kumar [12] and Isik and Sezer [13]. Bernstein series solutions with a priori error estimate for linear second-order partial differential equations with general conditions were given by Isik et al. [14]. Maleknejad et al. [15] proposed a numerical method for solving the systems of high order linear Volterra-Fredholm integro-differential equations by using Bernstein operational matrices. Rostamy and Karimi [16] presented a numerical method consists of the high-order derivative matrix of the Bernstein polynomials. Multistage Bernstein polynomials (MB-polynomials) method which is a modification of Bernstein polynomials method was developed by Alshbool and Hashim [17] for solving fractional-order stiff systems. Bernstein operational matrix of derivative was adapted to solve linear and non-linear fractional differential equations by Alshbool et al. [18]. Asgari and Ezzati [19] solved two-dimensional fractional integral equations by two-dimensional Bernstein polynomials operational matrix. An approximate solution method, called multistage Bernstein collocation method, to solve strongly nonlinear damped systems was given in [20]. Khataybeh et al. [21] demonstrated for the first time the applicability of the operational matrices of Bernstein polynomials method for solving directly third-order ODEs. Direct solution of second-order system of ODEs using Bernstein polynomials was presented in [22]. Bataineh et al. [23] presented a two-dimensional Bernstein polynomials method for solving time-dependent Emden-Fowler type of equations. The Bernstein polynomials method incorporating residual correcting procedure were applied to a system of secondorder BVPs, Brusselator system and nonlinear stiff system by Alshbool et al. [24]. Very recently, Alshbool et al. [25] solved a class of fractional diffusion equations by fractional Bersnsetin series solution.
In this study, we present two new methods, namely generalized Bernstein function (GBF) tau and GBF collocation methods, to numerically solve the systems of ODEs. The methods are obtained by a special generalization of m-th degree Bernstein polynomials and collocation or tau methods. To introduce the methods, we first give the definitions of Bernstein polynomials (BPs) and GBFs in Section 2. We give the approximation functions in the same section. In Section 3, we give the operational matrices for BPs and GBFs. Then, the new methods are formed by using tau method or collocation method. We also give Bernstein collocation method and Bernstein tau method to show how these new methods are generated. Residual correction procedure is modified for all methods. The numerical stabilities of the methods are also given in Section 3. Several examples are studied to demonstrate the accuracy and efficiency of the methods. We apply the methods for different values of m to show the dependency of m values.

Existence and Uniqueness Theorem
Let the function f (x, u 1 , . . . , u r ) be defined on the set Then we say f satisfies the Lipschitz condition on D in the variables u i for i = 1, 2, . . . , r if there exists a constant L > 0 such that As a result of the mean value theorem, f satisfies the Lipschitz condition on D in the variables u i for i = 1, 2, . . . , r, if f and its first partial derivatives are continuous on D and if for all (x, u 1 , . . . , u r ) ∈ D. The existence and uniqueness theorem can be found in [26,27].
Theorem 1. Suppose f i , i = 1, 2, ..., r be continuous and satisfy a Lipschitz condition on the set Then, the system of (2) subject to the initial conditions has a unique solution for 0 ≤ x ≤ T 0 .

Bernstein Polynomials
The BP of degree m are defined by where the binomial coefficient is There are m + 1 BPs which are degree m that are formed as a base for the n dimensional polynomial space. It is set B i,m = 0 in case of i < 0 or i > m.

Generalized Bernstein Functions
Let us define a GBF of degree m s similar to (3) as Again we get m + 1 functions for m s th-degree GBF and set B i,m = 0 for i < 0 or i > m. We should note thatB i,m for 0 < i < m have to be continuous on [0, ∞).

Approximation of Functions
First at all, we approximate the functions u j (x) and u j (x) for j = 1, 2, . . . , r with the mth-degree BP and the m s th-degree GBF respectively as follows where c 0,j , c 1,j , . . . , c m,j are to be determined. Let us call u j,m and u j,m as BP series solution obtained by collocation method (BPSSC) or tau method (BPSST) and GBF series solution obtained by collocation method (GBFSSC) or tau method (GBFSST).

Applications of Operational Matrices
In this section, we will obtain the approximate solutions of systems (2). We will first consider tau methods, i.e., we will obtain BPSST and GBFSST by using the operational matrices of mth-degree BP and the m s th-degree GBF, respectively. To solve (2) by means of the operational matrices, we employ Equations (5)-(8) and then we define the residuals (x) and (x) for Equation (2) respectively as

The Approximate Solutions Obtained by Tau Method
Multiplying and integrating (x) and (x) yields m equations sets Also, by imposing the initial conditions of Equation (2) into Equations (5) and (7) we have Equation (11) with Equation (13) or Equation (12) with (14) generate m + 1 sets of equations respectively. Solving these equations gives the unknown coefficients c 0,j , c 1,j , . . . , c m,j . Consequently, u j,m (x) and u j,m (x) given in Equations (5) and (7) can be calculated, respec-tively. Thus, BPSST and GBFSST are obtained.
Let us call the standard tau method with Bernstein polynomials which yields BPSST solution as Bernstein tau method. Similarly, let us call the new method depending on tau method and GBF which produces GBFSST as GBF tau method..

Residual Correction Procedure for Bernstein Tau Method and GBF Tau Method
We will constitute the residual correction procedure for Bernstein tau method. Let u 1,m , u 2,m , . . . , u r,m be the approximate solution set of the system (2). Adding the equality u j,m = u j,m into the both sides of (2) yields as where e j (x) := u j (x) − u j,m (x) and u j (x) is the exact solution. A similar argument for the initial conditions yields e j (0) = 0, j = 1, 2, . . . , r.
Let us approximate to e j by using the present method Let us define the residue e (x) as Then, the constants c e 0,j , c e 1,j , . . . , c e m,j can be obtained by constructing m + 1 sets of equations by applying a typical tau method with the initial conditions (16). Thus, e j,m (x) can be obtained by solving these sets of linear or nonlinear equations.
The following results are the same when the residual correction procedure is used for the error estimates. Let · be any norm defined on continuous function space. If e j − e j,m < , j = 1, 2, . . . , r where > 0 is sufficiently small, then the absolute errors e j can be estimated by e j,m for j = 1, 2, . . . , r, respectively. Hence, the optimal m for the absolute errors may be obtained measuring the error functions e j,m for different m values in any norm. If u j,m , j = 1, 2, . . . , r are the BPSST of (2), then u j,m + e j,m , j = 1, 2, . . . , r are also approximate solutions for (2). Moreover, their error functions are e j − e j,m , j = 1, 2, . . . , r.
Note that the approximate solution set u j,m + e j,m : j = 1, 2, . . . , r is a better approximation set than u j,m : j = 1, 2, . . . , r in the norm if e j − e j,m ≤ u j − u j,m where u j : j = 1, 2, . . . , r are the exact solution of (2). Let us call the approximate solutions u j,m + e j,m , j = 1, 2, . . . , r as corrected BPSST.
Similar arguments can be done to estimate the error obtained by GBF tau method. For j = 1, 2, . . . , r, adding the terms u j,m into the both side of the j-th equation in (2) gives e j (x) − f j ( e 1 + u 1,m , . . . , e r + u r,m ) = − u j,m (x) (19) where e j (x) := u j (x) − u j,m (x) and u j (x) is the exact solution with the conditions e j (0) = 0, j = 1, 2, . . . , r.
Approximating e j by using the method Finally, we get the c e 0,j , c e 1,j , . . . , c e m,j by solving tau method to the following system Hence, e j,m (x) can be obtained by solving these sets of equations. In case of e j − e j,m < , j = 1, 2, . . . , r where > 0 is sufficiently small, the absolute errors e j can be estimated by e j,m for j = 1, 2, . . . , r, respectively. Again, the optimal m for the absolute errors might be obtained by measuring the errors e j,m . We obtain another solutions, namely corrected GBFSSTs, by adding the error to the GBFSSTs u j,m + e j,m , j = 1, 2, . . . , r.

Approximate Solutions Obtained by Collocation Method
Let the collocation nodes be {0 ≤ 1]. By inserting the nodes into the (9) or (10) with impose the initial conditions (13) or (14), we get the residuals (x) or (x) defined in respectively Solving these equations yields the coefficients c 0,j , c 1,j , . . . , c m,j . Thus, u j,m (x) and u j,m (x) given in (5) and (8) are founded.
Let us call the standard collocation method with Bernstein polynomials which yields BPSSC solution as Bernstein collocation method. Similarly, let us call the new method depending on collocation method and GBF which produces GBFSSC as GBF collocation method. The collocation nodes using in this work are the roots of Chebyshev polynomials

Residual Correction Procedure for Bernstein Collocation Method and GBF Collocation Method
Let us constitute residual correction procedure for Bernstein collocation method. We omit the residual correction procedure for GBF collocation method. By using the same method described in Section 4.1.1, we can get the coefficients c e 0,j , c e 1,j , . . . , c e m,j of e j,m (x). To do this, we construct m sets of linear or nonlinear equations such that e (x i ) = 0, i = 1, . . . , m, with the zero initial conditions. Then, e j,m (x) and hence corrected BPSSC can be obtained by solving these sets of equations.

Numerical Experiments
We demonstrate the efficiency of the present methods on five test examples. The first three examples will be solved by using the Bernstein tau method and Bernstein collocation method. The last two examples will be presented by using GBF tau method and GBF collocation method.
The exact solution set is Let us perform both method to the problem. First, we use Bernstein tau method to obtain the approximate solutions. For m = 2, approximate solutions are of the forms Now, (11) gives Additionally, we have from (13) c 0,1 = 0, c 0,2 = 1.
In Figure 1, the absolute error, the estimation of absolute error and the corrected absolute error are given for Example 1 and for m = 6. From these figures we can say that the BPSST and BPSSC solutions are well fit to the exact solutions. On the other hand, we can specify the absolute errors by using residual correction procedure. Moreover, the corrected BPSST and BPSSC solutions are more accurate. Table 1 represents the maximum absolute error for different values of m. We can say from Table 1 that increasing m yields a decreasing on the errors.

Example 2
Consider the nonlinear stiff system of ODE [3] subject to the initial conditions The exact solution is By using Theorem 1, the f = ( f 1 , f 2 ) satisfies the Lipschitz condition since f and its derivative are continuous and bound on a rectangle D. Thus, the problem has a unique solution set on an interval [0, T 0 ].
The results obtained by tau method for m = 2 are given as follows Similarly, for collocation method, the results We also obtain the approximate solutions for m = 6 and give them in Figure 2 as absolute error, corrected absolute error and corrected BPSST and BPSSC solutions to Example 2. We can see from these figures, the methods will give again more accurate results. The procedure again works well. Table 2 represents the maximum absolute error for different values of m to display the impact of m.

Example 3
Let us consider the nonlinear Genesio system [3] subject to the initial conditions where a, b and c are positive constants, satisfying ab < c. The Genesio system includes a simple square part and three simple ordinary differential equations that depend on three positive real parameters [3]. The problem has a unique solution set on [0, T 0 ] since f and its first partial derivatives are continuous and bounded on a square region. Table 3 shows the differences between the present results and that of Maple's built-in RK4. We can see that the present results agree with RK4 (step-size 0.001) at least up to 16 decimal places. Figure 3 further reconfirms the accuracy of the present solutions as compared to RK4. The maximum values of the absolute error using estimate of the absolute error on the interval [0, 1] to display the convergence of the solutions as the order m of BPSST and BPSSC are increased from 5 to 15 are given in Table 4.     the exact solution set which are non-smooth, we obtain better approximation results by GBF tau method and GBF collocation method for each problem. Institutional Review Board Statement: The study did not involve humans or animals.