# A Collocation Method Based on the Bernoulli Operational Matrix for Solving Nonlinear BVPs Which Arise from the Problems in Calculus of Variation

Authors: Emran Tohidi; Adem Kılıçman
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A new collocation method is developed for solving BVPs which arise from the problems in calculus of variation. These BVPs result from the Euler-Lagrange equations, which are the necessary conditions of the extremums of problems in calculus of variation. The proposed method is based upon the Bernoulli polynomials approximation together with their operational matrix of differentiation. After imposing the collocation nodes to the main BVPs, we reduce the variational problems to the solution of algebraic equations. It should be noted that the robustness of operational matrices of differentiation with respect to the integration ones is shown through illustrative examples. Complete comparisons with other methods and superior results confirm the validity and applicability of the
presented method.


---

## Body

## 1. Introduction

In a large number of applied sciences problems such as analysis, mechanics, and geometry, it is necessary to determine the maximal and minimal of a certain function. Because of the important role of this subject in sciences and engineering, considerable attention has been provided on this kind of problems. Such problems are called variational problems. The various applications of variational problems such as industrial and biological are introduced in [1]. Since a huge size of such equations cannot be solved explicitly, it is often necessary to resort to the approximation and numerical techniques.In the recent years, the studies on variational problems were developed very rapidly and intensively. For instance, one can point out to the methods that based upon operational matrices of integration of a huge size of polynomials and functions. In the last four decades, numerical methods which are based on the operational matrices of integration (especially for orthogonal polynomials and functions) have received considerable attention for dealing with variational problems. The key idea of these methods is based on the integral expression(1)∫0x‍Φ(τ)dτ≈Φ(x)P,
where Φ(t)=[Φ1(x),Φ2(x),…,ΦN(x)] is an arbitrary basis vector and P is the (N+1)×(N+1) constant matrix called the operational matrix of integration. The matrix P has already been determined for many types of orthogonal (or nonorthogonal) bases such as the Walsh functions [2], the Laguerre polynomials [3], the Chebyshev polynomials [4], the Legendre polynomials [5], and the Fourier series [6]. However, methods that are based on the high-order Gauss quadrature rules [7, 8] could be applied for variational problems, but the need to more CPU times and ill-conditioning of the associated algebraic systems are some disadvantages of these approaches.On the other hand, since the beginning of 1994, the Laguerre, Chebyshev, Taylor, Legendre, Hermite, Fourier, and Bessel (matrix and collocation) methods have been used in the works [9–16] to solve the high-order linear and the nonlinear differential (including hyperbolic partial differential equations) Fredholm-Volterra integrodifferential difference delay equations and their systems. Also, the Bernoulli matrix method has been used to find the approximate solutions of differential and integrodifferential equations [17–20]. The main characteristic of these approaches is based on the operational matrices of differentiation instead of integration. The best advantage of these techniques with respect to the integration methods is that, in the fundamental matrix relations, there is not any approximation symbol; meanwhile, in the integration forms such as (1), the approximation symbol could be seen obviously. In other words,
(2)Φ′(x)=Φ(x)B,
where B is the operational matrix of differentiation for any selected basis such as the above-mentioned polynomials, functions, and truncated series. The readers can see that there is no approximation symbol in (2); meanwhile, this can be seen in (1) by using operational matrices of integration. For justifying this expression, one can refer to this subject that after differentiating an Nth degree polynomial, we usually reach to a polynomial which has less than Nth degree. However, in the integration processes, the degree of polynomials would be increased.In this paper, in the light of the above-mentioned methods and by means of the matrix relations between the Bernoulli polynomials and their derivatives, we develop a new method called the Bernoulli collocation method (BCM) for finding the extremum of the functional(3)J[y1(x),y2(x),…,yn(x)]=∫01‍G(yn′(x)x,y1(x),y2(x),…,yn(x),y1′(x),y2′(x),…,yn′(x))dx.
To find the extreme value of J, the boundary points of the admissible curves are known in the following form:
(4)yi(0)=γi,i=1,2,…,n,yi(1)=δi,i=1,2,…,n.
The necessary condition to extremize J[y1(x),y2(x),…,yn(x)] is that it should satisfy the Euler-Lagrange equations
(5)∂G∂yi-ddx(∂G∂yi′)=0,i=1,2,…,n,
with boundary conditions given in (4). The system of boundary value problems (5) does not always have a solution, and if the solution exists, it may not be unique. Note that in many variational problems, the existence of a solution is obvious from the physical or geometrical meaning of the problem, and if the solution of the Euler equation satisfies the boundary conditions, it is unique. Also, this unique extremal will be the solution of the given variational problem. Thus, another approach for solving the variational problem (3) is finding the solution of the system of ordinary differential equations (5) which satisfies the boundary conditions (4) which were called systems of boundary value problems (BVPs). The simplest form of the variational problem (3) is
(6)J[y(x)]=∫01‍G(x,y(x),y′(x))dx,
with the given boundary conditions
(7)y(0)=γ,y(1)=δ.
Here, the necessary condition for the extremum of the functional (6) is to satisfy the following second-order differential equation:
(8)∂G∂y-ddx(∂G∂y′)=0,
with boundary conditions given in (7).We again emphasize that our aim is to solve BVPs such as (5) and (8) by using our method which has superior results with respect to several methods in the literature. It should be noted that the handeling of BVPs needs more accuracy with respect to the initial value problems (IVPs). For instance, see [21–23] and the references therein. To the best of our knowledge, this is the first work concerning the Bernoulli collocation method for solving nonlinear BVPs. This partially motivated our interest in such method.The rest of this paper is organized as follows. In Section2, we review several properties of the Bernoulli polynomials. Section 3 is devoted to the basic idea of this paper (i.e., the Bernoulli collocation method). Error analysis and accuracy of the approximate solution is given in Section 4. Several illustrative examples are provided in Section 5 for confirming the effectiveness of the presented method. Section 6 contains some conclusions and notations about the future works.


## 2. The Bernoulli Polynomials and Their Operational Matrix

The Bernoulli polynomials play an important role in different areas of mathematics, including number theory and the theory of finite differences. The classical Bernoulli polynomialsBn(x) are usually defined by means of the exponential generating functions (see [18]):
(9)wexwew-1=∑k=0∞‍Bk(x)wkk!,(|w|<2π).
The following familiar expansion (see [17]):
(10)∑k=0n‍(n+1k)Bk(x)=(n+1)xn,
is the most primary property of the Bernoulli polynomials. The first few Bernoulli polynomials are
(11)B0(x)=1,B1(x)=x-12,B2(x)=x2-x+16,B3(x)=x3-32x2+12x,B4(x)=x4-2x3+x2-130.
The Bernoulli polynomials satisfy the well-known relations (see [18]):
(12)dBn(x)dx=nBn-1(x),(n≥1),∫01‍Bn(x)dx=0,(n≥1).
If we introduce the Bernoulli vector B(x) in the form B(x)=[B0(x),B1(x),…,BN(x)], then the derivative of the B(x), with the aid of the first property of (12), can be expressed in the matrix form by
(13)[B0′(x)B1′(x)B2′(x)⋮BN-1′(x)BN′(x)]︸B′(x)T=[000⋯000100⋯000020⋯000⋮⋮⋮⋱⋮⋮⋮000⋯N-100000⋯0N0]︸M[B0(x)B1(x)B2(x)⋮BN-1(x)BN(x)]︸BT(x),
where M is the (N+1)×(N+1) operational matrix of differentiation.Accordingly, thekth derivative of B(x) can be given by
(14)B′(x)T=MB(x)T⇒B(1)(x)=B(x)MT,B(2)(x)=B(1)(x)MT=B(x)(MT)2,B(3)(x)=B(1)(x)(MT)2=B(x)(MT)3,⋮B(k)(x)=B(x)(MT)k,
where M is defined in (13).


## 3. Basic Idea

Now, we consider the general form of the variational problem (3). Finding the solution of the problem (3) needs to solve the corresponding BVP (5) with boundary conditions (4). We assume that (y1(x),y2(x),…,yn(x)) is the exact solution of the BVP (5). Our aim is to approximate (y1(x),y2(x),…,yn(x)) over the interval [0,1] by the following linear combinations of the Bernoulli polynomials:
(15)y1(x)≈y1,N(x)=∑i=0N‍y1,iBi(x)=B(x)Y1,y2(x)≈y2,N(x)=∑i=0N‍y2,iBi(x)=B(x)Y2,⋮yn(x)≈yn,N(x)=∑i=0N‍yn,iBi(x)=B(x)Yn,
where B(x)=[B0(x),B1(x),…,BN(x)], Yk=[yk,0,yk,1,…,yk,N]T in which yk,i for i=0,1,…,N and k=1,2,…,n are the unknown coefficients and Bi(x) are the Bernoulli polynomials which are defined in the previous section. For convenience, consider the second-order BVP (5) as follows:
(16)F(y1(x),y2(x),…,yn(x),y1′(x),y2′(x),…,yn′(x),y1′′(x),y2′′(x),…,yn′′(x))=0.
By using (14), for k=1 and k=2, we have
(17)y1′(x)≈y1,N′(x)=B(x)MTY1,y1′′(x)≈y1,N′′(x)=B(x)(MT)2Y1,y2′(x)≈y2,N′(x)=B(x)MTY2,y2′′(x)≈y2,N′′(x)=B(x)(MT)2Y2,⋮yn′(x)≈yn,N′(x)=B(x)MTYn,yn′′(x)≈yn,N′′(x)=B(x)(MT)2Yn.
We consider the matrix vector forms of yk(x), yk′(x), and yk′′(x) for k=1,2,…,n which have been shown in (15) and (17) and then replace in (16) as follows:
(18)F(B(x)(MT)2YnB(x)Y1,…,B(x)Yn,B(x)MTY1,…,B(x)MTYn,B(x)(MT)2Y1,…,B(x)(MT)2Yn)=0.
Because of stability properties of the Gaussian points [7], we collocate the above BVP at the nodes xj for j=1,2,…,n(N-1) as the roots of the n(N-1)th degree shifted Legendre polynomial Pn(N-1)*(x)=Pn(N-1)(2x-1) as follows:
(19)F((MT)2YnB(xj)Y1,…,B(xj)Yn,B(xj)MTY1,…,B(xj)MTYn,B(xj)(MT)2Y1,…,B(xj)(MT)2Yn)=0,
where 1≤j≤n(N-1). Note that according to the properties of Gaussian points we have 0<xj<1 for all values of j. The above system consists of n(N-1) equations with n(N+1) unknowns. Now, consider the 2n equations from boundary conditions as follows
(20)y1(0)≈y1,N(0)=B(0)Y1=γ1,y1(1)≈y1,N(1)=B(1)Y1=δ1,y2(0)≈y2,N(0)=B(0)Y2=γ2,y2(1)≈y2,N(1)=B(1)Y2=δ2,⋮yn(0)≈yn,N(0)=B(0)Yn=γn,yn(1)≈yn,N(1)=B(1)Yn=δn.
The above equations together with (19) form a nonlinear algebraic system with n(N+1) equations and n(N+1) unknowns yk,i for i=0,1,…,N and k=1,2,…,n. After solving this algebraic system, we obtain the approximated solutions of (3).


## 4. Error Analysis and Accuracy of the Solution

This section is devoted to provide an error bound for the approximated solution which was presented in the previous section. Before presenting the main theorem of this section, we need to recall some useful corollaries and theorems. Then, we transform the basic equation (8) (or (5)) together with the boundary conditions (7) (or (4)) to a nonlinear Fredholm-Volterra integral equation (or system of nonlinear Fredholm-Volterra integral equations). Therefore, the main theorem could be stated which guarantees the convergence of the truncated Bernoulli series to the exact solution under several mild conditions.Now, suppose thatH=L2[0,1] and {B0(x),B1(x),…,BN(x)}⊂H is the set of the Bernoulli polynomials and
(21)Y=span{B0(x),B1(x),…,BN(x)},
and g is an arbitrary element in H. Since Y is a finite dimensional vector space, g has the unique best approximation belongs to Y such as g^∈Y, that is,
(22)∀y∈Y,∥g-g^∥≤∥g-y∥.
Since g^∈Y, there exist the unique coefficients g0,g1,…,gN such that
(23)g≈g^=∑n=0N‍gnBn(x)=B(x)GT,B(x)=[B0(x),B1(x),…,BN(x)],G=[g0,g1,…,gN].Corollary 1.
Assume thatg∈H=L2[0,1] is an arbitrary function and also is approximated by the truncated Bernoulli serie ∑n=0∞‍gnBn(x), then the coefficients gn for all n=0,1,…,∞ can be calculated from the following relation:
(24)gn=1n!∫01‍g(n)(x)dx.Proof.
(See [18]).Corollary 2.
Assume that one approximates the functiong on the interval [0,1] by the Bernoulli polynomials as discussed in Corollary 1. Then, the coefficients gn decay as follows:
(25)gn≤Gnn!,
where Gn denotes the maximum of g(n) in the interval [0,1].Proof.
Since it is trivial, we omit the proof.The above corollary implies that the Bernoulli coefficients are decayed rapidly as the increasing ofn.Consider Corollary1 again. We will provide the error of the associated approximation.Theorem 3 (see [17]).
Suppose thatg(x) is an enough smooth function in the interval [0,1] and is approximated by the Bernoulli polynomials as done in Corollary 1. With more details, assume that PN[g](x) is the approximate polynomial of g(x) in terms of the Bernoulli polynomials and RN[g](x) is the remainder term. Then, the associated formulas are stated as follows:
(26)g(x)=PN[g](x)+RN[g](x),x∈[0,1],PN[g](x)=∫01‍g(x)dx+∑j=1N‍Bj(x)j!(g(j-1)(1)-g(j-1)(0)),RN[g](x)=-1N!∫01‍BN*(x-t)g(N)(t)dt,
where BN*(x)=BN(x-[x]) and [x] denotes the largest integer not greater than x.Proof.
See [17].Trivially, the algebraic degree of exactness of the operatorPN[·] is N.Theorem 4.
Suppose thatg(x)∈C∞[0,1] and PN[g](x) is the approximate polynomial using the Bernoulli polynomials. Then, the error bound would be obtained as follows:
(27)∥error(g(x))∥∞≤1N!BNGN,
where BN and GN denote the maximum value of BN(x) and g(N)(x) in the interval [0,1], respectively.Proof.
By consideringRN[g](x)=-(1/N!)∫01‍BN*(x-t)g(N)(t)dt, the proof is clear.Corollary 5.
Assume thatK(x,t)∈H×H=L2[0,1]×L2[0,1] is an arbitrary function and also is approximated by the two variable truncated Bernoulli series ∑m=0N‍∑n=0N‍km,nBm(x)Bn(t), then the coefficients km,n for all m,n=0,1,…,N can be calculated from the following relation:
(28)km,n=1m!n!∬01‍∂m+nK(x,t)∂xm∂tndxdt,m,n=0,1,…,N.Proof.
See [18].Theorem 6.
Suppose thatK(x,t) is a smooth enough function and PN[K](x,t) is the approximate polynomial using the Bernoulli method. Then, the error bound would be obtained as follows:
(29)∥error(K(x,t))∥∞≤1N!BN2KN,N,
where BN and Kj,N,N denote the maximum value of BN(x) and Kj(N,N)(x,t) in the interval [0,1], respectively.Proof.
See [20].For the clarity of presentation, we only consider the equation (8) in the following simple form
(30)y′′(x)=f(x,y(x),y′(x)),y(0)=y0,y(1)=y1,
where f is a continuously differentiable function with respect to its arguments. Our aim is to transform the above BVP into a nonlinear Fredholm-Volterra equation. Therefore, we assume that y′′(x)=u(x). By integrating both sides of the mentioned relation in the interval [0,x], we have
(31)y′(x)=y′(0)+∫0x‍u(t)dt.
Again, one can integrate the above integral equation in the interval [0,x] as follows:
(32)y(x)=y0+xy′(0)+∫0x‍∫0t‍u(τ)dτdt=y0+xy′(0)+∫0x‍(x-t)u(t)dt.
Taking x=1 in (32) yields
(33)y′(0)=y1-y0-∫01‍(1-t)u(t)dt.
According to (31) and (32) we should have
(34)y(x)=y0+x(y1-y0-∫01‍(1-t)u(t)dt)+∫0x‍(x-t)u(t)dt,y′(x)=y1-y0-∫01‍(1-t)u(t)dt+∫0x‍u(t)dt.
Therefore, the nonlinear BVP (30) transformed into the following nonlinear Fredholm-Volterra integral equation:
(35)u(x)=g(x,∫01‍K1(x,t)u(t)dt,∫0x‍K2(x,t)u(t)dt),
where g, K1, and K2 are continuously differentiable functions with respect to their arguments.Theorem 7.
Assume that (35) has a uniqe solution u(x). Also, suppose that the kernels K1(x,t) and K2(x,t) are approximated by the Bernoulli truncated series as shown in Corollary 5, and, hence, by these approximations the equation (35) has the numerical solution uN(x) (in terms of the Bernoulli polynomials). Then, limN→∞uN(x)=u(x) provided by the conditions ||u(x)||∞≤U, ||K1,N(x,t)||∞≤K^1, ||K2,N(x,t)||∞≤K^2 and Lg(K^1+K^2)≪1.Proof.
According to the assumption, we should have(36)uN(x)=g(x,∫01‍K1,N(x,t)u(t)dt,∫0x‍K2,N(x,t)u(t)dt).
Subtracting (36) from (35) yields
(37)u(x)-uN(x)=g(x,∫01‍K1(x,t)u(t)dt,∫0x‍K2(x,t)u(t)dt)-g(x,∫01‍K1,N(x,t)uN(t)dt,∫0x‍K2,N(x,t)uN(t)dt).
Since g is a continuously differentiable function, it is jointly Lipschitz with respect to its second and third arguments, in other words,
(38)∥u(x)-uN(x)∥∞≤Lg(∥∫01‍K1(x,t)u(t)dt-∫01‍K1,N(x,t)uN(t)dt∥∞+∥∫0x‍K2(x,t)u(t)dt-∫0x‍K2,N(x,t)uN(t)dt∥∞)≤Lg(∥K1(x,t)-K1,N(x,t)∥∞∥u(x)∥∞+∥K1,N(x,t)∥∞∥u(x)-uN(x)∥∞+∥K2(x,t)-K2,N(x,t)∥∞∥u(x)∥∞+∥K2,N(x,t)∥∞∥u(x)-uN(x)∥∞)=Lg(E(K1)U+K^1∥u(x)-uN(x)∥∞+E(K1)U+K^2∥u(x)-uN(x)∥∞),
where ||u(x)||∞≤U, ||K1,N(x,t)||∞≤K^1, ||K2,N(x,t)||∞≤K^2, E(K1)=||K1(x,t)-K1,N(x,t)||∞ and E(K2)=||K2(x,t)-K2,N(x,t)||∞. By factorizing the coefficients of ||u(x)-uN(x)||∞, we have
(39)(1-Lg(K^1+K^2))∥u(x)-uN(x)∥∞≤ULg(E(K1)+E(K2)).
More precisely,
(40)∥u(x)-uN(x)∥∞≤ULg(E(K1)+E(K2))1-Lg(K^1+K^2).
Since K1 and K2 are enough differentiable, we can use Theorem 6, and hence limN→∞E(K1)=limN→∞E(K2)=0. On the other hand, because of Lg(K^1+K^2)≪1, one can deduce that limN→∞||u(x)-uN(x)||∞=0, and this completes the proof.


## 5. Numerical Examples

In this section, several numerical examples are given to illustrate the accuracy and effectiveness of the proposed method, and all of them are performed on a computer using programs written in MAPLE 13. In this regard, we have reported in tables and figures the values of the exact solutiony(x), the polynomial approximate solution yN(x), and the absolute error function eN(x)=|y(x)-yN(x)| at any selected points of the given interval [0,1]. It should be noted that, in the second example, we provide an interesting example in which our results are superior with respect to an integral operational matrix method which was based on the Legendre polynomials [5]. Moreover, our method obtains more accurate results with regard to the differentiation of operational matrix which were based on Laguerre, Hermite, and Bessel polynomials and series [3, 14, 15].Example 8 (see [3]).
As the first example, we consider the following variational problem with the exact solutiony(x)=exp(3x)(41)Min∫01‍(y(x)+y′(x)-4exp(3x))2dx,
subject to the boundary conditions
(42)y(0)=1,y(1)=exp(3).
According to (8), the associated Euler Lagrange equation is as follows:
(43)y′′(x)-y(x)-8exp(3x)=0.
We solve this problem by considering notations in the Section 3 by taking N=8,10,11, and 13. In the first figure (i.e., Figure 1), we compare our absolute errors with the Laguerre method [3] (LM) by taking N=8 and 10. Also, in the second figure (i.e., Figure 2), we compare our absolute errors with the Hermite method [14] (HM) by taking N=11 and 13. We should recall that the Laguerre method [3] is a matrix method which has less accuracy with respect to any global polynomial collocation method. This fact is obvious from Figure 1. Moreover, the Hermite method [14] is a collocation method which has low stability property with regard to our method. The interested author can refer to the reference [18] and see the associated ill-conditioned matrix related to the Hermite method.Figure 1
Comparisons of the presented method (PM) and Laguerre method (LM) of Example8.Figure 2
Comparisons of the presented method (PM) and the Hermite method (HM) of Example8.Example 9 (see [5]).
As the second example, we consider the following interesting variational problem with the exact solutiony(x)=(exp(-x)-1)(exp(1)-2exp(2)-2exp(x)+exp(x+1))/4(exp(2)-1)(44)Min∫01‍((y′(x))2+xy′(x)+(y(x))2)dx,
subject to the boundary conditions
(45)y(0)=0,y(1)=14.
According to (8), the associated Euler Lagrange equation is as follows:
(46)y(x)-y′′(x)-12=0.
Again, by using the basic idea of this paper, we solve this problem by taking N=3 and 7. In the third figure (i.e., Figure 3), we compare our absolute errors with the Legendre method [5] (LGM) by taking N=3 and 7. From this figure, one can see the robustness of collocation methods (for instance, the proposed method) with regard to integral operational matrices methods. One of the disadvantages of matrix methods is approximating the known functions by the selected bases (e.g., Legendre), and this may affect the numerical solution seriously and surely the errors increased. However, collocation methods do not approximate the known functions in the basic problem, and the corresponding errors do not increase, and by using a suitable choice of collocation (interpolation) nodes, we can reach to a stable numerical solution.Figure 3
Comparisons of the presented method (PM) and Legendre method (LGM) of Example9.Example 10 (see [5]).
In this example, we consider the following variational problem with the exact solutiony(x)=sinh(-0.4812118250x):
(47)Min∫01‍1+(y(x))2(y′(x))2dx,
subject to the boundary conditions
(48)y(0)=0,y(1)=12.
In this case, the Euler-Lagrange equation is written in the following form:
(49)y′′(x)+y′′(x)(y(x))2-y(x)(y′(x))2=0.
The above BVP seems to be more hard with respect to the other examples of this section. We also solve this BVP by our method and reach to high accurate results with regard to other collocation methods such as the Bessel collocation method [15] and the Taylor collocation method [11]. The absolute error comparisons with the Bessel method (BM) have been depicted in Figure 4, and the absolute error comparisons with the Laguerre method (LM) and Taylor method (TM) have been illustrated in Figure 5. Again, we can see that our results are superior even in comparison with collocation methods.Figure 4
Comparisons of the presented method (PM) and the Bessel method (BM) of Example10.Figure 5
Comparisons of the presented method (PM) and the Laguerre and the Taylor method (LM; TM) of Example10.


## 6.  Conclusions

This paper described an efficient method for finding the extremum of a function over the specified domain. The main goal is to find the solution of an ordinary differential equation (in the BVP form) which arises from the variational problem. Our approach was based on the Bernoulli polynomials and their operational matrix of differentiation together with a set of suitable collocation nodes. The approximation of the solution together with imposing the collocation nodes is utilized to reduce the computation of this problem to some algebraic equations. The method is computationally attractive, and applications are demonstrated through illustrative examples. The obtained results showed that this approach can solve the problem effectively. Moreover, the proposed method could be applied for fractional variational problems [24] with some modifications.



---

*Source: 757206-2013-04-09.xml*