Abstract

A numerical method for solving differential equations by approximating the solution in the Bernstein polynomial basis is proposed. At first, we demonstrate the relation between the Bernstein and Legendre polynomials. By using this relation, we derive the operational matrices of integration and product of the Bernstein polynomials. Then, we employ them for solving differential equations. The method converts the differential equation to a system of linear algebraic equations. Finally some examples and their numerical solutions are given; comparing the results with the numerical results obtained from the other methods, we show the high accuracy and efficiency of the proposed method.

1. Introduction

In recent years, the Bernstein polynomials (B-polynomials) have attracted the attention of many researchers. These polynomials have been utilized for solving different equations by using various approximate methods. For instance, B-polynomials have been used for solving Fredholm integral equations [1, 2], Volterra integral equations [3], differential equations [47], and integro-differential equations [8]. Singh et al. [6] and Yousefi and Behroozifar [7] have proposed an operational matrix in different ways for solving differential equations. In [6], the B-polynomials have been first orthonormalized by using Gram-schmidt orthonormalization process, and then the operational matrix of integration has been obtained. By the expansion of B-polynomials in terms of Taylor basis, Yousefi and Behroozifar have found the operational matrices of integration and product of B-polynomials. In this paper, firstly, we present operational matrices of integration 𝑃𝑏 and product 𝐶 for the B-polynomials, by the expansion of B-polynomials in terms of Legendre polynomials. Then, we use them for solving differential equation𝑠𝑗=0𝜌𝑗(𝑥)𝑦(𝑗)(𝑥)=𝑔(𝑥),0𝑥1,(1.1)

with the initial conditions𝑦(𝑘)(0)=𝑏𝑘,0𝑘𝑠1,(1.2)

where 𝑔(𝑥) and 𝜌𝑗(𝑥), 𝑗=0,,𝑠 are given functions and 𝑦(𝑥) is the unknown function to be determined. The main characteristic of this technique is that it reduces these equations to those of an easily soluble algebraic equation, thus greatly simplifying the equations. Special attention has been given to the applications of Legendre wavelets method [911], Homotopy perturbation method (HPM) [12], modified decomposition method (MDM) [13], Taylor matrix method [14], and Chebyshev wavelets method [15]. The organization of this paper is as follows: in Section 2, we introduce the B-polynomials and their properties. Section 3 is devoted to the function approximation by using B-polynomials basis. Section 4 introduces the expansion of B-polynomial in terms of Legendre basis, and vice versa. The operational matrices of integration and product will be derived in Section 5. Section 6 is devoted to the solution method of differential equations. In section 7, we present some numerical examples. Numerical solution of each equation based on the exact and approximate solutions are compared. And Section 8 offers our conclusion.

2. B-Polynomials and Their Properties

The B-polynomials of 𝑚th degree are defined on the interval [0,1] as [4]𝐵𝑖,𝑚𝑚𝑖𝑥(𝑥)=𝑖(1𝑥)𝑚𝑖,0𝑖𝑚,(2.1)

where𝑚𝑖=𝑚!𝑖!(𝑚𝑖)!.(2.2)

There are 𝑚+1,𝑚th degree B-polynomials. For mathematical convenience, we usually set 𝐵𝑖,𝑚(𝑥)=0, if 𝑖<0 or 𝑖>𝑚. These polynomials are quite easy to write down: the coefficients can be obtained from Pascal’s triangle. It can easily be shown that each of the B-polynomials is positive and also the sum of all the B-polynomials is unity for all real 𝑥[0,1], that is,𝑚𝑖=0𝐵𝑖,𝑚[](𝑥)=1,𝑥0,1.(2.3)

See [16] for complete details.

3. Function Approximation

B-polynomials defined above form a complete basis [1] over the interval [0,1]. It is easy to show that any given polynomial of degree 𝑚 can be expressed in terms of linear combination of the basis functions. A function 𝑓(𝑥) defined over [0,1] may be expanded as𝑓(𝑥)𝑃𝑚(𝑥)=𝑚𝑖=0𝑐𝑖𝐵𝑖,𝑚(𝑥),𝑚1.(3.1)

Equation (3.1) can be written as𝑃𝑚(𝑥)=𝐶𝑇𝜙(𝑥),(3.2)

where 𝐶 and 𝜙(𝑥) are (𝑚+1)×1 vectors given by𝑐𝐶=0,𝑐1,,𝑐𝑚𝑇,𝐵(3.3)𝜙(𝑥)=0,𝑚(𝑥),𝐵1,𝑚(𝑥),,𝐵𝑚,𝑚(𝑥)𝑇.(3.4)

The use of an orthogonal basis on [0,1] allows us to directly obtain the least-squares coefficients of 𝑃𝑚(𝑥) in that basis, and also ensures permanence of these coefficients with respect to the degree 𝑚 of the approximant, that is, all the coefficients of 𝑃𝑚+1 agree with those of 𝑃𝑚(𝑥), except for that of the newly introduced term. The B-polynomials are not orthogonal. But, these can be expressed in terms of some orthogonal polynomials, such as the Legendre polynomials. The Legendre polynomials constitute an orthogonal basis that is well suited [17, 18] to least-squares approximation.

4. Expansion of B-Polynomials in Terms of Legendre Basis and Vice Versa

To use the Legendre polynomials for our purposes, it is preferable to map this to [0,1]. A set of shifted Legendre polynomials, denoted by {𝐿𝑘(𝑥)} for 𝑘=0,1,, is orthogonal with respect to the weighting function 𝑤(𝑥)=1 over the interval [0,1]. These polynomials satisfy the recurrence relation [19](𝑘+1)𝐿𝑘+1(𝑥)=(2𝑘+1)(2𝑥1)𝐿𝑘(𝑥)𝑘𝐿𝑘1(𝑥),𝑘=1,2,,(4.1)

with𝐿0𝐿(𝑥)=1,1(𝑥)=2𝑥1.(4.2)

The orthogonality of these polynomials is expressed by the relation10𝐿𝑗(𝑥)𝐿𝑘1(𝑥)𝑑𝑥=2𝑘+1,𝑗=𝑘,0,𝑗𝑘,𝑗,𝑘=0,1,2,.(4.3)

When the approximant (3.1) is expressed in the Legendre form𝑃𝑚(𝑥)=𝑚𝑗=0𝑙𝑗𝐿𝑗(𝑥),(4.4)

by using (4.3), we can obtain the Legendre coefficients as𝑙𝑗=(2𝑗+1)10𝐿𝑗(𝑥)𝑓(𝑥)𝑑𝑥,𝑗=0,,𝑚.(4.5)

Lemma 4.1. The Legendre polynomial 𝐿𝑘(𝑥) can be expressed in the kth degree Bernstein basis 𝐵0,𝑘(𝑥),𝐵1,𝑘(𝑥),,𝐵𝑘,𝑘(𝑥) as [20] 𝐿𝑘(𝑥)=𝑘𝑖=0(1)𝑘+𝑖𝑘𝑖𝐵𝑖,𝑘(𝑥).(4.6)

Now consider a polynomial 𝑃𝑚(𝑥) of degree 𝑚, expressed in the 𝑚th degree Bernstein and Legendre bases on 𝑥[0,1]𝑃𝑚(𝑥)=𝑚𝑗=0𝑐𝑗𝐵𝑗,𝑚(𝑥)=𝑚𝑘=0𝑙𝑘𝐿𝑘(𝑥).(4.7)

We write the transformation of the Legendre polynomials on [0,1] into the 𝑚th degree Bernstein basis functions as [21]𝐵𝑘,𝑚(𝑥)=𝑚𝑖=0𝑤𝑘,𝑖𝐿𝑖(𝑥),𝑘=0,,𝑚.(4.8)

The elements 𝑤𝑘,𝑖, 𝑘,𝑖=0,1,,𝑚, form an (𝑚+1)×(𝑚+1) basis conversion matrix 𝑊. To compute them, we multiply (4.8) by 𝐿𝑗(𝑥), integrate over 𝑥[0,1], and use (4.3) to obtain𝑤𝑘,𝑗=(2𝑗+1)10𝐵𝑘,𝑚(𝑥)𝐿𝑗(𝑥)𝑑𝑥.(4.9)

We now replace (4.6) into (4.9) and obtain𝑤𝑘,𝑗=(2𝑗+1)𝑗𝑖=0(1)𝑗+𝑖𝑗𝑖10𝐵𝑘,𝑚(𝑥)𝐵𝑖,𝑗(𝑥)𝑑𝑥.(4.10)

The integrals of the products of Bernstein basis functions can be found using10(1𝑥)𝑟𝑥𝑖1𝑑𝑥=(𝑟+𝑖+1)𝑖𝑟+𝑖,𝑖,𝑟𝑁{0},(4.11)

as follows:10𝐵𝑘,𝑚(𝑥)𝐵𝑖,𝑗𝑚𝑘𝑗𝑖(𝑥)𝑑𝑥=10𝑥𝑘+𝑖(1𝑥)𝑚+𝑗𝑘𝑖(𝑑𝑥=𝑚𝑘)𝑗𝑖(𝑚+𝑗+1)𝑚+𝑗𝑘+𝑖.(4.12)

Therefore, we have the elements of 𝑊 as𝑤𝑘,𝑗=(2𝑗+1)𝑚𝑘𝑚+𝑗+1𝑗𝑖=0(1)𝑗+𝑖𝑗𝑖𝑗𝑖𝑚+𝑗𝑘+𝑖,𝑘,𝑗=0,,𝑚.(4.13)

Now, we write the transformation of the B-polynomials on [0,1] into 𝑚th degree Legendre basis functions as [21]𝐿𝑘(𝑥)=𝑚𝑗=0Λ𝑘,𝑗𝐵𝑗,𝑚(𝑥),𝑘=0,,𝑚.(4.14)

The elements Λ𝑘,𝑗 form an (𝑚+1)×(𝑚+1) basis conversion matrix Λ. Replacing (4.14) into (4.7) and rearranging the order of summation, we obtain𝑐𝑗=𝑚𝑘=0𝑙𝑘Λ𝑘,𝑗,𝑗=0,,𝑚.(4.15)

Since we can express each 𝑘th degree Bernstein basis function in the 𝑚th degree Bernstein basis as [21]𝐵𝑖,𝑘(𝑥)=𝑚𝑘+𝑖𝑗=𝑖𝑘𝑖𝑚𝑘𝑗𝑖𝑚𝑗𝐵𝑗,𝑚(𝑥),𝑖=0,,𝑘,(4.16)

replacing (4.16) into (4.6) and rearranging the order of summation, we find that the basis transformation (4.14) is defined by the elementsΛ𝑘,𝑗=1𝑚𝑗min{𝑗,𝑘}𝑖=𝑟(1)𝑘+𝑖𝑘𝑖𝑘𝑖𝑚𝑘𝑗𝑖,𝑟=max{0,𝑗+𝑘𝑚},(4.17)

of the matrix Λ for 𝑘,𝑗=0,,𝑚. If we denote the Legendre basis vector as𝐿𝐿(𝑥)=0(𝑥),𝐿1(𝑥),,𝐿𝑚(𝑥)𝑇,(4.18)

using (3.4), (4.8), (4.14), and (4.18), we have𝜙(𝑥)=𝑊𝐿(𝑥),(4.19)𝐿(𝑥)=Λ𝜙(𝑥).(4.20)

5. Operational Matrices of Integration and Product of B-Polynomials

5.1. B-Polynomials Operational Matrix of Integration

Let 𝑃𝑏 be an (𝑚+1)×(𝑚+1) operational matrix of integration, then𝑥0𝜙(𝑡)𝑑𝑡𝑃𝑏𝜙(𝑥),0𝑥1.(5.1)

By using (4.19), we have𝑥0𝜙(𝑡)𝑑𝑡=𝑊𝑥0(𝑡)𝑑𝑡=𝑊𝑃𝐿(𝑥),(5.2)

where the (𝑚+1)×(𝑚+1) matrix 𝑃 is the operational matrix of integration of the shifted Legendre polynomials on the interval [0,1] and can be obtained as [22]1𝑃=2110000013013000001501500000001012𝑚12𝑚100000102𝑚+1,(5.3)

and, therefore, by using (4.20)–(5.2), we have the operational matrix of integration as𝑃𝑏=𝑊𝑃Λ.(5.4)

5.2. B-Polynomials Operational Matrix of Product

In this subsection, we present a general formula for finding the operational matrix of product of 𝑚th degree B-polynomials. Suppose that 𝐶 is an arbitrary (𝑚+1)×1 vector, then 𝐶 is an (𝑚+1)×(𝑚+1) operational matrix of product whenever𝐶𝑇𝜙(𝑥)𝜙𝑇(𝑥)𝜙𝑇(𝑥)𝐶.(5.5)

Using (4.19) and since 𝐶𝑇𝜙(𝑥)=𝑚𝑖=0𝑐𝑖𝐵𝑖,𝑚, we have𝐶𝑇𝜙(𝑥)𝜙𝑇=𝐶(𝑥)𝑇𝐿𝜙(𝑥)𝑇(𝑥)𝑊𝑇=𝐿0𝐶(𝑥)𝑇𝜙(𝑥),𝐿1𝐶(𝑥)𝑇𝜙(𝑥),,𝐿𝑚𝐶(𝑥)𝑇𝑊𝜙(𝑥)𝑇=𝑚𝑖=0𝑐𝑖𝐿0(𝑥)𝐵𝑖,𝑚,(𝑥)𝑚𝑖=0𝑐𝑖𝐿1(𝑥)𝐵𝑖,𝑚(𝑥),,𝑚𝑖=0𝑐𝑖𝐿𝑚(𝑥)𝐵𝑖,𝑚𝑊(𝑥)𝑇.(5.6)

Now, we approximate all functions 𝐿𝑘(𝑥)𝐵𝑖,𝑚(𝑥) in terms of {𝐵𝑖,𝑚(𝑥)}𝑚𝑖=0 for 𝑖,𝑘=0,1,,𝑚. Let𝜂𝑘,𝑖=𝜂0𝑘,𝑖𝜂1𝑘,𝑖𝜂𝑚𝑘,𝑖,(5.7)

by (3.2), we have𝐿𝑘(𝑥)𝐵𝑖,𝑚(𝑥)𝜂𝑇𝑘,𝑖𝜙(𝑥),𝑖,𝑘=0,1,,𝑚.(5.8)

Using (4.15), we can obtain the elements of vector 𝜂𝑘,𝑖, for 𝑖,𝑘=0,1,,𝑚. Therefore,𝑚𝑖=0𝑐𝑖𝐿𝑘(𝑥)𝐵𝑖,𝑚(𝑥)𝑚𝑖=0𝑐𝑖𝑚𝑗=0𝜂𝑗𝑘,𝑖𝐵𝑗,𝑚=(𝑥)𝑚𝑗=0𝐵𝑗,𝑚(𝑥)𝑚𝑖=0𝑐𝑖𝜂𝑗𝑘,𝑖=𝜙𝑇(𝑥)𝑚𝑖=0𝑐𝑖𝜂0𝑚𝑘,𝑖𝑖=0𝑐𝑖𝜂1𝑘,𝑖𝑚𝑖=0𝑐𝑖𝜂𝑚𝑘,𝑖=𝜙𝑇𝜂(𝑥)𝑘,0,𝜂𝑘,1,,𝜂𝑘,𝑚𝐶=𝜙𝑇𝐶(𝑥)𝑘,(5.9)

where𝐶𝑘=𝜂𝑘,0,𝜂𝑘,1,,𝜂𝑘,𝑚𝐶,𝑘=0,1,,𝑚.(5.10)

If we define a (𝑚+1)×(𝑚+1) matrix 𝐶𝐶=[0,𝐶1𝐶,,𝑚], then by using (5.6) and (5.9), we have𝐶𝑇𝜙(𝑥)𝜙𝑇(𝑥)𝜙𝑇𝐶(𝑥)0,𝐶1𝐶,,𝑚𝑊𝑇=𝜙𝑇(𝑥)𝐶𝑊𝑇,(5.11)

and, therefore, we have the operational matrix of product as𝐶=𝐶𝑊𝑇.(5.12)

6. Solution of the Linear Differential Equation

Consider the linear differential equation (1.1) with the initial conditions (1.2). If we approximate 𝑔(𝑥), 𝜌𝑗(𝑥), 𝑗=0,,𝑠 and 𝑦(𝑠)(𝑥) as follows:𝑔(𝑥)=𝐺𝑇𝜌𝜙(𝑥),𝑗(𝑥)=𝑃𝑇𝑗𝑦𝜙(𝑥),𝑗=0,,𝑠,(6.1)(𝑠)(𝑥)=𝐶𝑇𝜙(𝑥),(6.2)

where 𝐺, 𝑃𝑗, 𝑗=0,,𝑠, and 𝐶 are the coefficients which are defined similarly to (3.3). With 𝑠-times integrating from (6.2) with respect to 𝑥 between 𝑥=0 to 𝑥=𝑥, using (5.1) and the initial conditions (1.2), we will have𝑦(𝑠1)(𝑥)=𝑏𝑠1+𝐶𝑇𝑃𝑏𝑦𝜙(𝑥),(𝑠2)(𝑥)=𝑏𝑠2+𝑏𝑠1𝑥+𝐶𝑇𝑃2𝑏𝜙(𝑥),𝑦(𝑥)=𝑏1+𝑏2𝑏𝑥+3𝑥2!2𝑏++𝑠1𝑥(𝑠2)!(𝑠2)+𝐶𝑇𝑃𝑏𝑠1𝜙(𝑥),𝑦(𝑥)=𝑏0+𝑏1𝑏𝑥+2𝑥2!2𝑏++𝑠1(𝑥𝑠1)!(𝑠1)+𝐶𝑇𝑃𝑠𝑏𝜙(𝑥).(6.3)

Let𝑥𝑖=𝑑𝑇𝑖𝑏𝜙(𝑥),𝑖=1,2,,𝑠1,𝑠𝑖=𝑏𝑠𝑖𝐸𝑇𝜙(𝑥),𝑖=1,2,,𝑠,(6.4)

where1=𝐸𝑇𝜙(𝑥).(6.5)

Substituting (6.4) into (6.3), we have𝑦(𝑠)(𝑥)=𝐶𝑇𝜙(𝑥)=𝑄𝑇𝑠𝑦𝜙(𝑥),(𝑠1)𝑏(𝑥)=𝑠1𝐸𝑇+𝐶𝑇𝑃𝑏𝜙(𝑥)=𝑄𝑇𝑠1𝑦𝜙(𝑥),(𝑠2)(𝑏𝑥)=𝑠2𝐸𝑇+𝑏𝑠1𝑑𝑇1+𝐶𝑇𝑃2𝑏𝜙(𝑥)=𝑄𝑇𝑠2𝑏𝜙(𝑥),𝑦(𝑥)=1𝐸𝑇+𝑏2𝑑𝑇1+𝑏3𝑑2!𝑇2𝑏++𝑠1𝑑(𝑠2)!𝑇𝑠2+𝐶𝑇𝑃𝑏𝑠1𝜙(𝑥)=𝑄𝑇1𝑏𝜙(𝑥),(6.6)𝑦(𝑥)=0𝐸𝑇+𝑏1𝑑𝑇1+𝑏2𝑑2!𝑇2𝑏++𝑠1(𝑑𝑠1)!𝑇𝑠1+𝐶𝑇𝑃𝑠𝑏𝜙(𝑥)=𝑄𝑇0𝜙(𝑥).(6.7)

Replacing (6.6) and (6.7) into (1.1), we obtain𝑠𝑗=0𝑃𝑇𝑗𝜙(𝑥)𝜙𝑇(𝑥)𝑄𝑗=𝐺𝑇𝜙(𝑥).(6.8)

Using (5.5), we have𝑠𝑗=0𝜙𝑇𝑃(𝑥)𝑗𝑄𝑗=𝜙𝑇(𝑥)𝐺.(6.9)

Therefore, we get𝑠𝑗=0𝑃𝑗𝑄𝑗=𝐺.(6.10)

The unknown vector 𝐶 can be obtained by solving (6.10). Once 𝐶 is known, 𝑦(𝑥) can be calculated from (6.7).

7. Illustrative Examples

Example 7.1. Consider the eighth-order linear differential equation given in [14] by 𝑦(𝑣𝑖𝑖𝑖)(𝑥)𝑦(𝑥)=8𝑒𝑥,0𝑥1,(7.1)with the initial conditions 𝑦(0)=1,𝑦(0)=0,𝑦𝑦(0)=1,(0)=2,𝑦(𝑖𝑣)(0)=3,𝑦(𝑣)(𝑦0)=4,(𝑣𝑖)(0)=5,𝑦(𝑣𝑖𝑖)(0)=6.(7.2)The exact solution for this example is 𝑦(𝑥)=(1𝑥)𝑒𝑥. Using the method described in Section 6, we assume that 𝑦(𝑣𝑖𝑖𝑖)(𝑥) is approximated by 𝑦(𝑣𝑖𝑖𝑖)(𝑥)=𝐶𝑇𝜙(𝑥).(7.3) By using (5.1) and the initial conditions (7.2), we have 𝑦(𝑥)=𝑄0𝜙(𝑥),(7.4) where 𝑄0=𝐸𝑇1𝑑2!𝑇22𝑑3!𝑇33𝑑4!𝑇44𝑑5!𝑇55𝑑6!𝑇66𝑑7!𝑇7+𝐶𝑇𝑃8𝑏=𝐴𝑇+𝐶𝑇𝑃8𝑏.(7.5) We can express function 8𝑒𝑥 as 8𝑒𝑥=𝐺𝑇𝜙(𝑥).(7.6) Substituting (7.3)–(7.6) into (7.1), we obtain 𝐶𝑇𝐴𝜙(𝑥)𝑇+𝐷𝑇𝑃8𝑏𝜙(𝑥)=𝐺𝑇𝜙(𝑥).(7.7) Therefore, we get 𝑃𝐶=𝐼8𝑏𝑇1(𝐴+𝐺),(7.8) where 𝐼 is the (𝑚+1)×(𝑚+1) identity matrix and Equation (7.8) is a set of algebraic equations which can be solved for 𝐶. Now, we apply the method presented in this paper with 𝑚=8 to solve (7.1) with the initial conditions (7.2). In Table 1, the numerical results obtained by the present method are compared with the results of the HPM [12] and MDM [13] and method in [14]. As we see from this table, it is clear that the result obtained by the present method is very superior to that by HPM, MDM methods, and method in [14]. The absolute difference between exact and approximate solutions is plotted in Figure 1. It is observed in this figure that the accuracy is of the order of 1010.

Example 7.2. Consider the Lane-Emden equation given in [10] by 𝑦2(𝑥)+𝑥𝑦(𝑥)=22𝑥2+3𝑦(𝑥),0𝑥1,(7.9) with the initial conditions 𝑦(0)=1,𝑦(0)=0.(7.10) The exact solution of this example is 𝑒𝑥2. We solve (7.9) with the initial conditions (7.10) by using the method in Section 6 with 𝑚=11. The comparison among the present method, Legendre wavelets solution [10], and analytic solution for 𝑚=11 is shown in Table 2. As we see from this Table, it is clear that the result obtained by the present method is very superior to that by Legendre wavelets method. It is noted that the mean square error for this example, obtained in Legendre wavelets method, is 2.07×105; but in the present method, the mean square error is 2.3174×1018. We display a plot of the approximate and exact solution of this example for 𝑚=11 in Figure 2.

Example 7.3. Consider the Bessel differential equation of order zero given in [15] by 𝑥𝑦+𝑦+𝑥𝑦=0,(7.11) with the initial conditions 𝑦(0)=1,𝑦(0)=0.(7.12) The exact solution of this example is 𝒥0(𝑥)=𝑞=0(1)𝑞(𝑞!)2𝑥22𝑞.(7.13) Now, we solve (7.11) with the initial conditions (7.12) by using the method in Section 6 with 𝑚=5. Table 3 shows the absolute difference between exact and approximate solutions of the present method and the methods in [11, 15] for equality basis functions. The results of [11] have been given in [15]. As we see from this table, the maximum error for this example, for the methods in [11, 15], is 104; but in the present method, the maximum error is 107. We display a plot of the approximate and exact solution of this example for 𝑚=5 in Figure 3.

8. Conclusion

In this article, at first, we demonstrate the relation between the Bernstein and Legendre polynomials. By using this relation, we derived the operational matrix of integration and product of B-polynomials. They are applied to solve ordinary differential equations. The present method reduces an ordinary differential equations into a set of algebraic equations. We applied the presented method on three test problems and compared the results with their exact solutions and the other methods, revealing that the present method is very effective and convenient.

Acknowledgment

The work was supported by Alzahra university.