Cubic Lidstone-Spline for numerical solution of BVPs

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation (

)

– www.elsevier.com/locate/matcom

Original articles

Cubic Lidstone-Spline for numerical solution of BVPs Francesco A. Costabile a,∗ , Maria Italia Gualtieri a , Giada Serafini b a Department of Mathematics and Computer Science, University of Calabria, Italy b Department of Mathematics, Computer Science and Economics, University of Basilicata, Italy

Received 7 December 2015; received in revised form 14 January 2017; accepted 19 January 2017 Available online xxxx

Abstract A new collocation method based on cubic Lidstone Splines is introduced for solving second order BVPs. It derives directly from piecewise Lidstone polynomials of degree 3 by requiring the continuity of the first derivative at the nodal points. For equally spaced nodes it reduces to a classical finite difference method of second order. The estimation of local and global error is given. Finally we solve some numerical problems and we compare the results with those of other methods. c 2017 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Boundary value problem; Cubic spline; Lidstone polynomials

1. Introduction This work deals with the global approximation of the solution of the following general two-point Boundary Value Problem (BVP)  ′′  y (x) = f (x, y(x)), x ∈ I = [a, b] y(a) = ya (1)  y(b) = yb . Problems of this kind arise in engineering and other branches of the physical sciences, for instance, boundary layer theory in fluid mechanics, heat power transmission, space technology, control and optimization theory. We assume that f (x, y) is a real defined and continuous function over the strip S = [a, b] × R. Moreover we assume that over the strip S the two following inequalities hold true – a Lipschitz constant L exists s.t. for any y1 , y2 ∈ R | f (x, y1 ) − f (x, y2 )| ≤ L|y1 − y2 |; – f y (x, y) is continuous and satisfies f y (x, y) ≥ 0, a ≤ x ≤ b, −∞ < y < ∞. ∗ Corresponding author.

E-mail addresses: [email protected] (F.A. Costabile), [email protected] (M.I. Gualtieri), [email protected] (G. Serafini). http://dx.doi.org/10.1016/j.matcom.2017.01.006 c 2017 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

2

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

This problem, of class M [20], has a unique solution and several numerical methods are known for solving it [8]. Among them we recall shooting methods, which transform the given boundary value problem into an initial value problem with estimated parameters, then adjust the parameters iteratively to reproduce the given boundary values [8]. Another approach for solving BVPs consists of spectral methods, the key elements of which are the basic functions (also called trial functions) for a truncated series expansion of the solution and the choice of any of the three discrete procedures most commonly used, that is, the Galerkin, collocation and Tau scheme [10]. Another class of methods for BVPs is the so-called “finite elements”, that are philosophically similar to spectral algorithms. Finite element methods chop the interval [a, b] into a number of subintervals, and choose the basic functions which are polynomials of fixed degree, non-zero only over a couple of subintervals. Finally, classical finite difference methods [24] approximate the solution, y(x), of a BVP by a sequence of overlapping polynomials which interpolate y(x) at a set of grid-points. Boyd [10] observed that “when many decimal places of accuracy are needed, the contest between pseudospectral algorithms and finite difference and finite element methods is not an even battle but a rout: pseudospectral methods win hands-down”. This is one of the reasons why physicists and quantum chemists have always preferred spectral methods. The possibility of using spline polynomials to obtain a smooth approximate solution of (1) was briefly discussed for the first time by Ahlberg [3] and, successively, by many other authors [5–7,4,9,11,12,17,18,21–23,26–28,30–33, and the references therein]. In particular, in [6] the authors obtain smooth approximations by using quadratic splines for solving problems of type (1); moreover the use of higher order collocation spline functions and collocation for solving second order BVPs has been considered in [27,32,33] and the references therein. In this note we introduce a piecewise collocation method based on cubic spline functions. We will consider the general nonlinear problem of class M on non-uniform meshes. We use the Lidstone polynomials as a basis and therefore the method is named Collocation Cubic Lidstone-Spline. The paper is organized as follows: in Section 2 we give the preliminaries; in Section 3 we consider our method; in Section 4 local and global error are introduced and in Section 5 we consider numerical examples. Finally we give some conclusions. 2. Preliminaries Lidstone polynomials [25] are defined by  Λ0 (x) = x, Λ′′ (x) = Λk−1 (x) x ∈ [0, 1], k ≥ 1.  k Λk (0) = Λk (1) = 0

(2)

For any k ≥ 0, Λk (x) is a polynomial of degree 2k + 1. Setting h := b − a L n ( f, x) =

n−1 

      x −a b−x f (2k) (a) + Λk f (2k) (b) , h 2k Λk h h k=0

(3)

is the Lidstone interpolating polynomial of a given function f ∈ C 2n ([a, b]) , n > 0 [2]. Theorem 1 (Lidstone Interpolation [2,1]). Let f (x) ∈ C 2n ([a, b]), with n > 0, then f (x) = L n ( f, x) + E 2n (x)h 2n f (2n) (ξx ),

∀x ∈ [a, b], ξx ∈ (a, b)

(4)

where E 2n (x) is the Euler polynomial of degree 2n in [a, b]. Furthermore, f (2k) (a) = L n(2k) ( f, a)

and

f (2k) (b) = L n(2k) ( f, b) ,

k = 0, 1, . . . , n − 1.

Let a ≡ x0 < x1 < x2 < · · · < xn−1 < xn ≡ b. Theorem 2 (Spline Interpolation [13]). A family of interpolant cubic spline, s (x), on set points (xν , yν ), ν = 0, . . . , n, is given by     (xν−1 − x)3 (xν − x)3 x − xν−1 h 2ν xν − x h 2ν pν (x) = Mν + Mν−1 + yν − Mν + yν−1 − Mν−1 , (5) −6h ν 6h ν hν 6 hν 6

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

3

–

x ∈ Iν = [xν−1 , xν ], ν = 0, . . . , n, where ν = 1, . . . , n

h ν = xν − xν−1 yν = s(xν ) Mν = s ′′ (xν )

(6)

ν = 0, . . . , n,

iff the following relations hold yν − yν−1 yν+1 − yν h ν + h ν+1 hν h ν+1 − = Mν−1 + Mν + Mν+1 , h ν+1 hν 6 3 6

ν = 1, . . . , n − 1.

(7)

Proof. The proof is similar to the construction of the natural interpolant cubic spline [13]. Remark 1. We stress explicitly that (7) is a linear system of n − 1 equations in the n + 1 unknowns, M0 , . . . , Mn . If we define L n,ν ( f, x) as the polynomial L n ( f, x) in [a, b] ≡ Iν , ν = 1, . . . , n, we observe that (5) and (7) can be obtained by (3) with n = 2, and L n,ν ( f, x) ≡ pν (x), assuming the continuity of the first derivative at the nodal points L ′n,ν ( f, xν ) = L ′n,ν+1 ( f, xν ) .

(8)

3. The Collocation Cubic Lidstone-Spline Denoting by y(x) the solution of the BVP in (1), we will approximate it by a spline of collocation on the grid-points {xν }, ν = 0, . . . , n. To this end we have Theorem 3 (Main Theorem). Let y(x) be the solution of the problem in (1) and for any x ∈ Iν = [xν−1 , xν ], ν = 1, . . . , n, let s(x) ∈ C 2 ([a, b]) defined as   (xν−1 − x)3 (xν − x)3 x − xν−1 h 2ν s(x) := p ν (x) = f ν + f ν−1 + yν − f ν −6h ν 6h ν hν 6   2 xν − x h + yν−1 − f ν−1 ν (9) hν 6 with s (xν ) = yν ,

s ′′ (xν ) = f ν ,

ν = 0, . . . , n − 1

and yν − yν−1 hν h ν + h ν+1 h ν+1 yν+1 − yν − = f ν−1 + fν + f ν+1 h ν+1 hν 6 3 6

ν = 1, . . . , n − 1.

(10)

If in (9) we consider f ν = f (xν , yν ) , ν = 0, . . . , n, with f as in (1), s (x) is a cubic collocation spline for (1) and we can assume y (xν ) ≈ s (xν ) = yν , y (x0 ) = y0 ≡ ya ,

ν = 1, . . . , n − 1

(11)

y (xn ) = yn ≡ yb .

Proof. It follows by Theorem 2 after observing that Mν = s ′′ (xν ) = f ν ,

s(xν ) = p ν (xν ) = yν

ν = 0, . . . , n.

Furthermore, (11) is justified by the collocation principle. Remark 2. We stress that, in the general case, (1) is a nonlinear BVP, therefore (10) is a nonlinear system of n − 1 equations in n − 1 unknowns, since f 0 = f (x0 , y0 ) and f n = f (xn , yn ) can be obtained by the boundary conditions in (1).

4

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

After Theorem 3 we call s(x) in (9) the Collocation Cubic Lidstone-Spline for (1). In order to compute the approximate solution, we have to solve the system in (10) which can be rewritten in the following matrix form AY = −B F(Y ) + K ,

(12)

where  h2 + h1  h2h1   1  −  h2 A=  . ..     0

1 h2 h3 + h2 h3h2 .. .

···

0



1 − h3 .. .

.. .

          

−

···

−

1

−

1

h n−1 h n + h n−1 h n h n−1  0   ..   .   h n−1    6  h n−1 + h n 3

h n−1 h + h h 1 2 2 ···  3 6  h2 + h3 h3  h2   6 3 6 B= ..  .. ..  . . .   h n−1 0 ··· 6  T Y = y1 , . . . , yn−1 ]T , F(Y ) = [ f (x1 , y1 ), . . . , f (xn−1 , yn−1 )   h1 hn y0 yn T K = M0 − , 0, . . . , 0, Mn − 6 h1 6 hn ′′ Mν = s (xν ) = f (xν , yν ), ν = 0, . . . , n. Theorem 4. The matrix A of the system in (12) is invertible. Proof. We observe that A is a square, tridiagonal and diagonally predominant matrix. Furthermore, the conditions ai,i−1 ̸= 0, i = 2, 3, . . . , n and ai,i+1 ̸= 0, i = 1, 2, . . . , n − 1 assure the irreducibility of the matrix A. Therefore, we can state that A is nonsingular. Theorem 5. Under the assumption ∥A−1 B∥L < 1, where L is a Lipschitz constant of f and ∥ · ∥ is the infinity norm, the system (12) is uniquely solvable. Proof. Since by Theorem 4 the matrix A is nonsingular, the system (12) can be written in the equivalent form Y = −A−1 B F(Y ) + A−1 K .

(13)

Setting G(Y ) = −A−1 B F(Y ) + A−1 K , we have Y = G(Y ). Under the considered hypotheses, it follows that G is a contractive map and the assertion follows by fixed point theorem. For the effective calculation of the solution, we observe that if f (x, y) is linear, i.e. f (x, y) = g(x)y(x) + k(x), the system (12) becomes AY = K

(14)

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

5

where  1 hi   g(xi−1 ) −   6 h  i   h i + h i+1 h i+1 + h i + A = (a i, j ) = g(xi )  3 h i+1 h i     1 h  g(xi+1 ) i+1 − 6 h i+1

j = i − 1; j = i;

(15)

j =i +1

and 

   h1 h2 h1 + h2 − k(x2 ) − k(x1 )   6 3 6       h2 h2 + h3 h3   −k(x1 ) − k(x2 ) − k(x3 )   6 3 6     ..     . .   K =   h k + h k+1 hk h k+1   −k(xk−1 ) − k(xk ) − k(xk+1 )   6 3 6     ..     .         h n−1 h n−1 + h n hn 1 hn −k(xn−2 ) − k(xn−1 ) − k(xn ) + − g(xn ) yn 6 3 6 hn 6 

1 h1 − g(x0 ) h1 6



y0 − k(x0 )

Since the matrix A of the linear system in (14) is tridiagonal, an efficient algorithm can be used to solve it. In the general case the system in (12) is nonlinear and it can be written in the form Y = B F(Y ) + K ,

(16)

where B = A−1 B, K = A−1 K , and we can apply an iterative Newton-type method, solving a linear tridiagonal system at each iteration. 4. Local and global error estimation First we observe that for equally spaced nodes, i.e. h ν = h for ν = 1, . . . , n, the system in (10) becomes   2 1 2 1 yν+1 − 2yν + yν−1 − h f (xν−1 , yν−1 ) + f (xν , yν ) + f (xν+1 , yν+1 ) = 0 6 3 6 which can be interpreted as a classical linear difference finite method [20]. In this case we have Theorem 6. For the linear operator L [y, h] := y(x + h) − 2y(x) + y(x − h) − h 2



 1 ′′ 2 1 y (x + h) + y ′′ (x) + y ′′ (x − h) , 6 3 6

associated to (17), assuming y ∈ C 4 ([a, b]), we have L [y, h] = −

1 4 IV h y (ξ ), 12

ξ ∈ (x − h, x + h).

Proof. The theorem follows by the classic Taylor formula. Remark 3. In the classical formula yn+1 − 2yn + yn−1 = h 2 f (xn , yn ) [20] L [y, h] := y(x + h) − 2y(x) + y(x − h) − h 2 y ′′ (x)

(17)

6

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

and we have L [y, h] =

1 4 IV h y (ξ ), 12

ξ ∈ (x − h, x + h).

In a similar way Theorem 7. For the linear operator    y, h ν , h ν+1 = h ν [y(x + h ν ) − y(x)] − h ν+1 [y(x) − y(x − h ν )] L   1 h ν + h ν+1 ′′ 1 − h ν h ν+1 h ν y ′′ (x − h) + y (x) + h ν+1 y ′′ (x + h) , 6 3 6 associated to (10), assuming y ∈ C 4 ([a, b]), we have      y, h ν , h ν+1 = − 1 h ν h ν+1 h 3 + h 3 y I V (ξ ), L ν ν+1 24

(18)

ξ ∈ (xν−1 , xν+1 ).

Remark 4. In particular, we observe that for h ν = h, ν = 1, . . . , n, we have [y, h, h] = hL [y, h] . L [y, h, h] will be the local error of discretization in (10). In what follows Tν [h] = L Remark 5. Furthermore, we observe that (10) can be considered a discretization of the operator y ′′ (x) on a nonuniform grid of points, if we assume yν = y(xν ) and f ν = y ′′ (xν ). In order to obtain a global error estimate, if Y = [y(x1 ), . . . , y(xn−1 )]T is the vector of exact solution, T = [T1 (h), . . . , Tn−1 (h)]T the vector of local error and Y = [y1 , . . . , yn−1 ]T the vector of approximate solution, the following theorem holds Theorem 8. Under the assumption ∥A−1 B∥L < 1, where L is a Lipschitz constant of f , we have ∥Y − Y ∥ ≤

∥A−1 ∥ ∥T ∥ . 1 − ∥A−1 B∥L

Proof. From (10) we have AY = −B F(Y ) + K

(19)

and by (18) AY = −B F(Y ) + T + K .

(20)

By subtracting (19) and (20), we have A(Y − Y ) = −B(F(Y ) − F(Y )) − T.

(21)

Since f (x, y) satisfies the Lipschitz condition with constant L, the theorem easily follows. 5. Numerical results We test the efficiency of the Collocation Cubic Lidstone-Spline method (shortly CCLS) in some numerical problems (linear and nonlinear) with equally spaced nodes. All these problems have been considered in the literature [15,16,19,29, and the references therein]. For the linear and nonlinear cases, we compare the CCLS-method with the Finite Difference method (shortly FD) of second order, with equally spaced nodes [20] or with other methods discussed in the literature [19,15]. For the nonlinear case, we had used the Newton iterative method with tolerance 10−14 and maximum number of iterations 100. In all problems, we report the maximum error on the nodes. The presented results were obtained using Matlab.

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

7

Table 1 Numerical results of Example 1 by FD-method and CCLS-method. n

FD-method

CCLS-method

10 50 150 600 2000

4.579e−004 1.860e−005 2.068e−006 1.292e−007 1.163e−008

3.108e−004 3.064e−006 3.824e−007 2.824e−008 2.567e−009

Example 1 ([16]).  1 1−x   y+  y ′′ = (1 + x)2 (1 + x)2  1   y(0) = 1, y(1) = 2 1 with exact solution y(x) = 1+x . Using (9) and the support of Mathematica, for n Collocation Cubic Lidstone-Spline  1.0 − 5.144x + x 2 − 0.703x 3    2 3  1.703 − 4.51x + 0.792x − 0.356x p ν (x) = 2.228 − 3.996x + 0.604x 2 − 0.200x 3   2.625 − 3.569x + 0.462x 2 − 0.121x 3   2.936 − 3.212x + 0.357x 2 − 0.077x 3

= 5 we have the following analytical expression for the x x x x x

∈ [0, 0.2] ∈ [0.2, 0.4] ∈ [0.4, 0.6] ∈ [0.6, 0.8] ∈ [0.8, 1].

We compare the results obtained by the FD-method [20] with ours (see Table 1). Example 2 ([29,19]).   y ′′ + g = rC (y − T ) K kA  y(0) = y(a) = T   γa with exact solution y(x) = T + A 1 − cosh γ x − 1−cosh sinh γ a sinh γ x , where g = 71.7, a = 2.9, T = 325, k =  g 0.01, C/A = 4, r = 0.035, γ = rC k A and A = kγ 2 . Using (9) and the support of Mathematica, for n = 5 we have the following analytical expression for the Collocation Cubic Lidstone-Spline  325.0 + 933.758x − 3585.x 2 + 1973.98x 3 x ∈ [0, 0.58]    2 3  x ∈ [0.58, 1.16] 1483.36 − 1008.44x − 293.999x + 82.596x p ν (x) = 2401.6 − 1341.86x − 6.565x 2 x ∈ [1.16, 1.74]  2 3   3567.89 − 2058.65x + 424.586x − 82.596x x ∈ [1.74, 2.32]   26 103.2 − 31 800.6x + 13 588.6x 2 − 1973.98x 3 x ∈ [2.32, 2.9]. We compare the results obtained in [19] with ours. In Table 2, n is the number of subintervals and m is the degree of the approximating polynomial obtained from the Newton interpolating series (for further details see [19]). Remark 6. We observe that in the above examples the error for CCLS-method compares favourably with other methods considered and in some cases is much better.

8

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

Table 2 Numerical results of Example 2 by [19] and CCLS-method. m

[19]-method

n

CCLS-method

10 12 14 16

0.210e−000 0.755e−001 0.748e−002 0.127e−003

50 150 600 650

0.368e−002 0.389e−003 0.395e−005 0.370e−006

Table 3 Numerical results of Example 3 by CN-method [15] and CCLS-method. N

CN-method

n

CCLS-method

2 3 4 5 6 7

0.0277 5.7029e−004 1.1403e−004 1.0672e−006 5.6117e−007 2.0498e−009

3 5 10 30 99 200

2.363e−002 9.785e−003 2.623e−003 2.935e−004 2.644e−005 6.610e−006

Example 3 ([15]).  ′′ y = −(1 + 0.01y 2 ) + 0.01 cos3 (x) y(−1) = cos(−1), y(1) = cos(1) with exact solution y(x) = cos(x). Using (9) and the support of Mathematica, for n = 5 we have the following analytical expression for the Collocation Cubic Lidstone-Spline  −1.056 − 2.096x − 0.500x 2 + 0.0001x 3 x ∈ [−1, −0.6]    2 3  −0.431 − 2.383x − 0.4997x + 0.0005x x ∈ [−0.6, −0.2]  p ν (x) = 0.518 − 2.360x − 0.500x 2 x ∈ [−0.2, 0.2]  2 3  1.384 − 1.970x − 0.4997x − 0.0005x x ∈ [0.2, 0.6]   1.726 − 1.211x − 0.500x 2 − 0.0001x 3 x ∈ [0.6, 1]. We recall that, in Table 3, N is the degree of the polynomials and n the number of subintervals. For the CN-method, we considered polynomials of degree N = {2, 3, 4, 5, 6, 7} and 100 equally spaced nodes in [−1, 1] (fixed!). For the CCLS-method, we considered a number of equally spaced subintervals in [−1, 1] equal to n = {3, 5, 10, 30, 99, 200}. We observe that, for polynomials of degree almost equal, the spline function s (x) gives better results. 6. Conclusions We have presented a new Collocation Cubic Lidstone-Spline method based on cubic splines for solving BVPs of second order, in linear and nonlinear cases. It derives directly from piecewise Lidstone polynomials of degree 3 by requiring the continuity of the first derivative at the nodal points. Further developments are possible, for example, we can use splines of degree greater than 3. Furthermore, for odd order BVP, we can use Lidstone polynomials of II type (Complementary Lidstone) [14]. References [1] [2] [3] [4] [5]

R.P. Agarwal, P.J.Y. Wong, Lidstone polynomials and boundary value problems, Comput. Math. Appl. 17 (1989) 1397–1421. R.P. Agarwal, P.J.Y. Wong, Error Inequalities in Polynomial Interpolation and their Applications, Springer Science & Business Media, 2012. J.H. Ahlberg, T. Ito, A collocation method for two-point boundary value problems, Math. Comp. 29 (1975) 761–776. E.L. Albasiny, W.D. Hoskins, Cubic spline solutions to two-point boundary value problems, Comput. J. 12 (1969) 151–153. E.A. Al-Said, Spline solutions for system of second-order boundary value problems, Int. J. Comput. Math. 62 (1996) 143–154.

F.A. Costabile et al. / Mathematics and Computers in Simulation (

)

–

9

[6] E.A. Al-Said, Quadratic spline solutions of two point boundary value problem, J. Natur. Geom. 12 (1997) 125–134. [7] E.A. Al-Said, Cubic spline method for solving two-point boundary value problems, Korean J. Comput. Appl. Math. 5 (1998) 669–680. [8] U.M. Ascher, R.M.M. Mattheij, R.D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Vol. 13, Siam, 1994. [9] W.G. Bickley, Piecewise cubic interpolation and two-point boundary problems, Comput. J. 11 (1968) 206–208. [10] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Courier Corporation, 2001. [11] J. Chang, Q. Yang, L. Zhao, Comparison of B-spline method and finite difference method to solve BVP of linear ODEs, J. Comput. 6 (2011) 2149–2155. [12] C.C. Christara, K.S. Ng, Optimal quadratic and cubic spline collocation on nonuniform partitions, Computing 76 (2006) 227–257. [13] F. Costabile, Appunti di Calcolo Numerico con Software Didattico, Liguori, Napoli, 1992. [14] F.A. Costabile, F. Dell’Accio, R. Luceri, Explicit polynomial expansions of regular real functions by means of even order Bernoulli polynomials and boundary values, J. Comput. Appl. Math. 176 (2005) 77–90. [15] F.A. Costabile, A. Napoli, A method for polynomial approximation of the solution of general second order BVPs, Far East J. Appl. Math. 25 (2006) 289–305. [16] F. Costabile, A. Napoli, A collocation method for global approximation of general second order BVPs, Comput. Lett. 3 (2007) 23–34. [17] H.D. Doctor, A.B. Bulsari, N.L. Kalthia, Spline collocation approach to boundary value problems, Internat. J. Numer. Methods Fluids 4 (1984) 511–517. [18] D.J. Fyfe, The use of cubic splines in the solution of two-point boundary value problems, Comput. J. 12 (1969) 188–192. [19] G. Groza, M. Jianu, Polynomial approximations of solutions of boundary value problems for ODEs which arise from engineering, in: Proc. of RIGA 2014, Riemannian Geometry and Applications to Engineering and Economics, Bucharest, Romania. [20] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley & Sons, Inc., Interscience Publishers Inc., 1962. [21] M.K. Jain, T. Aziz, Spline function approximation for differential equations, Comput. Methods Appl. Mech. Engrg. 26 (1981) 129–143. [22] A.S.V.R. Kanth, Y.N. Reddy, Cubic spline for a class of singular two-point boundary value problems, Appl. Math. Comput. 170 (2005) 733–740. [23] S.A. Khuri, A. Sayfy, Spline collocation approach for the numerical solution of a generalized system of second-order boundary value problems, Appl. Math. Sci. 3 (2009) 2227–2239. [24] R. Lakshmi, M. Muthuselvi, Numerical solution for boundary value problem using difference method, Int. J. Innov. Res. Sci., Eng. Technol. 2 (2013). [25] G.J. Lidstone, Notes on the extension of Aitken’s theorem (for polynomial interpolation) to the Everett types, in: Proc. of the Edinburgh Mathematical Society (Series 2), Vol. 2, Cambridge Univ. Press, 1930, pp. 16–19. [26] C. Manni, F. Mazzia, A. Sestini, H. Speleers, BS2 methods for semi-linear second order boundary value problems, Appl. Math. Comput. 255 (2015) 147–156. [27] G. M¨ullenheim, Solving two-point boundary value problems with spline functions, IMA J. Numer. Anal. 12 (1992) 503–518. [28] D.N. Pop, A collocation method using cubic B-splines functions for solving second order linear value problems with conditions inside the interval [0, 1], Studia Universitatis Babes-Bolyai, Mathematica, 2010. [29] D.L. Powers, Boundary Value Problems and Partial Differential Equations, Academic Press, 2009. [30] M. Sakai, Numerical solution of boundary value problems for second order functional differential equations by the use of cubic splines, Mem. Fac. Sci., Kyushu Univ. Ser. A Math. 29 (1975) 113–122. [31] S.S. Siddiqi, E.H. Twizell, Spline solutions of linear twelfth-order boundary value problems, J. Comput. Appl. Math. 78 (1997) 371–390. [32] R.A. Usmani, M. Sakai, A connection between quartic spline solution and numerov solution of a boundary value problem, Int. J. Comput. Math. 26 (1989) 263–273. [33] R.A. Usmani, S.A. Warsi, Quintic spline solutions of boundary value problems, Comput. Math. Appl. 6 (1980) 197–203.