Numerical Solutions of High-Order Differential Equations with Polynomial Coefficients Using a Bernstein Polynomial Basis

The paper presents a novel method that allows one to establish numerical solutions of linear and nonlinear ordinary differential equations—with polynomial coefficients—that contain any finite products of the unknown functions and/or their general derivatives. The presented algorithm provides numerical solutions of these differential equations subject to initial or boundary conditions. This algorithm proposes the desired solution in terms of B-polynomials (Bernstein polynomial basis) and then uses the orthonormal relation of B-polynomials with its weighted dual basis with respect to the Jacobi weight function to construct a linear/nonlinear system in the unknown expansion coefficients which can be solved using a suitable solver. The properties of B-polynomials provide greater flexibility in which to impose the initial or boundary conditions at the end points of the interval [0, R] and enable us to obtain exactly and explicitly some of the unknown expansion coefficients in the form of a suggested numerical solution. Consequently, the presented algorithm leads to a linear or nonlinear algebraic system in the unknown expansion coefficients that has a simpler form than that was obtained by the other algorithms. So that, this procedure is a powerful tool that we may utilize to overcome the difficulties associated with boundary and initial value problems with less computational effort than the other techniques. An accepted agreement is obtained between the exact and approximate solutions for the given examples. The error analysis was also studied, and the obtained numerical results clarified the validity of the theoretical results.


Introduction
A boundary-value problem for an ordinary differential equation (or system of equations) is obtained by requiring that the dependent variable (or variables) is satisfying subsidiary conditions at one or more distinct points.If these conditions are imposed at one or two points, the existence and uniqueness theory for initial/boundary-value problems can be expanded for many special classes of equations and systems of equations.Many approximate methods for obtaining numerical solutions for various classes of differential equations using orthogonal polynomials are available in the literature (see, for instance, [4][5][6]50]).
A great interest to use B-polynomials for solving different classes of differential equations (see, for instance, [7][8][9]12,16,35,36,[42][43][44][45]) due to the perfect properties like the recursive relation, the appropriate properties, and making partition of unity (see, for instance, [17,[24][25][26]28,40]).These advantages make them the most known used basis in approximation theory and computer-aided geometric design (CAGD) (see, [23,32,37]).Farouki and GoodMan [26] show that there is no other nonnegative basis that provides systematically smaller condition numbers than B-polynomials, in the sense that they have optimal stability.Also, they declare that although B-polynomials are not uniquely optimal, no other basis in popular use provides this distinction, and it is uncertain whether other optimally stable bases would present the advantages and algorithms that we associate with the Bernstein form.Farouki and Rajan [29] had a head start in developing Bernstein formulations for all the procedures of basic polynomials required in geometric modeling algorithms, and they declared that these formulations were as simple as their customary power formulations.
In the current paper, a new numerical algorithm to find numerical solutions of higher order linear and nonlinear ordinary differential equationswith polynomial coefficients-that contain any finite products of the unknown functions and/or their derivatives subject to initial and boundary conditions will be discussed.Our interest in such models stems from the fact that they appear as models of many practical problems in mechanics and other areas of mathematical physics.For instance, the study of hydrodynamics, hydromagnetics, fluid dynamics, astrophysics, astronomy, beam and long wave theory, and applied physics.To be precise, for example, fifth-order boundary-value problems (BVPs) can be used to model viscoelastic flows (see, [19,34]), sixthorder BVPs arise in the narrow convecting layers bounded by stable layers, which are believed to surround A-type stars [10,14,22,48], the seventh-order BVPs generally arise in modeling induction motors with two rotor circuits [41] and ordinary convection and overstability yield a tenth-order and a twelfthorder BVPs, respectively.
The suggested numerical algorithm approximates the solution in the B-polynomials and it is based on the orthonormal relation of B-polynomials with its weighted dual basis with respect to the Jacobi weight function to construct a linear/nonlinear system in the unknown expansion coefficients which can be solved using a suitable solver.It will also be demonstrated that this algorithm provides excellent agreement between exact and approximate solutions, with the maximum pointwise error of no more than O (10 −20 ).In this paper, we discuss formulae related to B-polynomials, which are useful in many ways.These formulae enable us with great flexibility to impose boundary conditions at x = 0, R or initial conditions at x = 0. Thus, the contribution of this paper can be summarized in the following items: • A new explicit formula expressing the general-order derivative of Bpolynomial in terms of the B-polynomials themselves is proved (see Lemma 2).• A formula for the B-polynomial coefficients of the moments of one single or a general-order derivative of B-polynomial is given (see Theorem 1).• An explicit formula, which expresses the B-polynomial coefficients of a general-order derivative of any polynomial in terms of its original Bpolynomial coefficients, is also proved (see Theorem 2).• Establishing an expression of the B-polynomial coefficients of the finite products of polynomials and/or their general-order derivatives in terms of its original B-polynomial coefficients is given (see Theorem 3 and Corollary 3).
• Establishing an expression of the B-polynomial coefficients of the moments of finite products of polynomials and/or their general-order derivatives in terms of its original B-polynomial coefficients is given (see Corollary 2).• Designing an approach based on the orthonormal relation of B-polynomials with its weighted dual basis with respect to the Jacobi weight function to treat two types of high-order initial and initialboundary problems.• Investigating the error analysis of the B-polynomial expansion.
• Examining the accuracy of the proposed method via the presentation of some illustrative examples.
Up to now, and to the best of our Knowledge, it is worthy to state that the proven formulae in Theorems 1-3, two Corollaries 2 and 3 and Lemma 2 for those previously mentioned are new and traceless in the literature.These formulae lead to the systematic character and simplicity of the proposed algorithm, which allow it to be implemented in any computer algebra (here, the Mathematica 12) symbolic language used.Also, the proposed algorithm provide exactly and explicitly some of the unknown expansion coefficients (see Eqs. (54) and (61)) in the form of a suggested numerical solution.Consequently, the presented algorithm leads to a linear or nonlinear algebraic system with unknown expansion coefficients that has a simpler form than that obtained by the other algorithms.So that, this procedure is a powerful tool that we may utilize to overcome the difficulties associated with boundary and initial value problems with less computational effort than the other techniques.

MJOM
The current paper is organized as follows: In Sect.2, computational implementations which have useful uses in introducing the proposed algorithm are given.In Sects. 3 and 4, the details of the proposed numerical method to solve linear and nonlinear differential equations with polynomial coefficients are provided.The error analysis is presented in Sect. 5. Numerical examples are given in Sect.6, and the obtained numerical results illustrate the validity of the theoretical results in Sect. 5.

Polynomial Basis
The general form of the B-polynomials of nth degree is defined by [12, pp.273-274] which constitute the entire basis of (n + 1) nth-degree polynomials.For convenience, we set A recursive definition can also be used to generate the B-polynomials over this interval, allowing us to write the ith nth-degree B-polynomial in the form The derivative of B i,n (x) is given by Each B-polynomial is positive and the sum of all the B-polynomials is unity for all x ∈ [0, R], that is The set of functions which has these properties is called a partition of unity on the interval [0, R].In addition, any given polynomial P (x) of degree n can be expanded as follows [27, p.3]: The dual basis to the Bernstein polynomials with respect to the Jacobi weight function w is characterized by the property [39, p.1587] where the dual basis functions d k,n (u) have explicit representations and the coefficients c j,k (n) are defined by where Therefore, the dual basis, (9) and has the explicit representations Proof.Using relation (4), we have then using the orthogonality relation (9), one can see that (11) is valid, and the proof of Lemma 1 is complete.
Remark 1. Formula ( 11) is a generalization for the result of Jani et al. [33, p.7666] Lemma 2. Let B i,n (x) be the ith nth-degree B-polynomial, and then, the pthderivatives can be written in the form H. M. Ahmed

MJOM
where λ Proof.To prove (12), we proceed by induction.In view of relation (3), we may write and this in turn shows that ( 12) is true for p = 1.Proceeding by induction, assuming that ( 12) is valid for p, we want to prove that From (3) and assuming the validity for p, we have Collecting similar terms and using the relations we get (13) and the proof of Lemma 2 is complete.

Note 1.
Since B i,n (x) = 0, if i < 0 or i > n, then formula (12) takes the form In particular, and for the special case R = 1, this formula coincides formula (3.1) in [20].
Corollary 1.Let B i,n (x) be the ith nth-degree B-polynomial, then it is not difficult to show that where (a) n = Γ(a + n)/Γ(a) is the pochhammer symbol.
Theorem 1.Let B i,n (x) be the ith nth degree B-polynomial, then Proof.Using Corollary 1 and Lemma 2, we obtain where By (n−k)-fold degree elevation (see, [28]), we can express each Bernstein basis function of degree k in the Bernstein basis of degree n as where Then, by the aid of (20), formula (18) can be written in the form Collecting similar terms and using B i,n (x) = 0,if i < 0 or i > n, then formula (22) takes the form Substitution of ( 19) and ( 21) into ( 23) gives (16), and this completes the proof of Theorem 1.
According to Theorem 1, we can state the following theorem which relates the B-polynomial basis coefficients of the moments of a general-order p-derivative of a polynomial f (x) in terms of its B-polynomial basis coefficients.

MJOM
Theorem 2. Assume that a polynomial f (x) of degree N at most, and the moments of a general-order p-derivative of f (x), x m f (p) (x), are formally expanded in a finite series of B-polynomial basis and then Proof.Using relation (16) gives Collecting similar terms, we get which completes the proof of Theorem 2.
Note 2. It is to be noted here that and then Also, it is worth to note that In particular, and for the special case m = r = 0, and using relation ( 29)-and after some rather manipulation-formula (26) takes the form (31) This shows that Theorem 3.3 in [20] is a direct consequence of Theorem 2 for R = 1.
Farouki and Rajan [29, pp.13-14] show that if the two functions f (x) and g(x) have the form then the product f (x)g(x) may be expressed as In view of this result and Theorem 2, the following theorem can be proved.
Theorem 3. Assume that the polynomials f (x) of degrees N ( = 1, 2, . . ., s) at most, respectively, are formally expanded in finite series of B-polynomial basis where C (p1,p2,...,ps) js . . . where Proof.In view of Theorem 2, we have By repeating the application of formula (32)-and after some rather manipulation-one can obtain (34).
It is worthy to note that formula (34) for the case of s = 1 can be written as where the coefficients C (p1) j1 can be defined as In view of two Theorems 2 and 3, we obtain the following corollary.
Corollary 2. Under the assumptions of two Theorems 2 and 3, we get where a (p1,p2,...,ps) js In view of Note 2, formula (32) can be generalized as a direct consequence of Corollary 2 as in the following corollary.

An Application for the Solution of High-Order Linear Differential Equations with Polynomial Coefficients
In this section, we aim to discuss an algorithm for approximating solutions to qth-order ordinary linear differential equations subject to the boundary conditions where q 1 + q 2 + 2 = q, or the initial conditions where p i (x), i = 0, 1, ..., q, are given polynomials.Without loss of any generality, we suppose that p i (x) = γ i x mi , i = 0, 1, ..., q, where m i are positive integers and γ i are constants, f (x) is a given source function and α j (j = 0, 1, ..., q 1 ), β j ( j = 0, 1, ..., q 2 ) are constants.An approximation to the solution of ( 43) may be written as In the case of boundary-value problem (43) and (44): We make y N (x) satisfies the following equations: for n = q 1 + 1, ..., N − q 2 − 1, and In this case, the coefficients a 0 , ..., a q1 , a N −q2 , ..., a N can be determined as follows: According to relation (12) and j <i≤ N, and then Eq. ( 49) takes the form and and It is not difficult to show that Now, the coefficients a q1+1 , ..., a N −q2−1 can be obtained as follows: In view of Theorem 2, we obtain then substituting (55) in (47) and using the formula (26) and the orthonormal relation where Any appearance of a 0 , ..., a q1 , a N −q2 , ..., a N in (57) is moved to the right-hand side, then the system of equations ( 57) is equivalent to the following matrix equation: where a

and the elements of
respectively.
In the case of the initial value problem (43) and (45): In similar way, it can be shown that the coefficients a 0 , ..., a q−1 have the form and the coefficients a q , ..., a N can be obtained by making y N (x) satisfies (47) for n = q, ..., N.Then, the obtained equations can be written in the matrix form (58) where a (q) N = [a q , ..., a N ] T , and the elements of T have the forms (59) and Note 4. It is worthy to note that, if the polynomials p j (x) (j = 0, 1, . . ., q) have the form N have the form

An Application for the Solution of High-Order Nonlinear Differential Equations with Polynomial Coefficients
In this section, we are interested in solving numerically the nonlinear differential equation with polynomial coefficients subject to the boundary conditions (44) or the initial conditions (45), where Assume the proposed approximated solution y N (x) as in (46).

In the case of initial value problem (63) and (45):
In similar way, the coefficients a 0 , ..., a q−1 have the form (61).While the coefficients a q , ..., a N can be obtained by making y N (x) satisfy (64) for n = q, ..., N, which can be written in the form of nonlinear equations (66) for n = q, ..., N.Then, these nonlinear equations can be solved using an appropriate solver to obtain these coefficients.

Error Analysis
In this section, we provide the studying of error resulting from the presented numerical method.For a positive integer N, consider the space S N defined by Let P N denote the orthogonal projection of the ω-weighted L 2 onto S N .This projection has the following approximation properties, whose proof is standard.

Proposition 1. There exists a constant C independent of N , such that for any
Lemma 3. The approximated solutions y N (x), in the form (46), of the differential Eq. ( 43) with the boundary conditions (44) or the initial conditions (45) satisfy the inequality Proof.We have Also, we have In view of relation ( 11), we get max From ( 69) and (70), we obtain (68).
The following lemma proves the stability of the presented solution (As an analog of the results in [8]).
Lemma 4. The obtained approximated solutions y N (x) satisfy the following two inequalities: where h Hence, for sufficient large N , we have Proof.If f is a function defined on [0, 1], then the Bernstein polynomial B N (f, x) of f is given by converges to f (x) uniformly on [0, 1] if f is continuous there [11,18].Regarding the rate of converges to f (x), Popviciu [38] has shown that 303 Page 16 of 31

MJOM
where ω f is the modulus of continuity of f in [0, 1], which have the form ω f (N −1/2 ) = C N −1/2 .Without loss of generality, we can consider [0, R] = [0, 1], then we obtain and Now, using (77), one can obtain then using ( 76) and (78) give From (79) and using that then relation (71) is proved, and this implies that for sufficient large N , we have and this completes the proof of (73).Also, we have then using that From (81) and using two relations (68) and (71) enable one to get then (72) is proved, and hence, for sufficient large N , we have which completes the proof of (74) and the proof of Lemma 4 is complete.

Numerical Solutions of High-Order Differential
Page 17 of 31 303 Theorem 4. Let y(x) ∈ H r [0, R], for some integer r ≥ 2, and y N (x) be the obtained approximate solution, then we have where λ is a constant independent of N. Hence, for sufficient large N , we have Proof.Let P N denote an orthogonal projection of the ω -weighted L 2 onto S N .We can define it as follows: where In view of Proposition 1, we get (87) However then the two inequalities (87) and (88) lead to The inequality (89) can be written as where λ = max(C1,C2) 1−M is a constant independent of N , then (82) is proved.Hence, for sufficient large N , (74) and (82) give However, we have (91) then using (90) and (91) lead to (83).Then, (84) and ( 85) are obtained directly from (83) which completes the proof of this theorem.The following two corollaries ensure that the stability of the solutions as N increases (see [1][2][3]8]).
Corollary 5.Under the assumptions of Theorem 4, we have for sufficient large N, and for minimum integer r ≥ 2, such that ( 83) is valid, we have and Proof.We have In view of (83), we get (92).Then, (93) is a direct consequence of (92).
Corollary 6.Under the assumptions of Theorem 4, we have for sufficient large N, and for minimum integer r ≥ 2, such that ( 83) is valid, we have and Proof.We have In view of Corollary 5, we get (94).Then, (95) is a direct consequence of (94).

Numerical Results
In this section, we use the algorithms presented in Sects.3 and 4 to solve several numerical examples.All computations presented were performed on a computer (Intel(R) Core(TM) i9-10850 CPU@ 3.60GHz, 3600 Mhz, 10 Core(s), 20 Logical Processor (s)) running Mathematica 12 .Comparisons between the maximum absolute errors obtained by the present method and other methods proposed in [20,21,31,46,47,49] are made.In view of Theorem 4, we have where is the upper bound of absolute errors.
Example 1.Consider the following boundary-value problem to demonstrate that the B-polynomials are powerful to approximate the solution to desired accuracy: whose exact solution is y(x) = e x .Computational results for E N ∞ , with various choices of α, β and N , are presented in Table 1 with displaying the computational time (CT(s)), in seconds, to get these results.In Table 2, a comparison between the error obtained using the presented method (using α = β = 2), the sinc-Galerkin, modified decomposition methods, BGM, and PBGM (see, [20,21]) is displayed which shows improved performance in the presented method.Figure 1a shows that the obtained absolute errors E 12 (x) and E 13 (x) at α = 2 and β = 2, to demonstrate the convergence of solutions.Additionally, Fig. 1b shows an excellent agreement between the exact and approximate solution y 14 (x) to the given problem (98).
Example 2. Consider the following seventh-order linear boundary-value problem: y(1) = 0, y (1) (1) = −e, y (2) (1) = −2e, ( whose exact solution is given by y(x) = (1 − x)e x .Computational results for E N ∞ -error norm with various choices of α, β and N, are presented  3.73140×10 −13 8.6685×10 −16   in Table 5.A comparison between the absolute errors obtained by using the present method (using N = 17, α = 0, and β = 0) and the methods in [46,47] at various values of x, is given in Table 6.This shows that the present method is more accurate.Figure 2a shows Log-errors for various N values, when α = 0 and β = 0, to demonstrate the stability of solutions.Additionally, Fig. 2b shows an excellent agreement between the exact and approximate solution y 17 (x) to the given problem (99).
where f (x) is selected, such that exact solution is such that the exact solution is y(x) = x 5 cosx.In Table 9, we list the obtained maximum absolute errors with different values of α, β at N = 9, 13, 17.Table 10 compares the obtained absolute errors from the presented method (using N = 17, α = −1/2 and β = −1/2) to those of [13] at various values of x. Figure 4a shows that the obtained absolute errors E 16 (x) and E 17 (x) at α = 1 and β = 0, to demonstrate the convergence of solutions.Additionally, Fig. 4b shows an excellent agreement between the exact and approximate solution y 9 (x) to the given problem (101).
Example 5. Consider the Abel equation of the second kind y(x)y (1)    where the exact solution is y(x) = e −x .In Table 11, we list the obtained maximum absolute errors with different values of α, β at N = 5, 10, 20.Table 12 compares the obtained absolute errors from the presented method (using N = 5, α = 0 and β = 0) to those of [30] at various values of x.
Figure 5a shows Log-errors for various N values, when α = 0 and β = 1, to demonstrate the stability of solutions.Additionally, Fig. 5b shows an excellent agreement between the exact and approximate solution y 6 (x) to the given problem (102).x y(x)y (2) (x) + x 2 y(x) + x y 2 (x) + x 3 y 3 (x) = f (x), 0 ≤ x ≤ 1, where f (x) is selected, such that the exact solution is y(x) = x 1.5 e −x .
In Table 13, we list the obtained maximum absolute errors with different values of α, β at N = 10, 15, 20. Figure 6a shows Log-errors for various N values, when α = 1 and β = 0, to demonstrate the stability of solutions.Additionally, Fig. 6b shows an excellent agreement between the exact and approximate solution y 10 (x) to the given problem (103).
Remark 3. In Example 6, the proposed algorithm is applied to a nonlinear differential equation with a nonsmooth solution.Also, Table 13 shows that the obtained accuracy is less than that obtained in Examples 1-5 where the exact solutions are smooth.Therefore, the effectiveness of the proposed algorithm in the case of problems that have nonsmooth solutions is lower but still accepted.

Conclusions
The proposed algorithm discussed some properties and formulae related to B-polynomials which are utilized with either the boundary conditions (44) or the initial conditions (45) to obtain explicitly the q unknown expansion coefficients (54) or (61), respectively.This algorithm leads to a linear algebraic system of (N − q + 1) equations in (N − q + 1) unknowns which has the matrix form (58) when the numerical solution of a linear differential equation with polynomial coefficients is discussed.While it leads to a nonlinear algebraic system of (N − q + 1) equations in (N − q + 1) unknowns when the numerical solution of a nonlinear differential equation with polynomial coefficients is discussed.In each case, a comparison between the absolute errors obtained by the presented method and other methods is given.Furthermore, in all of the numerical examples presented, agreements were reached for between 10 −16 and 10 −20 .All computations presented were performed on a computer running Mathematica 12 [Intel(R) Core(TM) i9-10850 CPU@ 3.60 GHz, 3600 Mhz, 10 Core(s), 20 Logical Processor (s)].We have shown that the presented B-polynomial method will return a valid solution for a differential equation and is a powerful tool that we may utilize to overcome the difficulties associated with boundary and initial value problems with less computational effort.We have also provided some theoretical results regarding the error resulting from the presented numerical method.These theoretical results have been validated through the presented numerical examples, and the corresponding explanations have been given.

MJOM
Numerical Solutions of High-Order Differential Page 13 of 31 303 b

Figure 4 .
Figure 4. Graph E N and approximate solution for Example 4

Figure 5 .
Figure 5. Graph Log 10 (E N ) and approximate solutions for Example 5

Table 3 .
Computations of y r and y N r for Example 1

Table 6 .
Comparison between the absolute errors of numerical results for Example 2 MJOM

Table 8 .
Comparison between the absolute errors of numerical results for Example 3 Example 4. Consider the following initial value problem to demonstrate that the B-polynomials are powerful to approximate the solution to desired accuracy:

Table 10 .
Comparison between the absolute errors of numerical results for Example 4

Table 12 .
Comparison between the absolute errors of numerical results for Example 5 MJOMExample 6.Consider the boundary value problem

Table 13 .
E N ∞ for N = 10, 15, 20, and various values of (α, β) for Example 6 Figure 6.Graph of Log 10 (E N ) and approximate solutions for Example 6 MJOMNumerical Solutions of High-Order Differential