Numerical Solution of High Order Bernoulli Boundary Value Problems

Abstract

For the numerical solution of high order boundary value problems with special boundary conditions a general procedure to determine collocation methods is derived and studied. Computation of the integrals which appear in the coefficients is generated by a recurrence formula and no integrals are involved in the calculation. Several numerical examples are presented to demonstrate the practical usefulness of the proposed method.

1. Introduction

Higher order differential equations arise in a variety of different areas of science, engineering, and technology (see [1, 2]) since they model a wide spectrum of phenomena.

Particularly, the solutions of fifth-order BVPs model viscoelastic flows [3] and the seventh-order BVPs model induction motors with two rotor circuits [4, 5]. Ordinary differential equations of sixth and eighth order arise in modeling instability when an infinite horizontal layer of fluid is heated from below and is subject to the action of rotation [6]. Moreover, high order boundary value problems arise in hydrodynamic, hydromagnetic stability [7], and other branches of applied sciences.

In [8] the authors presented a class of collocation methods for the numerical solution of high order boundary value problems: 𝑦^(𝑛)(𝑥)=𝑓(𝑥,𝐲(𝑥)),𝑥∈𝐼=[𝑎,𝑏],𝐵𝑥,𝑦=𝑔,𝑥∈𝜕𝐼, (1)(2) where 𝑛>1, 𝐲(𝑥)=(𝑦(𝑥),𝑦^′(𝑥),…,𝑦^(𝑞)(𝑥)), 0≤𝑞<𝑛, and 𝐵 is a linear operator on the boundary 𝜕𝐼, 𝑔∈ℝ^𝑛.

The idea in [8] is the following: the differential problem (1)-(2) is written in the following equivalent integral form: 𝑦(𝑥)=𝑃_𝑛−1𝑦,𝑥+^𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑓(𝑡,𝐲(𝑡))𝑑𝑡, (3) where 𝑃_𝑛−1[𝑦,𝑥] is the unique polynomial which satisfies the boundary conditions 𝐵𝑥,𝑃_𝑛−1=𝑔 (4) and 𝐺_𝑛−1(𝑥,𝑡) is a kernel (Green) function. 𝐺_𝑛−1(𝑥,𝑡) is such that 𝐵[𝑥,𝐺_𝑛−1]=0 and it is differentiable under the integral sign such that (3) satisfies (1).

Thus, from (3) and (4) we obtain a collocation polynomial which approximates the solution of problem (1)-(2).

In the present work the authors use this technique to derive collocation methods for the numerical solution of (1) with the particular boundary conditions 𝑦(𝑎)=𝛽₀,𝑦^(𝑘)(𝑏)−𝑦^(𝑘)(𝑎)=𝛽_𝑘+1,𝑘=0,…,𝑛−2 (5) with 𝛽_𝑘, 𝑘=0,…,𝑛−1, being real constants.

Conditions (5) are called the Bernoulli boundary conditions, since they are related to the Bernoulli interpolation problem [9]. They have physical and engineering interpretation [10], but to the authors’ knowledge, they have not been considered previously in the literature.

In [9, 10] the BVP (1)–(5) is considered: in [10] a nonconstructive proof of the existence and uniqueness of solution is given, and in [9] Picard’s method is applied in connection with Newton’s method for the numerical solution of the problem.

The present paper is organized as follows: in Section 2 we summarize some theoretical results on the existence and uniqueness of solution for problem (1)–(5). Then, in Section 3, we present the method for the numerical solution of this type of problems, which produces smooth, global approximations in the form of polynomial functions. In Section 4 we give an a priori estimation of error and, in Section 5, we present some particular cases. In Section 6 we propose an algorithm to compute the numerical solution of (1)–(5) in the nodal points and then, in Section 7, we present some numerical examples of both linear and nonlinear BVPs which confirm the theoretical results.

2. Preliminaries

Let 𝐵_𝑘(𝑥) be the Bernoulli polynomial of degree 𝑘 [11] and let us set 𝑆_𝑘(𝑡)=𝐵_𝑘(𝑡)−𝐵_𝑘(0). (6) Moreover, let ℎ=𝑏−𝑎,Δ𝑓_𝑎^(𝑘)=𝑓^(𝑘)(𝑏)−𝑓^(𝑘)(𝑎). (7) The following theorems hold.

Theorem 1 (see [12]). Let 𝑓∈𝐶^(𝑛+1)[𝑎,𝑏]. Then 𝑓(𝑥)=𝑓(𝑎)+^𝑛_𝑘=1𝑆_𝑘𝑥−𝑎ℎℎ^𝑘−1𝑘!Δ𝑓_𝑎^(𝑘−1)+𝑅_𝑛𝑓,𝑥, (8) where 𝑅_𝑛[𝑓,𝑥] is the remainder term 𝑅_𝑛𝑓,𝑥=^𝑏_𝑎𝐺_𝑛(𝑥,𝑡)𝑓^(𝑛+1)(𝑡)𝑑𝑡 (9) with 𝐺_𝑛(𝑥,𝑡) being Peano’s kernel: 𝐺_𝑛()=1𝑥,𝑡)𝑛!(𝑥−𝑡^𝑛₊−^𝑛_𝑘=1𝑆_𝑘𝑥−𝑎ℎℎ^𝑘−1𝑘×𝑛()𝑘−1𝑏−𝑡^{𝑛−𝑘+1}. (10)

Theorem 2 (see [12]). If 𝑓∈𝐶^(𝑛−1)[𝑎,𝑏], then the polynomial 𝑃_𝑛𝑓,𝑥=𝑓(𝑎)+^𝑛_𝑘=1𝑆_𝑘𝑥−𝑎ℎℎ^𝑘−1𝑘!Δ𝑓_𝑎^(𝑘−1) (11) satisfies the Bernoulli interpolation problem 𝑃_𝑛𝑓,𝑎=𝑓(𝑎),Δ𝑃_𝑛^(𝑘)=𝑃_𝑛^(𝑘)𝑓,𝑏−𝑃_𝑛^(𝑘)𝑓,𝑎=𝑓^(𝑘)(𝑏)−𝑓^(𝑘)(𝑎)=Δ𝑓_𝑎^(𝑘),𝑘=0,…,𝑛−1. (12)

The proof of the existence and uniqueness of solution of (1)–(5) is based on (3) [8], under the hypothesis that the function 𝑓(𝑥,𝐲) satisfies the Lipschitz condition 𝑓𝑥,𝐲₁(𝑥)−𝑓𝑥,𝐲₂(𝑥)≤^𝑞_𝑘=0𝐿_𝑘𝑦₁^(𝑘)(𝑥)−𝑦₂^(𝑘)(𝑥) (13) in a certain domain interval of [𝑎,𝑏]×ℝ^𝑞+1.

3. The Collocation Method

Let 𝑦(𝑥) be the solution of (1)–(5). If 𝑥_𝑖, 𝑖=1,…,𝑚, are 𝑚 distinct points in [𝑎,𝑏] and 𝑦(𝑥)∈𝐶^(𝑛+𝑚)[𝑎,𝑏], using Lagrange interpolation, we get 𝑦^(𝑛)(𝑥)=^𝑚_𝑖=1𝑙_𝑖(𝑥)𝑦^(𝑛)𝑥_𝑖+𝑅_𝑚𝑦,𝑥, (14) where 𝑅_𝑚𝑦,𝑥=𝑥−𝑥₁⋯𝑥−𝑥_𝑚𝑦𝑚!^(𝑛+𝑚)𝜉_𝑥,𝜉_𝑥∈(𝑎,𝑏) (15) and 𝑙_𝑖(𝑡) are the fundamental Lagrange polynomials on the 𝑚 points 𝑥_𝑖.

Inserting (14) into (3), in view of (1), we obtain 𝑦(𝑥)=𝑃_𝑛−1+𝑦,𝑥_𝑚𝑖=1𝑓𝑥_𝑖,𝐲𝑥_𝑖^𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑙_𝑖(𝑡)+𝑑𝑡_𝑏𝑎𝐺_𝑛−1(𝑥,𝑡)𝑅_𝑚𝑦,𝑡𝑑𝑡. (16) Hence the following identity holds: 𝑦(𝑥)=𝑃_𝑛−1𝑦,𝑥+^𝑚_𝑖=1𝑝_𝑛,𝑖(𝑥)𝑓𝑥_𝑖,𝐲𝑥_𝑖+𝑇_𝑛,𝑚𝑦,𝑥, (17) where 𝑝_𝑛,𝑖(𝑥)=^𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑙_𝑖(𝑡)𝑑𝑡,𝑖=1,…,𝑚,𝑇_𝑛,𝑚𝑦,𝑥=^𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑅_𝑚𝑦,𝑡𝑑𝑡. (18)(19)

This suggests defining the polynomials 𝑦_𝑛,𝑚(𝑥)=𝑃_𝑛−1𝑦,𝑥+^𝑚_𝑖=1𝑝_𝑛,𝑖(𝑥)𝑓𝑥_𝑖,𝐲_𝑛,𝑚𝑥_𝑖, (20) where 𝐲_𝑛,𝑚(𝑥)=(𝑦_𝑛,𝑚(𝑥),𝑦^′_𝑛,𝑚(𝑥),…,𝑦^(𝑞)_𝑛,𝑚(𝑥)), 0≤𝑞≤𝑛−1.

Theorem 3. The polynomial of degree 𝑛+𝑚 implicitly defined by (20) satisfies the relations 𝑦_𝑛,𝑚(𝑎)=𝑦(𝑎),𝑦^(𝑘)_𝑛,𝑚(𝑏)−𝑦^(𝑘)_𝑛,𝑚(𝑎)=𝛽_𝑘+1,𝑘=0,…,𝑛−2,𝑦^(𝑛)_𝑛,𝑚𝑥_𝑖=𝑓𝑥_𝑖,𝐲_𝑛,𝑚𝑥_𝑖,𝑖=1,…,𝑚; (21)(22)(23) that is, 𝑦_𝑛,𝑚(𝑥) is a collocation polynomial for (1)–(5) on the nodes 𝑥_𝑖, 𝑖=1,…,𝑚.

Proof. From (18), 𝑝_𝑛,𝑖(𝑎)=𝑝_𝑛,𝑖(𝑏)=0, 𝑖=0,…,𝑚, and thus relations (21) follow from direct computation. To prove (22) we derive 𝐺_𝑛−1(𝑥,𝑡)  𝑘 times, 𝑘=1,…,𝑛−2, with respect to 𝑥, and using the well-known relation [11] 𝐵^′_𝑠(𝑥)=𝑠𝐵_𝑠−1(𝑥), 𝑠>0, we get 𝜕^𝑘𝜕𝑥^𝑘𝐺_𝑛−1()=𝑔𝑥,𝑡₁()=()𝑥,𝑡𝑥−𝑡^𝑛−𝑘()!−𝑛−𝑘_{𝑛𝑗=𝑘}𝐵_𝑗−𝑘𝑥−𝑎ℎℎ^{𝑗−𝑘−1}(𝑏−𝑡)^{𝑛−𝑗−1}!!𝑔𝑗−𝑘𝑛−𝑗−1𝑥≥𝑡₂(𝑥,𝑡)=−^𝑛_𝑗=𝑘𝐵_𝑗−𝑘𝑥−𝑎ℎℎ^{𝑗−𝑘−1}(𝑏−𝑡)^{𝑛−𝑗−1}𝑗−𝑘!𝑛−𝑗−1!𝑥<𝑡. (24) From the property of Bernoulli polynomials 𝐵_𝑠(1)=(−1)^𝑠𝐵_𝑠(0), we have 𝑔₁(𝑎,𝑡)=𝑔₂(𝑏,𝑡); thus 𝑝^(𝑘)_𝑛,𝑖(𝑎)=^𝑏_𝑎𝑔₁(𝑎,𝑡)𝑙_𝑖(𝑡)=𝑑𝑡_𝑏𝑎𝑔₂(𝑏,𝑡)𝑙_𝑖(𝑡)𝑑𝑡=𝑝^(𝑘)_𝑛,𝑖(𝑏). (25) Hence 𝑦^(𝑘)_𝑛,𝑚(𝑏)−𝑦^(𝑘)_𝑛,𝑚(𝑎)=𝑦^(𝑘)(𝑏)−𝑦^(𝑘)(𝑎)+^𝑚_𝑖=1𝑝^(𝑘)_𝑛,𝑖(𝑏)−𝑝^(𝑘)_𝑛,𝑖(𝑎)𝑓𝑥_𝑖,𝐲_𝑛,𝑚𝑥_𝑖=𝑦^(𝑘)(𝑏)−𝑦^(𝑘)(𝑎). (26) From this, (22) follows. Next, by deriving 𝑦_𝑛,𝑚(𝑥)  𝑛 times, we obtain 𝑦^(𝑛)_𝑛,𝑚(𝑥)=𝑃^(𝑛)_𝑛−1𝑦,𝑥+^𝑛−1_𝑘=1𝑝^(𝑛)_𝑛,𝑘(𝑥)𝑓𝑥_𝑘,𝐲_𝑛,𝑚𝑥_𝑘=^𝑛−1_𝑘=1𝑙_𝑘(𝑥)𝑓𝑥_𝑘,𝐲_𝑛,𝑚𝑥_𝑘 (27) and this implies (23).

4. The Error

In what follows for all 𝑦∈𝐶^(𝑞)[𝑎,𝑏] we define the norm ‖𝑦‖=max_{0≤𝑠≤𝑞}{max_{𝑎≤𝑡≤𝑏}|𝑦^(𝑠)(𝑡)|} [14] and the constants 𝐿=^𝑞_𝑘=0𝐿_𝑘,𝑅=max_{𝑎≤𝑥≤𝑏}𝑅_𝑚𝑦,𝑥. (28) Further, we define 𝑄_𝑚=max_{0≤𝑠≤𝑞}max^𝑚_{𝑎≤𝑥≤𝑏𝑖=1}𝑝^(𝑠)_𝑛𝑖(𝑥𝐷_𝑛,𝑠)=max_{𝑎≤𝑥≤𝑏}^𝑛−1_𝑘=𝑠𝐵_𝑘−𝑠((𝑥−𝑎)/ℎ)()!()!,𝑘−𝑠𝑛−𝑘Δ=max_{0≤𝑠≤𝑞}ℎ^{𝑛−𝑠+1}(𝑛−𝑠+1)!+ℎ^{𝑛−𝑠−1}𝐷_𝑛,𝑠. (29)

An a priori estimation of the global error is possible.

Theorem 4. With the previous notations, suppose that 𝐿𝑄_𝑚<1. Then 𝑦−𝑦_𝑛,𝑚≤𝑅Δ1−𝐿𝑄_𝑚. (30)

Proof. By deriving (17) and (20) we get 𝑦^(𝑠)(𝑥)−𝑦^(𝑠)_𝑛,𝑚(𝑥)=^𝑛−1_𝑖=1𝑝^(𝑠)_𝑛𝑖(𝑥)𝑓𝑥_𝑖,𝐲𝑥_𝑖−𝑓𝑥_𝑖,𝐲_𝑛,𝑚𝑥_𝑖+𝜕^𝑠𝜕𝑥^𝑠𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑅_𝑛,𝑚𝑦,𝑡𝑑𝑡. (31) Now, since 𝜕^𝑠𝜕𝑥^𝑠𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑅_𝑛,𝑚=1𝑦,𝑡𝑑𝑡(𝑛−𝑠)!^𝑥_𝑎(𝑥−𝑡)^𝑛−𝑠𝑅_𝑛,𝑚−𝑦,𝑡𝑑𝑡_{𝑛−1𝑘=𝑠}𝐵_𝑘−𝑠𝑥−𝑎ℎℎ^{𝑘−𝑠−1}()!()!×𝑘−𝑠𝑛−𝑘+1_𝑏𝑎(𝑏−𝑡)^{𝑛−𝑘+1}𝑅_𝑛,𝑚𝑦,𝑡𝑑𝑡, (32) we have 𝜕^𝑠𝜕𝑥^𝑠𝑏_𝑎𝐺_𝑛−1(𝑥,𝑡)𝑅_𝑛≤ℎ𝑦,𝑡𝑑𝑡_{𝑛−𝑠+1}(𝑛−𝑠+1)!𝑅+ℎ^{𝑛−𝑠−1}𝑅𝐷_𝑛,𝑠. (33) Thus 𝑦^(𝑠)(𝑥)−𝑦^(𝑠)_𝑛,𝑚(𝑥≤_{𝑛−1𝑖=1})𝑝^(𝑠)_𝑛𝑖(𝑥)^𝑞_𝑘=0𝐿_𝑘𝑦^(𝑠)𝑥_𝑖−𝑦_𝑛^(𝑠)𝑥_𝑖+ℎ^{𝑛−𝑠+1}(𝑛−𝑠+1)!𝑅+ℎ^{𝑛−𝑠−1}𝑅𝐷_𝑛𝑠≤𝐿𝑦−𝑦_𝑛𝑄_𝑚+𝑅Δ (34) and inequality (30) follows.

5. Particular Cases

Now we present explicitly the formulas for some values of 𝑛.

For the computation of 𝑝_𝑛,𝑖(𝑥) we need ∫^𝑥_𝑎𝑡^𝑘𝑙_𝑖(𝑡)𝑑𝑡 and ∫^𝑏_𝑥𝑡^𝑘𝑙_𝑖(𝑡)𝑑𝑡. Letting 𝐹_𝑖1(𝑥)=^𝑥_𝑎𝑙_𝑖(𝑡)𝑑𝑡,𝑀_𝑖1(𝑥)=^𝑏_𝑥𝑙_𝑖(𝑡)𝐹𝑑𝑡,_𝑖𝑘(𝑥)=^𝑥_𝑎𝐹_{𝑖,𝑘−1}(𝑡)𝑑𝑡,𝑀_𝑖𝑘(𝑥)=^𝑏_𝑥𝑀_{𝑖,𝑘−1}(𝑡)𝑑𝑡,𝑘≥2, (35) and integrating by parts 𝑘 times, we obtain ^𝑥_𝑎𝑡^𝑘𝑙_𝑖(𝑡)𝑑𝑡=^𝑘_𝑗=0(−1)^𝑗𝑘!!𝑥𝑘−𝑗^𝑘−𝑗𝐹_𝑖,𝑗+1(𝑥),_𝑏𝑥𝑡^𝑘𝑙_𝑖(𝑡)𝑑𝑡=^𝑘_𝑗=0𝑘!!𝑥𝑘−𝑗^𝑘−𝑗𝑀_𝑖,𝑗+1(𝑥). (36)

5.1. The Fifth-Order Case

Now we consider the case of the fifth-order BVP 𝑦^(𝑣)(𝑥)((𝑥[],𝑦(0)=𝑓𝑥,𝐲)),𝑥∈0,1=𝛽₀,𝑦^(𝑘)(1)−𝑦^(𝑘)(0)=𝛽_𝑘+1,𝑘=0,…,3. (37) In this case Green’s function is 𝐺₄()=1𝑥,𝑡𝑡4!⁴(1−𝑥)+𝑡³2𝑥²−2𝑥+𝑡²−2𝑥³+3𝑥²𝑥−𝑥+𝑡⁴−2𝑥³+𝑥²1𝑡≤𝑥4!−𝑡⁴𝑥+𝑡³2𝑥²+2𝑥−𝑡²2𝑥³+3𝑥²𝑥+𝑥+𝑡⁴+2𝑥³+𝑥²−𝑥⁴𝑥≤𝑡,4!𝑝_5,𝑖(𝑥)=𝑥⁴−2𝑥³+𝑥²𝐹_𝑖2(𝑥)−𝑀_𝑖2(𝑥)−2𝑥2𝑥²−3𝑥+1𝐹_𝑖3(𝑥)+𝑀_𝑖3(𝑥(𝐹)+12𝑥𝑥−1)_𝑖4(𝑥)−𝑀_𝑖4(𝑥()𝐹)+241−𝑥_𝑖5(𝑥)−24𝑥𝑀_𝑖5(𝑥). (38) Hence 𝑦_5,𝑚(𝑥)=𝑃₄𝑦,𝑥+^𝑚_𝑖=1𝑝_5,𝑖(𝑥)𝑓𝑥_𝑖,𝐲_5,𝑚𝑥_𝑖. (39) By deriving (44) we get 𝑦^(𝑠)_5,𝑚(𝑥)=𝑃₄^(𝑠)𝑦,𝑥+^𝑚_𝑖=1𝑝^(𝑠)_5,𝑖(𝑥)𝑓𝑥_𝑖,𝐲_5,𝑚𝑥_𝑖𝑠=1,…,5, (40) where 𝑝^(𝑠)_5,𝑖(𝑥) can be easily computed using the same technique as for 𝑝_5,𝑖(𝑥).

5.2. The Seventh-Order Case

Consider 𝑦^{(𝑣𝑖𝑖)}(𝑥)((𝑥[],𝑦(0)=𝑓𝑥,𝐲)),𝑥∈0,1=𝛽₀,𝑦^(𝑘)(1)−𝑦^(𝑘)(0)=𝛽_𝑘+1,𝑘=0,…,5. (41)

In this case Green’s function is 𝐺₆()=1𝑥,𝑡𝑡6!⁶(1−𝑥)+3𝑡⁵𝑥²−𝑥+5𝑡⁴−𝑥³+32𝑥²−𝑥2+5𝑡³5𝑥⁴−2𝑥³+𝑥²+𝑡²−3𝑥⁵+152𝑥⁴−5𝑥³+𝑥2+𝑡𝑥⁶−3𝑥⁵+52𝑥⁴−𝑥²2𝑡≤𝑥−𝑡⁶𝑥+3𝑡⁵𝑥²+𝑥+5𝑡⁴−𝑥³−32𝑥²−𝑥2+5𝑡³5𝑥⁴+2𝑥³+𝑥²+𝑡²−3𝑥⁵+152𝑥⁴−5𝑥³+𝑥2+𝑡𝑥⁶−3𝑥⁵+52𝑥⁴−𝑥²2−𝑥⁶𝑥≤𝑡. (42) Hence 6!𝑝_7,𝑖(𝑥)=𝑥⁶−3𝑥⁵+52𝑥⁴−𝑥²2𝐹_𝑖2(𝑥)−𝑀_𝑖2(𝑥)−𝑥6𝑥⁴−15𝑥³+10𝑥²−1𝐹_𝑖3(𝑥)+𝑀_𝑖3(𝑥())+30𝑥−1²𝐹_𝑖4(𝑥)−𝑀_𝑖4(𝑥)−60𝑥2𝑥²−3𝑥+1𝐹_𝑖5(𝑥)+𝑀_𝑖5(𝑥(𝐹)+360𝑥𝑥−1)_𝑖6(𝑥)−𝑀_𝑖6(𝑥()𝐹)+6!1−𝑥_𝑖7(𝑥)−6!𝑥𝑀_𝑖7(𝑥),𝑦_7,𝑚(𝑥)=𝑃₆𝑦,𝑥+^𝑚_𝑖=1𝑝_7,𝑖(𝑥)𝑓𝑥_𝑖,𝐲_7,𝑚𝑥_𝑖. (43)(44)

5.3. Order 𝑛=8,9,10

For 𝑛=8, we have 7!𝑝_8,𝑖(𝑥)=𝑥⁸4−𝑥⁷+76𝑥⁶−7𝑥12⁴+𝑥²6𝐹_𝑖1(𝑥)+𝑀_𝑖1(𝑥+𝑥₇)−72𝑥⁶+72𝑥⁵−76𝑥³+𝑥6𝐹_𝑖2(𝑥)−𝑀_𝑖2(𝑥𝑥)−7𝑥⁵−3𝑥⁴+52𝑥³−𝑥2𝐹_𝑖3(𝑥)+𝑀_𝑖3(𝑥)−7𝑥6𝑥⁴−15𝑥³+10𝑥²−1𝐹_𝑖4(𝑥)−𝑀_𝑖4(𝑥)+210𝑥²(𝑥−1)²𝐹_𝑖5(𝑥)+𝑀_𝑖5(𝑥)−420𝑥2𝑥²−3𝑥+1𝐹_𝑖6(𝑥)−𝑀_𝑖6(𝑥(𝐹)+2520𝑥𝑥−1)_𝑖7(𝑥)+𝑀_𝑖7(𝑥()𝐹)+7!1−𝑥_𝑖8(𝑥)+7!𝑥𝑀_𝑖8(𝑥). (45)

For 𝑛=9, we get 8!𝑝_9,𝑖(𝑥)=𝑥⁸−4𝑥⁷+143𝑥⁶−73𝑥⁴+23𝑥²𝐹_𝑖2(𝑥)−𝑀_𝑖2(𝑥)+4𝑥−2𝑥⁶+7𝑥⁵−7𝑥⁴+73𝑥²−𝑥3𝐹_𝑖3(𝑥)+𝑀_𝑖3(𝑥)+28𝑥²2𝑥⁴−6𝑥³+5𝑥⁴−1𝐹_𝑖4(𝑥)−𝑀_𝑖4(𝑥)+56𝑥−5𝑥⁴+15𝑥³−10𝑥²+1𝐹_𝑖5(𝑥)+𝑀_𝑖5(𝑥)+1680𝑥²(𝑥−1)²𝐹_𝑖6(𝑥)−𝑀_𝑖6(𝑥)+3360𝑥−2𝑥²+3𝑥−1𝐹_𝑖7(𝑥)+𝑀_𝑖7(𝑥(𝐹)+20160𝑥𝑥−1)_𝑖8(𝑥)−𝑀_𝑖8(𝑥()𝐹)+8!1−𝑥_𝑖9(𝑥)−8!𝑥𝑀_𝑖9(𝑥). (46)

For 𝑛=10, we obtain 9!𝑝_10,𝑖(𝑥)=𝑥¹⁰5−𝑥⁹+32𝑥⁸−75𝑥⁶+𝑥⁴−3𝑥10³×𝐹_𝑖1(𝑥)+𝑀_𝑖1(𝑥𝑥)−𝑥⁸−92𝑥⁷+6𝑥⁶−215𝑥⁴+2𝑥²−3×𝐹10_𝑖2(𝑥)−𝑀_𝑖2(𝑥)+3𝑥²3𝑥⁶−12𝑥⁵+14𝑥⁴−7𝑥²+2𝐹_𝑖3(𝑥)+𝑀_𝑖3(𝑥)+12𝑥−6𝑥⁶+21𝑥⁵−21𝑥⁴+7𝑥²−1𝐹_𝑖4(𝑥)−𝑀_𝑖4(𝑥)+252𝑥²2𝑥⁴−6𝑥³+5𝑥²−1𝐹_𝑖5(𝑥)+𝑀_𝑖5(𝑥)+504𝑥−6𝑥⁴+15𝑥³−10𝑥²+1𝐹_𝑖6(𝑥)−𝑀_𝑖6(𝑥)+15120𝑥²(𝑥−1)²𝐹_𝑖7(𝑥)+𝑀_𝑖7(𝑥)+30240𝑥−2𝑥²+3𝑥−1𝐹_𝑖8(𝑥)−𝑀_𝑖8(𝑥(𝐹)+181440𝑥𝑥−1)_𝑖9(𝑥)+𝑀_𝑖9(𝑥()𝐹)+9!1−𝑥_𝑖10(𝑥)+9!𝑥𝑀_𝑖10(𝑥). (47)

6. Algorithms

To calculate the approximate solution of problem (1)–(5) by (20) at 𝑥∈[𝑎,𝑏], we need the values 𝑦_𝑗^(𝑘)=𝑦^(𝑘)_𝑛,𝑚(𝑥_𝑗), 𝑗=1,…,𝑚, and 𝑘=0,…,𝑞. These values can be calculated by solving the following system: 𝑦_𝑖^(𝑘)=𝑃^(𝑘)_𝑛−1𝑦,𝑥_𝑖+^𝑚_𝑗=1𝑝^(𝑘)_𝑛𝑗𝑥_𝑖𝑓𝑥_𝑗,𝐲_𝑗𝑖=1,…,𝑚,𝑘=0,…,𝑞 (48) with 𝐲_𝑗=(𝑦_𝑗,𝑦^′_𝑗,…,𝑦_𝑗^(𝑞)), 0≤𝑞≤𝑛−1.

To solve it, if we put 𝑌=𝑌₀,…,𝑌_𝑞^𝑇,𝑌_𝑘=𝑦₁^(𝑘),…,𝑦_𝑚^(𝑘),𝐹(𝑌)=𝐹_𝑚,…,𝐹_𝑚^𝑇∈ℝ^𝑞+1,𝐹_𝑚=𝑓₁,…,𝑓_𝑚,𝑓_𝑖=𝑓𝑥_𝑖,𝐲_𝑖,𝐶=𝐵₀,…,𝐵_𝑞^𝑇,𝐵_𝑘=𝑃_𝑚^(𝑘)𝑦,𝑥₁,…,𝑃^(𝑘)_𝑛−1𝑦,𝑥_𝑚,𝐴𝐴=⁰0⋯00⋱⋮⋮⋱00⋯0𝐴_𝑞,𝐴_𝑘=𝑝^(𝑘)_1,1⋯𝑝^(𝑘)_1,𝑚𝑝⋮⋮_(𝑘)𝑚,1⋯𝑝^(𝑘)_𝑚,𝑚 (49) with 𝑝^(𝑘)_𝑖,𝑗=𝑝^(𝑘)_𝑛,𝑗(𝑥_𝑖), 𝑘=0,…,𝑞, we write (48) as 𝑌−𝐴𝐹(𝑌)=𝐶 (50) or, equivalently, 𝑌=𝐺(𝑌), where 𝐺(𝑍)=𝐴𝐹(𝑍)+𝐶. (51)

For the existence and the uniqueness of the solution of (50) the following result holds.

Proposition 5. Let 𝐿 be defined as in (28). If 𝑇=𝐿‖𝐴‖_∞<1, the system (50) has a unique solution which can be calculated by an iterative method 𝑌^(𝜈+1)=𝐺𝑌^(𝜈),𝜈=1,2,… (52) with a fixed 𝑌⁽⁰⁾∈ℝ^𝑠, 𝑠=𝑚(𝑞+1), and 𝐺 defined as in (51). Moreover, if 𝑌 is the exact solution of the system, 𝑌^(𝜈)−𝑌_∞≤𝑇^𝜈𝑌1−𝑇⁽¹⁾−𝑌⁽⁰⁾_∞. (53)

Proof. If 𝑉=(𝑉₀,…,𝑉_𝑞)^𝑇, 𝑉_𝑘=(𝑣₁^(𝑘),…,𝑣_𝑚^(𝑘)) and 𝑊=(𝑊₀,…,𝑊_𝑞)^𝑇, 𝑊_𝑘=(𝑤₁^(𝑘),…,𝑤_𝑚^(𝑘)), then ‖𝐺(𝑉)−𝐺(𝑊)‖_∞≤‖𝐴‖_∞𝐿‖𝑉−𝑊‖_∞; hence 𝐺 is contractive. Thus the result follows from the well-known contraction mapping theorem.

To calculate the elements 𝐴₀,…,𝐴_𝑞 of the matrix 𝐴 we need the values 𝐹_𝑖𝑠(𝑥_𝑗) and 𝑀_𝑖𝑠(𝑥_𝑗), 𝑠=1,…,𝑛, 𝑖,𝑗=1,…𝑚, where 𝐹_𝑖𝑠(𝑥) and 𝑀_𝑖𝑠(𝑥) are defined in (35).

Since 𝑙_𝑖(𝑡)=∏^𝑚_{𝑘=1,𝑘≠𝑖}((𝑡−𝑥_𝑘)/(𝑥_𝑖−𝑥_𝑘)), it suffices to compute ^{𝑥𝑗=𝑡𝑘}_𝑐^𝑡𝑘−1_𝑐⋯^𝑡1_𝑐𝑟_𝑚,𝑖(𝑡)𝑑𝑡𝑑𝑡₁⋯𝑑𝑡_𝑘−1, (54) where 𝑐=𝑎 or 𝑐=𝑏, 𝑟_0,0(𝑡)=1, and 𝑟_𝑚,𝑖(𝑡)=𝑡−𝑥₁⋯𝑡−𝑥_𝑖−1𝑡−𝑥_𝑖+1⋯𝑡−𝑥_𝑚𝑖=1,2,…,𝑚. (55) Let us define 𝑔^(𝑖)_0,1,𝑐(𝑥)=𝑥−𝑐 (56) and, for 𝑠=1,…,𝑚−1, 𝑔^(𝑖)_{𝑠,𝑗,𝑐}(𝑥)=^{𝑥=𝑡𝑗}_𝑐^𝑡𝑗−1_𝑐⋯^𝑡1_𝑐𝑡−𝑧₁^(𝑖)𝑡−𝑧₂^(𝑖)⋯𝑡−𝑧_𝑠^(𝑖)𝑑𝑡𝑑𝑡₁⋯𝑑𝑡_𝑗−1, (57) where 𝑧_𝑗^(𝑖)=𝑥_𝑗𝑥if𝑗<𝑖_𝑗+1if𝑗≥𝑖𝑗=1,…,𝑚−1. (58) We can easily compute 𝑔^(𝑖)_0,𝑗,𝑐(𝑥)=(𝑥−𝑐)^𝑗.𝑗! (59) For the computation of (57) we use the recursive algorithm [15] 𝑔^(𝑖)_{𝑠,𝑗,𝑐}(𝑥)=𝑥−𝑧_𝑠^(𝑖)𝑔^(𝑖)_{𝑠−1,𝑗,𝑐}(𝑥)−𝑗𝑔^(𝑖)_{𝑠−1,𝑗+1,𝑐}(𝑥). (60) Thus, if 𝑊_𝑖=∏^𝑚_{𝑘=1,𝑘≠𝑖}(𝑥_𝑖−𝑥_𝑘), we get 𝐹_𝑖𝑘𝑥_𝑗=𝑔^(𝑖)_{𝑚−1,𝑘,𝑎}𝑥_𝑗𝑊_𝑖,𝑀_𝑖𝑘𝑥_𝑗=(−1)^𝑘𝑔^(𝑖)_{𝑚−1,𝑘,𝑏}𝑥_𝑗𝑊_𝑖. (61)

7. Numerical Examples

Now we present some numerical results obtained by applying method (20) to find numerical approximations of the solutions of some test problems. As the true solutions are known, we considered the error function 𝐸(𝑥)=|𝑦(𝑥)−𝑦_𝑛,𝑚(𝑥)|. To solve the nonlinear system (48) we used the so-called modified Newton method [16] (the same Jacobian matrix is used for more than one iteration) and algorithm (60) for the computation of the entries of the matrix. Equidistant points are used as nodal points. Analogous results are obtained in the considered examples by using as nodes the zeros of Chebyshev polynomials of first and second kind.

Example 1. Consider the following 𝑦^(𝑣)(𝑥)=𝑦(𝑥)−(15+10𝑥)𝑒^𝑥[],𝑦(0)(1)𝑦,𝑥∈0,1=0,𝑦=0,_′(1)−𝑦^′(0)=−(𝑒+1),𝑦^′′(1)−𝑦^′′(0)𝑦=−4𝑒,_′′′(1)−𝑦^′′′(0)=3−9𝑒, (62) with solution 𝑦(𝑥)=𝑥(1−𝑥)𝑒^𝑥. Figure 1 shows the graph of the error function 𝐸(𝑥) for two different values of 𝑚.

(a)

(b)

Figure 1 

Error function of problem (62) for 𝑚=4 (a) and for 𝑚=6 (b).

Example 2. Consider the following 𝑦^(𝑣)(𝑥)=−24𝑒^−5𝑦+481+𝑥⁵[],𝑦(0)(1)𝑦,𝑥∈0,1=0,𝑦=log2,_′(1)−𝑦^′(0)1=−2,𝑦^′′(1)−𝑦^′′(0)=34,𝑦^′′′(1)−𝑦^′′′(0)7=−4, (63) with solution 𝑦(𝑥)=log(𝑥+1). The graph of 𝐸(𝑥), for two different numbers of nodes, is plotted in Figure 2.

(a)

(b)

Figure 2 

Error function of problem (63) for 𝑚=4 (a) and for 𝑚=7 (b).

Example 3. Consider 𝑦^{(𝑣𝑖𝑖)}(𝑥)=−𝑦−𝑒^𝑥2𝑥²[],𝑦(0)(1)𝑦+12𝑥+35,𝑥∈0,1=𝑦=0,_′(1)−𝑦^′(0)=−(𝑒+1),𝑦^′′(1)−𝑦^′′(0)𝑦=−4𝑒,_′′′(1)−𝑦^′′′(0)(),𝑦=31−3𝑒_(𝑖𝑣)(1)−𝑦^(𝑖𝑣)(0)(),𝑦=81−2𝑒_(𝑣)(1)−𝑦^(𝑣)(0)=5(3−5𝑒), (64) with solution 𝑦(𝑥)=𝑥(1−𝑥)𝑒^𝑥. Figure 3 shows the graph of 𝐸(𝑥).

(a)

(b)

Figure 3 

Error function of problem (64) for 𝑚=4 (a) and for 𝑚=8 (b).

Note that the equation in (64) is the same as that in Example 2 of [1], but the boundary conditions are different.

Example 4. Consider 𝑦^(𝑖𝑥)(𝑥)=𝑦(𝑥)−9𝑒^𝑥[],𝑦(0)(1)𝑦,𝑥∈0,1=1,𝑦=0,_′(1)−𝑦^′(0)=−𝑒,𝑦^′′(1)−𝑦^′′(0)𝑦=1−2𝑒,_′′′(1)−𝑦^′′′(0)𝑦=2−3𝑒,_(𝑖𝑣)(1)−𝑦^(𝑖𝑣)(0)𝑦=3−4𝑒,_(𝑣)(1)−𝑦^(𝑣)(0)𝑦=4−5𝑒,_(𝑣𝑖)(1)−𝑦^(𝑣𝑖)(0)𝑦=5−6𝑒,_{(𝑣𝑖𝑖)}(1)−𝑦^{(𝑣𝑖𝑖)}(0)=6−7𝑒, (65) with solution 𝑦(𝑥)=(1−𝑥)𝑒^𝑥. Figure 4 shows the graph of 𝐸(𝑥).

(a)

(b)

Figure 4 

Error function of problem (65) for 𝑚=4 (a) and for 𝑚=8 (b).

The equation in (65) is the same as that in [1, 2], but the boundary conditions are different.

Example 5. Consider 𝑦^(𝑥)(𝑥)=𝑒^−𝑥𝑦²(𝑥)[],𝑦(0)𝑦,𝑥∈0,1=1,_(𝑘)(1)−𝑦^(𝑘)(0)=𝑒−1,𝑘=0,…,8 (66) with solution 𝑦(𝑥)=𝑒^𝑥. Figure 5 shows the graph of 𝐸(𝑥).

(a)

(b)

Figure 5 

Error function of problem (66) for 𝑚=3 (a) and for 𝑚=6 (b).

Note that the equation in (66) is the same as that in [2], but the boundary conditions are different. The conditions in [2] are the so-called Lidstone-type conditions. Problems of this type have been analyzed in [13] using a similar technique.

8. Conclusions

This paper presents a class of collocation methods for 𝑛th order differential equations with Bernoulli boundary conditions. For two positive integers 𝑛,𝑚 a polynomial of degree 𝑛+𝑚 approximating the exact solution is given explicitly. Numerical experiments support theoretical results. Further developments can be done, concerning particularly numerical estimates of the error.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.