Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Fibonacci, Golden Ratio, and Vector Bundles
Next Article in Special Issue
Fuzzy Numerical Solution via Finite Difference Scheme of Wave Equation in Double Parametrical Fuzzy Number Form
Previous Article in Journal
Optimal Operation of the Voltage-Doubler Boost Converter through an Evolutionary Algorithm
Previous Article in Special Issue
Double Parametric Fuzzy Numbers Approximate Scheme for Solving One-Dimensional Fuzzy Heat-Like and Wave-Like Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Enhanced Adaptive Bernstein Collocation Method for Solving Systems of ODEs

by
Ahmad Sami Bataineh
1,2,†,
Osman Rasit Isik
3,†,
Moa’ath Oqielat
1,† and
Ishak Hashim
2,*,†
1
Department of Mathematics, Faculty of Science, Al-Balqa Applied University, Al Salt 19117, Jordan
2
Department of Mathematical Sciences, Faculty of Science & Technology, Universiti Kebangsaan Malaysia, UKM, Bangi, Selangor 43600, Malaysia
3
Elementary Mathematics Education Program, Faculty of Education, Mugla Sitki Kocman University, Mugla 48000, Turkey
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2021, 9(4), 425; https://doi.org/10.3390/math9040425
Submission received: 25 January 2021 / Revised: 10 February 2021 / Accepted: 11 February 2021 / Published: 21 February 2021
(This article belongs to the Special Issue Computational Mathematics and Neural Systems)

Abstract

:
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.

1. 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 second-order 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.

2. Existence and Uniqueness Theorem

Let the function f ( x , u 1 , , u r ) be defined on the set
D = ( x , u 1 , , u r ) : a x b , < u i < for each i = 1 , 2 , , r .
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
f ( x , u 1 , , u r ) f ( x , z 1 , , z r ) L j = 1 r u j z j
for all ( x , u 1 , , u r ) , ( x , z 1 , , z r ) D .
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
f ( x , u 1 , , u r ) u i L
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
D = ( x , u 1 , , u r ) : a x b , < u i < for each i = 1 , 2 , , r .
Then, the system of (2) subject to the initial conditions has a unique solution for 0 x T 0 .

3. Bernstein Formulas and Their Operational Matrices

3.1. Bernstein Polynomials

The BP of degree m are defined by
B i , m ( x ) = m i ( x ) i ( 1 x ) m i , i = 0 , 1 , , m , x [ 0 , 1 ] ,
where the binomial coefficient is
m i = m ! i ! ( m i ) ! .
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 .

3.2. Generalized Bernstein Functions

Let us define a GBF of degree m s similar to (3) as
B ^ i , m ( x ) = m i ( x 1 s ) i ( 1 x 1 s ) m i , i = 0 , 1 , , m , s = 1 , 2 , , x [ 0 , 1 ] .
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 that B ^ i , m for 0 < i < m have to be continuous on [ 0 , ) .

3.3. 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
u j ( x ) u j , m ( x ) = C j T Φ ( x ) ,
u j ( x ) u j , m ( x ) = C j T Φ ( x ) ,
u j ( x ) u ^ j , m ( x ) = C j T Φ ^ ( x ) ,
u j ( x ) u ^ j , m ( x ) = C j T Φ ^ ( x ) ,
where C j T , Φ ( x ) , Φ ( x ) , Φ ^ ( x ) and Φ ^ ( x ) are an arbitrary ( m + 1 ) × 1 matrices defined as
C j T = [ c 0 , j , c 1 , j , , c m , j ] , Φ ( x ) = [ B 0 , m ( x ) , B 1 , m ( x ) , , B m , m ( x ) ] T , Φ ( x ) = [ B 0 , m ( x ) , B 1 , m ( x ) , , B m , m ( x ) ] T , Φ ^ ( x ) = [ B ^ 0 , m ( x ) , B ^ 1 , m ( x ) , , B ^ m , m ( x ) ] T , Φ ^ ( x ) = [ B ^ 0 , m ( x ) , B ^ 1 , m ( x ) , , B ^ m , m ( x ) ] T ,
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).

4. 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
( x ) = C j T Φ ( x ) f j x , C 1 T Φ ( x ) , , C m T Φ ( x ) .
^ ( x ) = C j T Φ ^ ( x ) f j x , C 1 T Φ ^ ( x ) , , C m T Φ ^ ( x ) .

4.1. The Approximate Solutions Obtained by Tau Method

Multiplying and integrating ( x ) and ^ ( x ) yields m equations sets
0 1 ( x ) B i , m 1 ( x ) d x = 0 , i = 0 , , m 1 .
0 1 ^ ( x ) B i , m 1 ( x ) d x = 0 , i = 0 , , m 1 .
Also, by imposing the initial conditions of Equation (2) into Equations (5) and (7) we have
u j , m ( 0 ) = C j T Φ ( 0 ) = α j ,
u ^ j , m ( 0 ) = C j T Φ ^ ( 0 ) = α j .
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, respectively. 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..

4.1.1. 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
e j ( x ) f j ( e 1 + u 1 , m , , e r + u r , m ) = u j , m
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
e j ( x ) e j , m ( x ) = C j e , T Φ ( x )
where
C j e , T = [ c 0 , j e , c 1 , j e , , c m , j e ] .
Let us define the residue e ( x ) as
e ( x ) = C j e , T Φ ( x ) f j C 1 e , T Φ ( x ) + u 1 , m , , C r e , T Φ ( x ) + u j , m + u j , m .
Then, the constants c 0 , j e , c 1 , j e , , c m , j e can be obtained by constructing m + 1 sets of equations by applying a typical tau method
0 1 e ( x ) B i , m 1 ( x ) d x = 0 , i = 0 , 1 , , m 1 ,
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 )
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
e ^ j ( x ) e ^ j , m ( x ) = C j e , T Φ ^ ( x ) .
Finally, we get the c 0 , j e , c 1 , j e , , c m , j e by solving tau method to the following system
^ e ( x ) = C j e , T Φ ^ ( x ) f j C 1 e , T Φ ^ ( x ) + u ^ 1 , m , , C r e , T Φ ^ ( x ) + u ^ r , m + u ^ j , m ( x ) ( x ) .
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 .

4.2. Approximate Solutions Obtained by Collocation Method

Let the collocation nodes be 0 x 1 < x 2 < < x m = 1 [ 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
( x i ) = C j T Φ ( x i ) f j x i , C 1 T Φ ( x i ) , , C m T Φ ( x i ) , i = 0 , 1 , 2 , , m 1 ,
^ ( x i ) = C j T Φ ^ ( x i ) f j x i , C 1 T Φ ^ ( x i ) , , C m e T Φ ^ ( x i ) , i = 0 , 1 , 2 , , m 1 .
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
x i = 1 2 + 1 2 cos ( 2 i + 1 ) π 2 m , i = 0 , 1 , , m 1 .

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 0 , j e , c 1 , j e , , c m , j e 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.

5. 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.

5.1. Example 1

Let us consider the following linear system of ODEs [7]:
u 1 ( x ) = u 1 ( x ) + u 2 ( x ) ,
u 2 ( x ) = u 1 ( x ) + u 2 ( x ) ,
with initial condition
u 1 ( 0 ) = 0 , u 2 ( 0 ) = 1 .
The exact solution set is
u 1 ( x ) = e x sin ( x ) , u 2 ( x ) = e x cos ( x ) .
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
u j , 2 ( x ) = c 0 , j B 0 , 2 ( x ) + c 1 , j B 1 , 2 ( x ) + c 2 , j B 2 , 2 ( x ) = C j T Φ ( x ) , j = 1 , 2 .
Now, (11) gives
11 12 c 0 , 1 + 1 4 c 2 , 1 + 1 6 c 1 , 1 1 4 c 0 , 2 1 6 c 1 , 2 1 12 c 2 , 2 = 0 , 11 12 c 0 , 2 + 1 4 c 2 , 2 + 1 6 c 1 , 1 + 1 12 c 2 , 1 + 1 4 c 0 , 1 + 1 6 c 1 , 2 = 0 , 5 12 c 0 , 1 + 5 12 c 2 , 1 1 2 c 1 , 1 1 12 c 0 , 2 1 6 c 1 , 2 1 4 c 2 , 2 = 0 , 5 12 c 0 , 2 + 5 12 c 2 , 2 + 1 6 c 1 , 1 + 1 4 c 2 , 1 + 1 12 c 0 , 1 1 2 c 1 , 2 = 0 .
Additionally, we have from (13)
c 0 , 1 = 0 , c 0 , 2 = 1 .
Finally by solving Equations (27) and (28) we get
c 0 , 1 = 0 , c 1 , 1 = 6 13 , c 2 , 1 = 30 13 , c 0 , 2 = 1 , c 1 , 2 = 22 13 , c 2 , 2 = 19 13 .
Thus
u 1 , 2 ( x ) u 2 , 2 ( x ) = c 0 , 1 c 1 , 1 c 2 , 1 c 0 , 2 c 1 , 2 c 2 , 2 1 2 x + x 2 2 x 2 x 2 x 2 12 13 x + 18 13 x 2 1 + 18 13 x 12 13 x 2 .
Now, let us constitute the procedure for the problem. The errors are the solutions of the following equations for j = 1 , 2
e j , 2 ( x ) = C j e , T Φ ( x ) = c 0 , j e B 0 , 2 ( x ) + c 1 , j e B 1 , 2 ( x ) + c 2 , j e B 2 , 2 ( x ) ,
or
11 12 c 0 , 1 e + 1 4 c 2 , 1 e + 1 6 c 1 , 1 e 1 4 c 0 , 2 e 1 6 c 1 , 2 e 1 12 c 2 , 2 e = 1 6 × 10 11 , 11 12 c 0 , 2 e + 1 4 c 2 , 2 e + 1 6 c 1 , 1 e + 1 12 c 2 , 1 e + 1 4 c 0 , 1 e + 1 6 c 1 , 2 e = 1 3 × 10 10 , 5 12 c 0 , 1 e + 5 12 c 2 , 1 e 1 2 c 1 , 1 e 1 12 c 0 , 2 e 1 6 c 1 , 2 e 1 4 c 2 , 2 e = 1 6 × 10 11 , 5 12 c 0 , 2 e + 5 12 c 2 , 2 e + 1 6 c 1 , 1 e + 1 4 c 2 , 1 e + 1 12 c 0 , 1 e 1 2 c 1 , 2 e = 1 3 × 10 10 , c 0 , 1 e = 0 , c 0 , 2 e = 0 .
The last row comes from (16). Therefore, the estimations of the errors are found as
c 0 , 1 e = 0 , c 1 , 1 e = 0.5769230763 × 10 11 , c 2 , 1 e = 0.6615384615 × 10 10 , c 0 , 2 e = 0 , c 1 , 2 e = 0.3884615384 × 10 10 , c 2 , 2 e = 0.8923076921 × 10 10 ,
that is e 1 , 2 ( x ) e 2 , 2 ( x ) 0.1153846153 × 10 10 x + 0.7769230768 × 10 10 x 2 0.7769230768 × 10 10 x + 0.1153846153 × 10 10 x 2 .
Second, we will find the approximate solutions by using Bernstein collocation method. If the steps in Section 4.2 are performed to the problem, the BPSSC solution set of Equations (24)–(26) will be obtained as
u 1 , 2 ( x ) u 2 , 2 ( x ) 0.96 x + 1.28 x 2 1 + 1.28 x 0.96 x 2 .
A similar argument to the first case yields the estimations of the errors as
e 1 , 2 ( x ) e 2 , 2 ( x ) 0.1500553308 × 10 9 x + 0.1157274735 × 10 8 x 2 0.1786389502 × 10 9 x + 0.26797339 × 10 9 x 2 .
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.

5.2. Example 2

Consider the nonlinear stiff system of ODE [3]
u 1 ( x ) = 1002 u 1 ( x ) + 1000 u 2 2 ( x ) ,
u 2 ( x ) = u 1 ( x ) u 2 ( x ) u 2 2 ( x ) ,
subject to the initial conditions
u 1 ( 0 ) = 1 , u 2 ( 0 ) = 1 .
The exact solution is
u 1 ( x ) = e 2 x , u 2 ( x ) = e x .
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
u 1 , 2 ( x ) u 2 , 2 ( x ) 1 1.75 x + 0.93 x 2 1 0.95 x + 0.31 x 2 ,
e 1 , 2 ( x ) e 2 , 2 ( x ) 0.19 × 10 8 x + 0.25 × 10 8 x 2 0.36 × 10 9 x + 0.63 × 10 9 x 2 .
Similarly, for collocation method, the results
u 1 , 2 ( x ) u 2 , 2 ( x ) 1 1.85 x + 1.03 x 2 1 0.95 x + 0.31 x 2 ,
e 1 , 2 ( x ) e 2 , 2 ( x ) 0.11 × 10 8 x + 0.26 × 10 8 x 2 0.16 × 10 9 x + 0.36 × 10 10 x 2 .
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.

5.3. Example 3

Let us consider the nonlinear Genesio system [3]
u 1 ( x ) = u 2 ( x ) ,
u 2 ( x ) = u 3 ( x ) ,
u 3 ( x ) = c u 1 ( x ) b u 2 ( x ) a u 3 ( x ) + u 1 ( x ) 2 ,
subject to the initial conditions
u 1 ( 0 ) = 0.2 u 2 ( 0 ) = 0.3 , u 3 ( 0 ) = 0.1 ,
where a, b and c are positive constants, satisfying a b < 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.

5.4. Example 4

Let us consider the following nonhomogeneous linear systems of ODEs:
u 1 ( x ) u 1 ( x ) u 2 ( x ) = 2 x 3 + 1 3 x 4 3 3 x 3 x 3 + 3 x 11 3 3 x 2 3 ,
u 2 ( x ) + u 1 ( x ) u 2 ( x ) = 7 x 4 3 3 3 x 2 + x 2 3 + x 3 x 7 3 + x 3 ,
with initial condition
u 1 ( 0 ) = 0 , u 2 ( 0 ) = 0 .
The exact solution set is
u 1 ( x ) = x 2 3 + x 3 , u 2 ( x ) = x 7 3 x 3 .
Let us perform the GBF tau method and GBF collocation method to obtain the approximate solutions. Now for m = 9 and s = 3 , approximate solutions are of the forms
u ^ j , 9 ( x ) = i = 0 m c i , j B ^ i , m ( x ) , j = 1 , 2 .
Now, (12) gives
333 665 c 0 , 1 + 1 1 , 175 , 720 c 9 , 2 = 1 , 510 , 963 12 , 932 , 920 1 1330 c 0 , 1 + + 3 152 , 152 c 9 , 2 = 36 , 551 25 , 865 , 840 , 9 1 , 293 , 292 c 0 , 1 + + 477 2380 c 9 , 2 = 55 , 089 180 , 880 .
Also, we have from (14)
c 0 , 1 = 0 , c 0 , 2 = .
Finally by solving Equations (39) and (40) we get
c 0 , 1 = 0 , c 1 , 1 = 1 9 , c 2 , 1 = 1 4 , c 3 , 1 = 5 12 , c 4 , 1 = 11 18 , c 5 , 1 = 5 6 , c 6 , 1 = 13 12 , c 7 , 1 = 49 36 , c 8 , 1 = 5 3 , c 9 , 1 = 2 , c 0 , 2 = 0 , c 1 , 2 = 0 , c 2 , 2 = 0 , c 3 , 2 = 0 , c 4 , 2 = 0 , c 5 , 2 = 0 , c 6 , 2 = 0 , c 7 , 2 = 1 36 , c 8 , 2 = 2 9 , c 9 , 2 = 0 .
Thus
u ^ 1 , 2 ( x ) u ^ 2 , 2 ( x ) = x 2 3 + x 3 x 7 3 x 3 .
which is the exact solution (38). Note that the exact solution (38) can be obtained for any m 9 .

5.5. Example 5

Finally we consider the non-homogeneous nonlinear systems of ODEs:
u 1 ( x ) 1002 u 1 ( x ) 1000 u 2 2 ( x ) = 1 4 x 6004 x + 2000 x 3 2 x ,
u 2 ( x ) u 1 ( x ) + u 2 ( x ) + u 2 2 ( x ) = 1 8 x + 2 x + 2 x 3 2 x ,
subject to the initial conditions
u 1 ( 0 ) = 1 , u 2 ( 0 ) = 1 .
The exact solution is
u 1 ( x ) = 1 + x , u 2 ( x ) = 1 x .
The results obtained by tau method for m = 3 and s = 2 are given as follows
u ^ 1 , 3 ( x ) u ^ 2 , 3 ( x ) 1 + x 0.3 × 10 17 x + 0.31 × 10 17 x 3 1 x 0.2 × 10 18 x + 0.3 × 10 18 x 3 ,
which is almost the exact solution (43).

6. Conclusions

We proposed two methods to numerically solve systems of ODEs. These methods are direct methods. We gave a detailed error analysis for both methods. We also found an upper bound of the absolute error for collocation method. As seen from the numerical examples, the methods give more accurate approximate solutions for linear and nonlinear cases of the problem with the stiff problem. Increasing the number of nodes yields a decrease of the absolute errors. Residual correction procedure estimates the errors for Example 1 and Example 2 with high accuracy. The corrected approximate solutions are better than the BPSST and BPSSC. On the other hand, the results are consistent with the results of RK4 method for Example 3. Even though the last two nonlinear problems have the exact solution set which are non-smooth, we obtain better approximation results by GBF tau method and GBF collocation method for each problem.

Author Contributions

Conceptualization, A.S.B., O.R.I., M.O. and I.H.; methodology, A.S.B. and O.R.I.; software, A.S.B., O.R.I., and M.O.; validation, A.S.B., O.R.I., M.O. and I.H.; formal analysis, A.S.B. and O.R.I.; writing—original draft preparation, A.S.B., O.R.I., M.O. and I.H.; writing—review and editing, A.S.B., O.R.I., M.O. and I.H.; funding acquisition, I.H. All authors have read and agreed to the published version of the manuscript.

Funding

The APC was funded by I.H.’s Universiti Kebangsaan Grant # GP-2019-6388.

Institutional Review Board Statement

The study did not involve humans or animals.

Informed Consent Statement

Not applicable.

Data Availability Statement

The study did not report any data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Shawagfeh, N.; Kaya, D. Comparing numerical methods for the solutions of systems of ordinary differential equations. Appl. Math. Lett. 2004, 17, 323–328. [Google Scholar] [CrossRef] [Green Version]
  2. Hairer, E.; Wanner, G. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems; Springer: Berlin/Heidelberg, Germany, 1991. [Google Scholar]
  3. Bataineh, A.S.; Noorani, M.S.M.; Hashim, I. Solving systems of ODEs by homotopy analysis method. Commun. Nonlinear Sci. Numer. Simul. 2008, 13, 2060–2070. [Google Scholar] [CrossRef]
  4. Park, J.H. GCS of a class of chaotic dynamic systems. Chaos Solitons Fract. 2005, 26, 1429–1435. [Google Scholar] [CrossRef]
  5. Park, J.H. Chaos synchronization between two different chaotic dynamical systems. Chaos Solitons Fract. 2006, 27, 549–554. [Google Scholar] [CrossRef]
  6. Genesio, R.; Tesi, A. A harmonic balance methods for the analysis of chaotic dynamics in nonlinear systems. Automatica 1992, 28, 531–548. [Google Scholar] [CrossRef]
  7. Hashim, I.; Chowdhury, M.S.H. Adaptation of homotopy-perturbation method for numeric-analytic solution of system of ODEs. Phys. Lett. A 2008, 372, 470–481. [Google Scholar] [CrossRef]
  8. Isik, O.R.; Sezer, M.; Guney, Z. A rational approximation based on Bernstein polynomials for high order initial and boundary values problems. Appl. Math. Comput. 2011, 217, 9438–9450. [Google Scholar]
  9. Rostamy, D.; Karimi, K. Bernstein polynomials for solving fractional heat- and wave-like equations. Fract. Calc. Appl. Anal. 2012, 15, 556–571. [Google Scholar] [CrossRef]
  10. Yuzbasi, S. Numerical solutions of fractional Riccati type differential equations by means of the Bernstein polynomials. Appl. Math. Comput. 2013, 219, 6328–6343. [Google Scholar]
  11. Baleanu, D.; Alipour, M.; Jafari, H. The Bernstein operational matrices for solving the fractional quadratic Riccati differential equations with the Riemann-Liouville derivative. Abstr. Appl. Anal. 2013, 2013, 461970. [Google Scholar] [CrossRef]
  12. Pandey, R.K.; Kumar, N. Solution of Lane-Emden type equations using Bernstein operational matrix of differentiation. New Astrom. 2012, 17, 303–308. [Google Scholar] [CrossRef]
  13. Isik, O.R.; Sezer, M. Bernstein series solution of a class of Lane-Emden type equations. Math. Probl. Eng. 2013, 2013, 423797. [Google Scholar] [CrossRef]
  14. Isik, O.R.; Sezer, M.; Guney, Z. Bernstein series solution of linear second-order partial differential equations with mixed conditions. Math. Meth. Appl. Sci. 2014, 37, 609–619. [Google Scholar] [CrossRef]
  15. Maleknejad, K.; Basirat, B.; Hashemizadeh, E. A Bernstein operational matrix approach for solving a system of high order linear Volterra-Fredholm integro-differential equations. Math. Comput. Model. 2012, 55, 1363–1372. [Google Scholar] [CrossRef]
  16. Rostamy, D.; Karimi, K. A new operational matrix method based on the Bernstein polynomials for solving the backward inverse heat conduction problems. Int. J. Numer. Meth. Heat Fluid Flow 2014, 24, 669–678. [Google Scholar] [CrossRef]
  17. Alshbool, M.H.T.; Hashim, I. Multistage Bernstein polynomials for the solutions of the fractional order stiff systems. J. King Saud Univ. Sci. 2016, 28, 280–285. [Google Scholar] [CrossRef] [Green Version]
  18. Alshbool, M.H.T.; Bataineh, A.S.; Hashim, I.; Isik, O.R. Solution of fractional-order differential equations based on the operational matrices of new fractional Bernstein functions. J. King Saud Univ. Sci. 2017, 29, 1–18. [Google Scholar] [CrossRef] [Green Version]
  19. Asgari, M.; Ezzati, R. Using operational matrix of two-dimensional Bernstein polynomials for solving two-dimensional integral equations of fractional order. Appl. Math. Comput. 2017, 307, 290–298. [Google Scholar] [CrossRef]
  20. Bataineh, A.S.; Alomari, A.K.; Isik, O.R.; Hashim, I. Multistage Bernstein collocation method for solving strongly nonlinear damped systems. J. Vib. Control 2019, 25, 122–131. [Google Scholar] [CrossRef]
  21. Khataybeh, S.; Hashim, I.; Alshbool, M. Solving directly third-order ODEs using operational matrices of Bernstein polynomials method with applications to fluid flow equations. J. King Saud Univ. Sci. 2019, 31, 822–826. [Google Scholar] [CrossRef]
  22. Khataybeh, S.; Hashim, I. Direct solution of second-order system of ODEs using Bernstein polynomials on an arbitrary interval. Int. J. Math. Comput. Sci. 2019, 14, 343–357. [Google Scholar]
  23. Bataineh, A.S.; Isik, O.R.; Alomari, A.K.; Shatnawi, W.; Hashim, I. An efficient scheme for time-dependent Emden-Fowler type equations based on two-dimensional Bernstein polynomials. Mathematics 2020, 8, 1473. [Google Scholar] [CrossRef]
  24. Alshbool, M.H.T.; Shatanawi, W.; Hashim, I.; Sarr, M. Residual correction procedure with Bernstein polynomials for solving important systems of ordinary differential equations. Comput. Math. Contin. 2020, 64, 63–80. [Google Scholar] [CrossRef]
  25. Alshbool, M.H.; Isik, O.; Hashim, I. Fractional Bernstein series solution of fractional diffusion equations with error estimate. Axioms 2021, 10, 6. [Google Scholar] [CrossRef]
  26. Burden, R.L.; Faires, J.D. Numerical Analysis, 9th ed.; Brooks/Cole: Boston, MA, USA, 2007. [Google Scholar]
  27. Birkhoff, G.; Rota, G.-C. Ordinary Differential Equations, 4th ed.; John Wiley & Sons: New York, NY, USA, 1989. [Google Scholar]
Figure 1. The absolute error, estimation of absolute error and the corrected absolute error to Example 1 for m = 6 .
Figure 1. The absolute error, estimation of absolute error and the corrected absolute error to Example 1 for m = 6 .
Mathematics 09 00425 g001aMathematics 09 00425 g001b
Figure 2. The absolute error, the corrected absolute error and corrected approximate solutions to Example 2 for m = 6 .
Figure 2. The absolute error, the corrected absolute error and corrected approximate solutions to Example 2 for m = 6 .
Mathematics 09 00425 g002
Figure 3. The corrected approximate solutions to Example 3 for m = 6 .
Figure 3. The corrected approximate solutions to Example 3 for m = 6 .
Mathematics 09 00425 g003
Table 1. The maximum absolute errors on the interval [ 0 , 1 ] for different values of m for different m values and Example 1.
Table 1. The maximum absolute errors on the interval [ 0 , 1 ] for different values of m for different m values and Example 1.
Methodm51015
Tau method u 1 u 1 , m 1.2 × 10 5 3.5 × 10 13 6.6 × 10 21
Tau method u 2 u 2 , m 6.8 × 10 6 1.3 × 10 12 1.2 × 10 20
Coll. method u 1 u 1 , m 2.0 × 10 5 6.8 × 10 13 1.1 × 10 20
Coll. method u 2 u 2 , m 1.2 × 10 5 2.2 × 10 12 1.9 × 10 20
Table 2. The maximum absolute errors on the interval [ 0 , 1 ] for different values of m and Example 2.
Table 2. The maximum absolute errors on the interval [ 0 , 1 ] for different values of m and Example 2.
Methodm51015
Tau method u 1 u 1 , m 6.9 × 10 5 4.8 × 10 11 7.2 × 10 16
Tau method u 2 u 2 , m 6.4 × 10 7 4.8 × 10 14 3.3 × 10 16
Coll. method u 1 u 1 , m 6.1 × 10 5 3.5 × 10 11 8.1 × 10 16
Coll. method u 2 u 2 , m 1.0 × 10 6 4.3 × 10 14 3.3 × 10 16
Table 3. Differences between Bernstein series solution and RK4 solutions in the case a = 1.2 , b = 2.92 , c = 6 , for i = 1 , 2 , 3 .
Table 3. Differences between Bernstein series solution and RK4 solutions in the case a = 1.2 , b = 2.92 , c = 6 , for i = 1 , 2 , 3 .
x Δ = | u i , 6 RK 4 0.001 | Δ = | u i , 6 + e i , 29 RK 4 0.001 |
Δ u 1 Δ u 2 Δ u 3 Δ u 1 Δ u 2 Δ u 3
0.0000000
0.10.14 × 10 6 0.50 × 10 6 0.23 × 10 6 0.15 × 10 18 0.13 × 10 18 0.21 × 10 18
0.20.18 × 10 6 0.42 × 10 6 0.46 × 10 6 0.33 × 10 18 0.32 × 10 18 0.67 × 10 18
0.30.12 × 10 7 0.42 × 10 6 0.19 × 10 6 0.56 × 10 18 0.50 × 10 18 0.30 × 10 18
0.40.13 × 10 6 0.64 × 10 6 0.11 × 10 6 0.81 × 10 18 0.65 × 10 18 0.20 × 10 18
0.50.61 × 10 8 0.50 × 10 7 0.30 × 10 7 0.11 × 10 17 0.73 × 10 18 0.11 × 10 17
0.60.22 × 10 6 0.73 × 10 6 0.30 × 10 6 0.12 × 10 17 0.67 × 10 18 0.19 × 10 17
0.70.23 × 10 6 0.46 × 10 6 0.60 × 10 6 0.17 × 10 17 0.48 × 10 18 0.32 × 10 17
0.80.39 × 10 7 0.40 × 10 6 0.40 × 10 6 0.19 × 10 17 0.10 × 10 18 0.45 × 10 17
0.90.24 × 10 7 0.47 × 10 6 0.20 × 10 6 0.22 × 10 17 0.47 × 10 18 0.58 × 10 17
1.00.87 × 10 7 0.10 × 10 7 0.30 × 10 6 0.23 × 10 17 0.12 × 10 17 0.70 × 10 17
Table 4. The maximum values of the absolute errors by using the estimations of the absolute errors on the interval [ 0 , 1 ] for Example 3.
Table 4. The maximum values of the absolute errors by using the estimations of the absolute errors on the interval [ 0 , 1 ] for Example 3.
Methodm51015
Tau method e 1 , m 1.8 × 10 6 1.2 × 10 12 4.1 × 10 18
Tau method e 2 , m 3.6 × 10 6 7.8 × 10 12 1.6 × 10 18
Tau method e 3 , m 1.4 × 10 5 3.8 × 10 11 4.5 × 10 17
Coll. method e 1 , m 3.4 × 10 6 1.9 × 10 12 6.2 × 10 18
Coll. method e 2 , m 7.2 × 10 6 1.2 × 10 11 4.1 × 10 18
Coll. method e 3 , m 2.3 × 10 5 6.1 × 10 11 6.8 × 10 17
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bataineh, A.S.; Isik, O.R.; Oqielat, M.; Hashim, I. An Enhanced Adaptive Bernstein Collocation Method for Solving Systems of ODEs. Mathematics 2021, 9, 425. https://doi.org/10.3390/math9040425

AMA Style

Bataineh AS, Isik OR, Oqielat M, Hashim I. An Enhanced Adaptive Bernstein Collocation Method for Solving Systems of ODEs. Mathematics. 2021; 9(4):425. https://doi.org/10.3390/math9040425

Chicago/Turabian Style

Bataineh, Ahmad Sami, Osman Rasit Isik, Moa’ath Oqielat, and Ishak Hashim. 2021. "An Enhanced Adaptive Bernstein Collocation Method for Solving Systems of ODEs" Mathematics 9, no. 4: 425. https://doi.org/10.3390/math9040425

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop