A class of Birkhoff–Lagrange-collocation methods for high order boundary value problems

Applied Numerical Mathematics 116 (2017) 129–140

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A class of Birkhoff–Lagrange-collocation methods for high order boundary value problems Francesco A. Costabile, Anna Napoli ∗ Department of Mathematics and Informatics, University of Calabria, 87036 Rende (CS), Italy

a r t i c l e

i n f o

a b s t r a c t

Article history: Available online 16 December 2016 Keywords: Boundary value problem Collocation methods Interpolation

A general procedure to determine collocation methods for high order boundary value problems is presented. These methods provide globally continuous differentiable solution in the form of polynomial functions and also numerical solution on a set of discrete points. Some special cases are considered. Numerical experiments are presented which support theoretical results and provide favorable comparisons with other existing methods. © 2016 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Let us consider the general r-th order ordinary differential equation (ODE)

y (r ) (x) = f (x, y(x)) , x ∈ I = [a, b] , r > 1   where y(x) = y (x), y  (x), . . . , y (q) (x) , 0 ≤ q < r. We associate with equation (1) the r boundary conditions B [ y ](x) = g (x) ,

g ∈ Rr , x ∈ ∂ I

(1)

(2)

where B is a linear operator. We suppose that f : [a, b] ×Rq+1 → R is continuous at least in the interior of the domain of interest and it satisfies a uniform Lipschitz condition in y, which means that there exists a nonnegative constant L such that, whenever (x, y 0 , y 1 , . . . , y q ) and (x, y 0 , y 1 , . . . , yq ) are in the domain of f , the inequality q          f x, y 0 , y 1 , . . . , yq − f x, y 0 , y 1 , . . . , yq  ≤ L  yk − yk 

(3)

k =0

holds. Moreover, we assume that the conditions for the existence and uniqueness of solution of problem (1)–(2) in a certain appropriate domain of [a, b] × Rq+1 are satisfied [1] and that the solution y (x) is differentiable with continuity up to what is necessary. Problems of these kinds model a wide spectrum of nonlinear phenomena. In fact, ordinary (and also partial) differential equations arise in the study of astrophysics, hydrodynamic and hydro magnetic stability, fluid dynamics, astronomy, beam and long wave theory, engineering and applied physics (see [7,21,23,30,31] and references therein).

*

Corresponding author. E-mail addresses: [email protected] (F.A. Costabile), [email protected] (A. Napoli).

http://dx.doi.org/10.1016/j.apnum.2016.12.003 0168-9274/© 2016 IMACS. Published by Elsevier B.V. All rights reserved.

130

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

Finite difference and finite element methods (see for instance [5,8,22]) have long histories as particularly flexible and powerful methods for the numerical solution of differential equations. In the last decades spectral methods (see [6,24,28] and references therein) have emerged as a valid alternative to those methods, being very much successful and extremely accurate for the numerical solution of ordinary or partial differential equations. They make use of global representations, usually by high-order polynomials or Fourier series. Spectral methods have been intensively and increasingly studied, especially, since the development of fast transform methods, with applications in problems where high accuracy is required. Theoretical studies and numerical experience have confirmed that for problems with smooth solutions they converge faster than finite difference and finite element methods and with a degree of accuracy that local methods cannot match. Summary of some applications is given in [24] (see also references therein). The basic idea of spectral methods is to use a set of basis functions φi , i = 0, . . . , N, also called trial or expansion approximating functions  N (very smooth and global such as polynomials), to represent the solution to a problem as a truncated series y N (x) = i =0 ai φi (x), x ∈ (a, b) where ai are the unknown coefficients. A spectral method is characterized by a specific way to determine these coefficients. For instance, in collocation methods (also called pseudo-spectral) the numerical approximation y N (x) of the solution y (x) to a problem is required to satisfy the given differential equation in a discrete set of points, called collocation nodes, of the interval (a, b). The choice of appropriate spectral representation depends on the kind of initial or boundary conditions involved in the problem. More recently, in [32] the authors, taking up an idea in [10–18,20], proved that the method, based on a Birkhoff type interpolation [15], is well conditioned. In this paper we will give a unified theory of collocation methods based on Birkhoff–Lagrange interpolation, already proposed by the authors and others (see [10,12–14,16–18,20,32]). The outline of the paper is the following: in Section 2 we derive the new class of methods and we give an a priori estimation of the global error. In Section 3, in order to implement the presented methods, we propose an algorithm to compute the numerical solution of (1)–(2) in a set of nodes. Then, in Section 4 we consider some special cases of the considered problem which include methods already present in the literature. The analysis is accompanied by numerical examples (Section 5), which support theoretical results. Comparison with other existing methods is also given. 2. The Birkhoff–Lagrange collocation methods Preliminary, we prove that problem (1)–(2) is equivalent to a Fredholm integral equation. We observe that (2) is a linear interpolation problem. Suppose that it has a unique solution, that is, there exists a unique polynomial of degree r − 1, P r −1 [ y ](x) ∈ Pr −1 which satisfies

B [P r −1 ] (x) = g ,

x ∈ ∂ I.

(4)

Moreover, if y ∈ C ( I ), then r

b K rx (x, t ) y (r ) (t )) dt

y (x) = P r −1 [ y ](x) +

(5)

a

where

K rx (x, t ) =

1

(r − 1)!

  (x − t )r+−1 − P r −1 (x − t )r+−1 (x)

(6)

is the Peano kernel of the remainder in the considered interpolation problem, and (x − t )+ = x − t if x ≥ t, (x − t )+ = 0 if x < t.





Remark 1. From (5), (2) and (4) it follows that B K rx (x) = 0, x ∈ ∂ I . Theorem 1. In the above hypothesis and notations the boundary value problem (1)–(2) is equivalent to the Fredholm integral equation

b K rx (x, t ) f (t , y(t )) dt

y (x) = P r −1 [ y ](x) +

(7)

a

where P r −1 [ y ](x) and K rx (x, t ) are, respectively, the unique polynomial solution and the corresponding Peano kernel in the interpolatory problem (2). Proof. Let y (x) be the solution of (1)–(2). Hence the linear interpolatory problem (2) has a unique solution, that is, there exists a unique polynomial P r −1 [ y ](x) of degree r − 1 and the related Peano kernel K rx (x, t ) such that (5) holds. The (7) follows tacking into account (1). Now, suppose that y (x) is solution of (7). Under the hypothesis on P r −1 [ y ] and K rx (x, t ), relations (2) and (5) hold. Therefore, from (5), (7) and the uniqueness of the Peano Kernel, the (1) follows.

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

131

Let {xi }m i =1 be m distinct points in [a, b ] and denote by l i (t ), i = 1, . . . , m, the fundamental Lagrange polynomials on the m ωm (t ) nodes xi , li (t ) = , where ωm (t ) = k=1 (t − xk ).  t − x ω ( x ) ( i) m i Theorem 2. If the solution y (x) of (1)–(2) is in C r +m ( I ), then

y (x) = P r −1 [ y ](x) +

m 

p r ,i (x) f (xi , y(xi )) + T r ,m ( y , x),

(8)

i =1

where

b K rx (x, t )li (t ) dt

p r ,i (x) =

i = 1, . . . , m ,

(9)

a

and the remainder term T r ,m ( y , x) is given by

b

1

T r ,m ( y , x) =

m!

K rx (x, t )ωm (t ) y (r +m) (ξx ) dt

a

being ξx a point in the smallest interval containing x and all xi , i = 1, . . . , m. Proof. From Lagrange interpolation

y (r ) (x) =

m 

li (x) y (r ) (xi ) + R m ( y , x)

(10)

i =1

where

1

R m ( y , x) =

m!

ωm (t ) y (r +m) (ξx )

(11)

is the remainder term. From (1), f (x, y (x)) = y (r ) (x). Then, by applying Theorem 1, if we insert (10) into (7), the (8) follows. Theorem 2 suggests to consider the implicitly defined polynomial

y r ,m (x) = P r −1 [ y ](x) +

m 



i =1



p r ,i (x) f xi , yr ,m (xi )

(12)



where yr ,m (xi ) = y r ,m (xi ), y r ,m (xi ), . . . , y r ,m (xi ) and p r ,i (x) is defined as in (9). (q)

Theorem 3 (The main Theorem). The polynomial y r ,m (x) of degree r + m − 1 defined in (12) satisfies the relations





B y r ,m (x) = g ,

x ∈ ∂I

  (r ) y r ,m (xi ) = f xi , yr ,m (xi )

(13) i = 1, . . . , m ,

(14)

that is, y r ,m (x) is a collocation polynomial for (1)–(2) at nodes xi .





Proof. It follows from (9) and Remark 1 that B p r ,i (x) = 0. This proves relation (13). In order to prove (14), by differentiating (12) r-times, we get (r )

(r )

y r ,m (x) = P r −1 [ y ](x) +

(r )

m 

(r )





p r ,k (x) f xk , yr ,m (xk ) .

k =1

(r )

From (9), p r ,k (x) = lk (x), k = 1, . . . , m, hence y r ,m (xi ) = fundamental Lagrange polynomials, (14) is proved.

m  k=1





lk (xi ) f xk , yr ,m (xk ) . From this and the known properties of

132

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

Remark 2. (Birkhoff-type interpolation). Theorem 3 is equivalent to the general Birkhoff interpolation problem [15]: given g ∈ Rr and αi ∈ R, i = 1, . . . , m, determine, if there exists, the polynomial Q (x) ∈ Pm+r −1 such that

B [ Q ](x) = g ,

x ∈ ∂ I,

Q (r ) (xi ) = αi ,

i = 1, . . . , m , x i ∈ I .

2.1. A-priori estimation of error In what follows, for all y ∈ C (q) ( I ) we define the norm

 y  = max

a≤t ≤b

q     (k)   y (t ) k =0

and the constants Q m =

m       pr ,i  and H = max  R m ( y , t ), where R m ( y , t ) is defined as in (11) and L is the Lipschitz i =1

constant of f . Further, we define F (x) =

b

a≤t ≤b

K rx (x, t )dt. a

Theorem 4. With the previous notations, if L Q m < 1, then

 y − yr ,m  ≤

HF 

(15)

.

1 − LQm

Proof. By deriving (8) and (12) s times, s = 0, . . . , q, we get (s)

y (s) (x) − y r ,m (x) =

m 





(s)



p r ,i (x) f (xi , y(xi )) − f xi , yr ,m (xi )

i =1

+

∂s ∂ xs

b K rx (x, t ) R m ( y , t ) dt . a

It follows that m     q          (s)  (s)    (k)  (s) (k)  y (x) − yr ,m (x) ≤  pr ,i (x) L  y (xi ) − yr ,m (xi ) + H  F (s) (x) i =1

k =0

m         (s)   ≤ L  y − yr ,m   pr ,i (x) + H  F (s) (x) . i =1

Hence m        y − yr ,m  ≤ L  y − yr ,m   p r ,i  + H  F  . i =1

This proves inequality (15). 3. Algorithms and implementation (s)

(s)

To calculate the approximate solution of problem (1)–(2) by (12) at x ∈ [a, b], we need the values yk = y r ,m (xk ), k = 1, . . . , m, s = 0, . . . , q. To do this we solve the following system (s)

(s)

yk = P r −1 [ y r ,m ](xk ) +

m  i =1

(s)

p r ,i (xk ) f (xi , yi )

s = 0, . . . , q, k = 1, . . . , m, where yi = y i , y i , . . . , y i

(q)

(16)

 .

We observe that, putting

⎛

A0 0 · · · 0

⎜ ⎜ 0 ... . .. ⎝ .. . 0 ··· 0

A=⎜ ⎜

⎞

.. ⎟ . ⎟ ⎟ ⎟ 0 ⎠

Aq

m(q+1)×m(q+1)

(17)

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

with

⎛

(s)

a˜ 1,1

⎜ . .. (s) a˜ m,1

As = ⎜ ⎝

⎞ (s) · · · a˜ 1,m .. ⎟ ⎟ . ⎠ (s) · · · a˜ m,m m×m

T Y = (Y 0 , . . . , Y q )m (q+1)×1 ,



(s)

(s)

s = 0, . . . , q ,

 ,

F m = ( f 1 , . . . , f m )T ,



q +1

(s)

Y s = y 1 , . . . , ym

F (Y ) = ( F m , . . . , F m )T ,



(s)

a˜ i , j = p r , j (xi ),

(s)

(q)



f i = f xi , y r ,m (xi ), . . . , y r ,m (xi ) ,



(s)

133



B s = P r −1 [ y r ,m ](x1 ), . . . , P r −1 [ y r ,m ](xm ) ,

C = B0, . . . , Bq

T m(q+1)×1

,

system (16) can be written in the form

Y − A F (Y ) = C .

(18)

For the existence and uniqueness of solution of (18) we have Theorem 5. If T = L  A ∞ < 1, system (18) has a unique solution which can be calculated by an iterative method

  (Y m )ν +1 = G (Y m )(ν ) ,

ν ≥0

(19)

with a fixed (Y m )0 ∈ Rm(q+1) and G (Y m ) = A F (Y m ) + C . Moreover, if Y is the exact solution of system (18), then

   ( Y m ) ν +1 − Y 

∞

≤

Tν 1−T

(Y m )1 − (Y m )0 ∞ .

Proof. In the above hypothesis and notations it’s easy to prove, with standard technique, that G is contractive if T < 1. Then, by applying the contraction mapping theorem [25], the result follows. Remark 3. Iterations (19) can be obtained by Picard’s iterations for the integral equation (7), that is,

b





K rx (x, t ) f t , (y(t ))ν dt .

( y (x))ν +1 = P r −1 [ y ](x) +

(20)

a

In fact, if we use Lagrange interpolation on the nodes xi , i = 1, . . . , m,





f t , (y(t ))ν =

m 





li (t ) f xi , (y(xi ))ν + R m ( f , t )

i =1

we have

( y (x))ν +1 = P r −1 [ y ](x) +

m  

f xi , (y(xi ))ν

i =1

= P r −1 [ y ](x) +

m 



b

b K rx (x, t )li (t ) dt

a



K rx (x, t ) R m ( f , t ) dt

+ a



p r ,i (x) f xi , (y(xi ))ν + S m r

(21)

i =1

with S m r =

b a

K rx (x, t ) R m ( f , t ) dt. For x = xk , k = 1, . . . , m, (21) coincides with (19).

Remark 4. If f is linear, that is, f (x, y(x)) =

q 

γ j (x) y ( j) (x) + g (x), then system (16) is a linear system:

j =0

BY =C,

where Y = (Y 1 , . . . , Y m ) T , Y i = y i , y i . . . , y i

(q)

(22)

 ,

134

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

⎛ ⎜

⎞ (k,s) (k,s) β1,1 · · · β1,q+1 ⎜ . .. ⎟ ⎟ B sk = ⎜ . ⎠ ⎝ .. (k,s) (k,s) βm,1 · · · βm,q+1 ⎛

⎞ · · · B 0m .. ⎟ .. . . ⎠ · · · B qm m(q+1)×m(q+1)

B 01

. B = ⎝ .. B q1

m×(q+1)

with s) βk(k,s,+ = pr(s,k) (xk )γs (xk ) − 1 1

and, for all (l, j ) = (k, s + 1), l = 1, . . . , m, j = 1, . . . , q + 1,

βl(,kj,s) = pr(s,k) (xl )γ j −1 (xk ) 

and C = B 0 , . . . , B q

T

with



(s) − P r −1 [ y 1 ](x1 ) −

Bs =

m  i =1

(s) (s) p r ,i (x1 ) g (xi ), . . . , − P r −1 [ ym ](xm ) −

m 

 (s) p r ,i (xm ) g (xi )

.

i =1

3.1. Numerical computation of the entries of matrix A (s)

(s)

To calculate the elements a˜ i ,k of the matrix A in (18), or the elements βi ,k of B in (22), we have to compute the integrals

b

ds (s) p r ,k (x) = s dx

K rx (x, t )lk (t ) dt

(23)

a

for x = xi . Integrating by parts, it remain to solve the problem of the computation of

x j

 

F i1 x j =

x j

 

F ik x j =

li (t ) dt , a

b

 

k = 2, . . . , n

F i ,k−1 (t )dt , a

M i1 x j =

b

 

M ik x j =

li (t ) dt , xj

M i ,k−1 (t )dt ,

k = 2, . . . , n

xj

i , j = 1, . . . m. To this aim it suffices to compute x j =tk tk−1





t1

··· c

c

rm,i (t ) dt dt 1 · · · dtk−1

(24)

c

where c = a or c = b, r0,0 (t ) = 1,

rm,i (t ) = (t − x1 ) · · · (t − xi −1 )(t − xi +1 ) · · · (t − xm ),

i = 1, 2, . . . , m .

For the computation of the integral (24), we use the recursive algorithm introduced in [11]: for each i = 1, . . . , m, let (i ) (i ) (i ) us consider the new points z j = x j if j < i, and z j = x j +1 if j ≥ i. Moreover, let us define g 0,1,c (x) = x − c and, for s = 1, . . . , m − 1, x=t j t j −1

(i ) g s, j ,c (x) =





··· c

c

t1

(i )

(i )

g s, j ,c (x) = x − z s Thus, if W i =

m 

F ik x j =

(i )

t − z2



 (i ) · · · t − zs dt dt 1 · · · dt j −1 .

(25)



(i )

(x − c ) j . For the computation of (25) the following recurrence formula [11] holds j! (i )

g s−1, j ,c (x) − jg s−1, j +1,c (x) .

(26)

(xi − xk ), then

k=1,k =i

 



c

We can easily compute g 0, j ,c (x) = (i )

(i )

t − z1

 

(i )

gm−1,k,a x j Wi

,

 

 

(i )

k

M ik x j = (−1)

gm−1,k,b x j Wi

.

(27)

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

135

Remark 5. An alternative approach for the exact computation of integrals (23) is to use a quadrature formula with a suitable degree of precision. Summarizing, the proposed method consists of the following steps: 1. determine the interpolation polynomial P r −1 [ y ](x) satisfying the boundary conditions, and compute the Peano remainder; 2. approximate y (r ) (x) by Lagrange interpolation on a set of fixed nodal point; 3. compute the elements of matrix A (17) and solve system (16); 4. obtain polynomial (12). 4. Special cases Now we consider some special cases of problem (1)-(2), that is, we associate with equation (1) particular kinds of boundary conditions and for each problem we determine the polynomials P r −1 [ y ](x) and K rx (x, t ). We also show how the proposed class of methods include the methods presented in previous works [10,12–14,17,18,20,32]. 4.1. Case r = 2 For the exact solution y (x) of the second order BVP





y  (x) = f x, y (x) , y  (x) ,

y (−1) = y 0 , y (1) = y 1

(28)

x ∈ [−1, 1], it is known that

y (x) =

y1 + y0 2

where

K 2x (x, t ) =

+x

1

y1 − y0

+

2





K 2x (x, t ) f x, y (x) , y  (x) dt

−1

⎧ (t + 1)(x − 1) ⎪ ⎨

t≤x

⎪ ⎩ (x + 1)(t − 1)

x
2 2

By applying method (12) we get [16]

y 2,m (x) =

y1 + y0 2

+x

y1 − y0 2

+

m 



p 2,i (x) f xi , y (xi ) , y  (xi )



i =1

1

with p 2,i (x) = −1 K 2 (x, t )li (t )dt . If xi = cos mπ+i1 , i = 1, . . . , m, we obtain [10]

p r ,i (x) =

1 m+1

sin

πi

"

m+1

m  G k (x)

where

G k (x) =

T k+1 (x) k+1

−

k

k =2

T k−1 (x) k−1

⎧ ⎪ ⎪ ⎨

sin

2x

2 + k −1 ⎪ 2 ⎪

⎩

k2 − 1

kπ i m+1



+ x2 − 1 sin

πi

#

m+1

even k odd k

and T k−1 (x), T k+1 (x) are the Chebyshev polynomials of the first kind and degree k − 1 and k + 1 respectively. The same method has been presented in [32] (without citations of [10] and [16]), where also stability has been studied. 4.2. Case r = 2n, n > 1 Assume [a, b] = [0, 1]. Several types of boundary conditions can be considered. – Hermite boundary conditions [17]:

y (h) (0) = αh , with

y (h) (1) = βh ,

h = 0, . . . , n − 1

αh , βh , h = 0, . . . , n − 1 real constants.

136

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

In this case P r −1 [ y ] is the Hermite interpolation polynomial of degree 2n − 1 [3]

P 2n−1 [ y ](x) =

n −1 



y (i ) (0) H i1 (x) + y (i ) (1) H i2 (x)

i =0

with

$

H i1 (x) =

H i2 (x) =

%

n − i −1 xi (1 − x)n  n + s − 1

i!

s =0

x (1 − x) n

i −1 $ i n− 

i!

The kernel is

K 2n (x, t ) =

n−1

%

n+s−1 n−1

s =0

⎧ n −1  (−t )2n−i −1 ⎪ ⎪ ⎪ H i1 (x) ⎪ ⎨ (2n − i − 1)! i =0 n −1 

⎪ ⎪ ⎪ ⎪ ⎩−

i =0

xs

(1 − x)s .

0≤t ≤x

(1 − t )2n−i −1 H i2 (x) x ≤ t ≤ 1 . (2n − i − 1)!

Hermite boundary conditions have been also considered by Agarwal [1], but for the numerical solution of the problem he uses Picard iterations and quasilinearization method. – Lidstone boundary conditions [12]:

y (2h) (0) = αh ,

y (2h) (1) = βh ,

h = 0, . . . , n − 1

(29)

where αh , βh , h = 0, . . . , n are real constants. In this case P r −1 [ y ] is the Lidstone interpolating polynomial [4] of degree 2n − 1

P 2n−1 [ y ](x) =

n −1  



y (2k) (0)k (1 − x) + y (2k) (1)k (x)

k =0

where k (x) are the Lidstone polynomials of degree 2k − 1 [4] and the function K 2n (x, t ) is

K 2n (x, t )

⎧ n −1 ⎪  ⎪ t 2n−2k−1 ⎪ ⎪ k (1 − x) t ≤ x ⎪ ⎨ (2n − 2k − 1)! k =0

n −1 ⎪  (1 − t )2n−2k−1 ⎪ ⎪ ⎪ k (x) ⎪ ⎩ (2n − 2k − 1)!

x≤t.

k =0

Also Lidstone boundary conditions have been considered by Agarwal [1], but he doesn’t construct a collocation method. 4.3. Case r = 2n + 1 If we consider the following boundary conditions

y (0) = α0 ,

y (2h−1) (0) = αh ,

y (2h−1) (1) = βh ,

h = 1, . . . , n

(30)

with α0 , αh , βh , h = 1, . . . , n real constants, then P r −1 [ y ] is the complementary Lidstone interpolating polynomial [9] of degree 2n

P 2n [ y ](x) = y (0) +

n  



y (2k−1) (0) ( v k (1) − v k (1 − x)) + y (2k−1) (1) ( v k (x) − v k (0)) ,

k =1

where v k (x) are the complementary Lidstone polynomials of degree k [2,9]. The kernel is

K 2n−1 (x, t ) =

⎧ n  ⎪ t 2n t 2n−2k+1 ⎪ ⎪ + ( v k (1 − x) − v k (1)) t ≤ x ⎪ ⎪ (2n − 2k + 1)! ⎨ (2n)! k =1

n ⎪  ⎪ (1 − t )2n−2k+1 ⎪ ⎪ − ( v k (x) − v k (0)) ⎪ ⎩ (2n − 2k + 1)!

x≤t.

k =1

In [12] the proposed method applied to equation (1) with conditions (29) and (30) respectively, is examined in detail.

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

137

4.4. Other special boundary conditions If r = n − 1, n > 2, and [a, b] = [0, 1], we can consider Bernoulli BCs [19]

y (0) = β0 ,

y (k) (1) − y (k) (0) = βk+1 ,

y (1) = β1 ,

k = 1, . . . , n − 2

with βk , k = 0, . . . , n − 1 real constants. The method has been examined in [13]. If r = n + 1 and [a, b] = [−1, 1], different Birkhoff -type boundary conditions can be examined, for example:

–

y (−1) = ω0 ,

y (1) = ωn ,

y  (xi ) = ωi ,

i = 1, . . . , n − 1

(31)

where ωi , i = 0, . . . , n are real constants and xi , i = 1, . . . , n − 1, are n − 1 distinct points in (−1, 1) with −1 < x1 < · · · < xn−1 < 1 [16];

–

y (k) (−1) = αk ,

y ( s ) ( xi ) = ωi

k = 0, . . . , s − 1 ,

i = 1, . . . , n − s + 1

(32)

where s ≤ n; αk , k = 0, . . . , s − 1 and ωi , i = 1, . . . , n − s + 1, are real constants; xi , i = 1, . . . , n − s + 1, are distinct points in (−1, 1) with −1 < x1 < · · · < xn−s+1 < 1. In [14] and [18] there is the description of the proposed method applied to equation (1) with conditions respectively (31) and (32). 5. Numerical examples Now we present some numerical results obtained by applying method (12), which we call BLC method, to find numerical approximations y r ,m (x) to the solution of some test problems. We compare the results with the ones obtained with other existing methods. In order to solve the nonlinear system (16), we use the so-called modified Newton method [29] (the same Jacobian matrix is used for more than one iteration) and we use algorithm (26) for the computation of the entries of the matrix, when the polynomials pr ,i (x) are not explicitly known. Since the true solutions are known, we consider the error function em (x) =  y (x) − y r ,m (x). If not specified, the collocation points are chosen to be equidistant points in the interval [a, b]. In the considered cases analogous results are obtained using as nodes the Chebyshev–Gauss–Lobatto points and the Legendre–Gauss–Lobatto points. The maximum values of em (x) over the interval [a, b] have also been calculated by using Matlab, particularly the following solvers:

• the powerful tool Chebfun; • the built-in solver bvp4c (with optional parameters RelTol=AbsTol= 1e −17 and a mesh of 200 points), a finite difference code that implements the three-stage Lobatto IIIa formula;

• the built-in solver bvp5c (with RelTol=AbsTol= 1e −17), a finite difference code that implements the four-stage Lobatto IIIa formula. Example 1. We consider the linear ninth-order BVP [34]

⎧ (9 ) ⎨ y (x) = −9e x + y (x) x ∈ [0, 1] y ( j ) (0) = 1 − j j = 0, . . . , 4 ⎩ ( j) j = 0, . . . , 3. y (1) = − j e

(33)

The exact solution is y (x) = (1 − x)e x . ( j) The unique polynomial P 8 (x) = P 8 [ y ](x) of degree 8 which satisfies the boundary conditions P 8 (0) = 1 − j for j = ( j)

0, . . . , 4, and P 8 (1) = − j e j = 0, . . . , 3 is

P 8 (x) = 1 −

1 2

$ + From (6) we get

x2 −

1321 12

1 3

−

x3 − 81 2

1 8

$ x4 +

%

1 e x6 +

31 2

$

1e −

71 2

e−

253 6 193 2

% x5

%

$ x7 +

685 24

−

21 2

% e x8 .

138

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

Table 1 Absolute error em (x) in MDM and BLC methods for problem (33). x

MDM [34]

BLC m=4

BLC m=6

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

2.0e−10 2.0e−10 2.0e−10 2.0e−10 2.0e−10 6.0e−10 1.0e−9 2.0e−9 3.4e−9

1.45e−14 3.93e−13 2.16e−12 5.70e−12 9.27e−12 1.00e−11 7.04e−12 2.70e−12 2.98e−13

0.00 1.11e−16 9.99e−15 2.00e−15 2.55e−15 2.66e−15 2.44e−15 2.83e−15 4.91e−15

Table 2 Maximum absolute error in problem (33) using Matlab build-in functions.

Chebfun 1.46

bvp4c 1.55e−12

bvp5c 4.44e−16

⎧ 4 4 4x5 + 6x6 − 4x7 + x8 )+ ⎪ ⎪ 70t 5 (x − ⎪ 3 ⎪ 56t (− x + 10x5 − 20x6 + 15x7 − 4x8 )+ ⎪ ⎪ ⎪ 6 2 5 ⎪ 28t ( x − 20x + 45x6 − 36x7 + 10x8 )+ ⎪ ⎪ ⎪ 7 5 ⎪ 8t (−x + 35x − 84x6 + 70x7 − 20x8 )+ ⎪ ⎪ ⎪ ⎪ 0≤t ≤x ⎪ t 8 (1 − 56x5 + 140x6 − 120x7 + 35x8 ) 1 ⎨ K 9 (x, t ) = · −x8 + 8tx7 − 28t 2 x6 + 56t 3 x5 + 8! ⎪ ⎪ ⎪ 70t 4 (−4x5 + 6x6 − 4x7 + x8 )+ ⎪ ⎪ ⎪ ⎪ 56t 5 (10x5 − 20x6 + 15x7 − 4x8 )+ ⎪ ⎪ ⎪ ⎪ ⎪ 28t 6 (−20x5 + 45x6 − 36x7 + 10x8 )+ ⎪ ⎪ ⎪ ⎪ 8t 7 (35x5 − 84x6 + 70x7 − 20x8 )+ ⎪ ⎩ 8 t (−56x5 + 140x6 − 120x7 + 35x8 ) x ≤ t ≤ 1. Now, using (27), we obtain the values of the integrals (23). Hence we can solve system (16) and, by substituting in (12), we obtain the approximating solution to problem (33). Table 1 shows the numerical results. The absolute errors computed are compared with those obtained in [34], where a modified decomposition method (MDM) is applied for the solution of problem (33). The second and third columns of the table show the error respectively in the MDM and BLC methods, using in both cases polynomials of degree 12. The last column contains the error in the approximation by a polynomial of degree 14 using BLC method. The maximum absolute error has also been calculated by using Matlab. The results are displayed in Table 2. Example 2 (Hermite boundary conditions). Consider now the following nonlinear problem [27]

⎧   ⎨ y (i v ) (x) = sin x + sin2 x − y  (x) 2 y (0) = 0 y  (0) = 1 ⎩ y (1) = sin(1) y  (1) = cos(1).

x ∈ [0, 1] (34)

The solution for this problem is y (x) = sin(x). Problems of this kind furnish a model to study, for example, traveling waves in suspension bridges [26] or the bending of an elastic beam [35]. Suspension bridges are generally susceptible to visible oscillations, which, if not controlled, can lead to failure of the bridge. In this case f represents the forcing term acting on the bridge (including the force due to the cables which are considered as a spring with a one-sided restoring, the gravitation force and the external force due to the wind or other external sources); the solution y represents the vertical displacement when the bridge is bending. In the case of elastic beam f represents the force exerted on the beam by the supports. x measures the position along the beam (x = 0 is the left-hand endpoint of the beam), y is the height of the beam, y  indicates the slope of the beam at x, y  measures the curvature of the graph of y (x). In physical terms, y  measures the bending moment of the beam at x, that is the torque that the load places on the beam at x. The considered boundary conditions state that the beam has both endpoints simply supported, that is, it is placed upon two supporting structures. Moreover, the beam at the wall is not horizontal, in fact the derivative of the deflection function are not zero at that points. Table 3 shows the comparison between the NMD method [27] and the BLC method with m = 5

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

139

Table 3 Error of NMD and BLC methods – problem (34). t

NMD [27]

BLC m=5

BLC m=9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

7.78e−8 2.72e−7 5.24e−7 7.77e−7 9.71e−7 1.05e−6 9.63e−7 6.84e−7 2.71e−7

4.45e−10 5.54e−10 8.95e−11 2.03e−10 3.32e−11 1.53e−10 9.48e−11 5.18e−10 4.15e−10

1.53e−15 3.02e−15 7.77e−16 6.66e−16 5.55e−17 0 0 1.11e−16 0

Table 4 Maximum absolute error in problem (34) using Matlab build-in functions.

Chebfun 1.67e−16

bvp4c 1.22e−8

bvp5c 8.88e−16

Table 5 Maximum absolute error of problem (35). m

BLC

h

Method in [33]

5 7 8 9 10 12

1.97e−4 1.65e−6 8.16e−7 1.17e−8 6.45e−9 3.67e−11

1/10 1/20 1/40 1/80 1/160 1/320

1.22e−4 7.62e−6 4.76e−7 2.97e−8 1.86e−9 1.16e−10

Table 6 Maximum absolute error in problem (35) using Matlab build-in functions.

Chebfun 8.14e−7

bvp4c

bvp5c

8.54e−7

8.13e−7

and m = 9 respectively. The approximating polynomial of NMD method has degree 11, while the polynomial considered in BLC method for m = 5 has degree 8. The maximum absolute errors calculated by using Matlab are displayed in Table 4. Example 3 (Lidstone boundary conditions). Consider the nonlinear sixth-order boundary value problem [33]

⎧ ⎪ ⎨ ⎪ ⎩

where

$ − y ( vi )

p

= σ1 (x) y + σ2 (x)

y (2i ) (0) = y (2i ) (1) = 0

y 

π2



%r + σ3 (x)

y (i v )

π4

s + q(x)

x ∈ [0, 1]

(35)

i = 0, 1 , 2

σi (x) are continuous functions,

 q(x) = sin(π x) π 6 − σ1 (x) sin p −1 (π x) + (−1)r +1 σ2 (x) sinr −1 (π x) − σ3 (x) sins−1 (π x)

and p , r , s are some positive integers ≥ 1. Then y (x) = sin(π x) is the solution of (35). Let p = 1, r = 2, s = 3 and σi (x) = 0.1 cos(π x), i = 1, 2, 3. In Table 5 we compare the results obtained by method (12) for several values of m, with the results obtained by the method in [33] for different values of the step h. To demonstrate the accuracy of our method we show the maximum absolute error in the two cases. Further, method (12) provides an explicit expression of the approximate solution. The maximum errors (in absolute value) by the Matlab solvers are in Table 6. 6. Conclusions In this paper we have presented a general procedure to determine collocation methods for the numerical solution of boundary value problems of r-th order. For two positive integers r , m a polynomial of degree r + m − 1 approximating the

140

F.A. Costabile, A. Napoli / Applied Numerical Mathematics 116 (2017) 129–140

exact solution is given explicitly. It is a collocation polynomial for the considered boundary value problem. The approximate solution contains the unique interpolating polynomial which satisfies the boundary conditions. The approximation of f does not depend on the nature of the function under consideration. In fact, we use Lagrange polynomials in an arbitrary set of nodes to approximate f . Moreover the determined unknowns provide the approximate solution in the nodal points. An a-priori estimation of error is given. We showed that the presented procedure is general and includes many special cases. Numerical experiments and comparisons with other existing methods are presented and show that the proposed method is highly competitive with other methods already present in the literature. Our future direction of work is to study the rate of convergence and stability properties of the method, and the extension to multidimensional case. Acknowledgements The authors would like to thank one of the anonymous reviewers for the helpful and constructive comments that contributed to improve the final version of the paper. This research was supported by INDAM-GNCS project 2016. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

R.P. Agarwal, Boundary Value Problems for Higher Order Differential Equations, World Scientific Publishing Co., Inc., Teaneck, NJ, 1986. R.P. Agarwal, S. Pinelas, P.J.Y. Wong, Complementary Lidstone interpolation and boundary value problems, J. Inequal. Appl. 2009 (2009) 624631. R.P. Agarwal, P.J. Wong, Error Inequalities in Polynomial Interpolation and Their Applications, vol. 262, Springer Science & Business Media, 2012. R.P. Agarwal, P.J.Y. Wong, Lidstone polynomials and boundary value problems, Comput. Math. Appl. 17 (1989) 1397–1421. S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, third edition, Texts in Applied Mathematics, vol. 15, Springer, New York, 2008. C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Technical Report, Springer Series in Computational Physics, Springer-Verlag, New York, 1988. S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, International Series of Monographs on Physics, Clarendon, Oxford, 1961. M. Chawla, C. Katti, Finite difference methods for two-point boundary value problems involving high order differential equations, BIT Numer. Math. 19 (1979) 27–33. F. 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. F. 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. F. Costabile, A. Napoli, A class of collocation methods for numerical integration of initial value problems, Comput. Math. Appl. 62 (2011) 3221–3235. F. Costabile, A. Napoli, Collocation for high-order differential equations with Lidstone boundary conditions, J. Appl. Math. 2012 (2012) 120792. F. Costabile, A. Napoli, Numerical solution of high order Bernoulli boundary value problems, J. Appl. Math. 2014 (2014) 276585. F. Costabile, A. Napoli, A method for high-order multipoint boundary value problems with Birkhoff-type conditions, Int. J. Comput. Math. 92 (2015) 192–200. F.A. Costabile, E. Longo, A Birkhoff interpolation problem and application, Calcolo 47 (2010) 49–63. F.A. Costabile, E. Longo, A new collocation method for a BVP, in: Applied and Industrial Mathematics in Italy III, in: Ser. Adv. Math. Appl. Sci., vol. 82, World Sci. Publ., Hackensack, NJ, 2010, pp. 289–297. F.A. Costabile, A. Napoli, Collocation for high order differential equations with two-points Hermite boundary conditions, Appl. Numer. Math. 87 (2015) 157–167. F.A. Costabile, A. Napoli, A multipoint Birkhoff type boundary value problem, J. Numer. Math. 23 (2015) 1–11. F.A. Costabile, A. Serpe, A. Bruzio, No classic boundary conditions, in: Proceedings of World Congress on Engineering 2007, London, July 2–4, 2007, pp. 918–921. M. Dehghan, S. Aryanmehr, M. Eslahchi, A technique for the numerical solution of initial-value problems based on a class of Birkhoff-type interpolation method, J. Comput. Appl. Math. 244 (2013) 125–139. M. Dehghan, F. Shakeri, Solution of an integro-differential equation arising in oscillating magnetic fields using He’s homotopy perturbation method, Prog. Electromagn. Res. 78 (2008) 361–376. E.J. Doedel, Finite difference collocation methods for nonlinear two point boundary value problems, SIAM J. Numer. Anal. 16 (1979) 173–185. M. Eslahchi, M. Dehghan, S. Ahmadi_Asl, The general Jacobi matrix method for solving some nonlinear ordinary differential equations, Appl. Math. Model. 36 (2012) 3387–3398. D. Gottlieb, S.A. Orszag, Numerical analysis of spectral methods: theory and applications, in: CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 26, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1977. A. Granas, J. Dugundji, Fixed Point Theory, Springer Science & Business Media, 2013. A. Lazer, P. McKenna, Large-amplitude periodic oscillations in suspension bridges: some new connections with nonlinear analysis, SIAM Rev. 32 (1990) 537–578. M.A. Noor, S.T. Mohyud-Din, An efficient method for fourth-order boundary value problems, Comput. Math. Appl. 54 (2007) 1101–1111. B. Ogundare, On the pseudo-spectral method of solving linear ordinary differential equations, J. Math. Stat. 5 (2009) 136. A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, second edition, Texts in Applied Mathematics, vol. 37, Springer-Verlag, Berlin, 2007. A. Saadatmandi, M. Dehghan, Variational iteration method for solving a generalized pantograph equation, Comput. Math. Appl. 58 (2009) 2190–2196. F. Shakeri, M. Dehghan, Solution of delay differential equations via a homotopy perturbation method, Math. Comput. Model. 48 (2008) 486–498. L.L. Wang, M.D. Samson, X. Zhao, A well-conditioned collocation method using a pseudospectral integration matrix, SIAM J. Sci. Comput. 36 (2014) 907–929. Y.M. Wang, H.Y. Jiang, R.P. Agarwal, A fourth-order compact finite difference method for higher-order Lidstone boundary value problems, Comput. Math. Appl. 56 (2008) 499–521. A.M. Wazwaz, Approximate solutions to boundary value problems of higher order by the modified decomposition method, Comput. Math. Appl. 40 (2000) 679–691. B. Yang, Estimates of positive solutions to a boundary value problem for the beam equation, Commun. Math. Anal. 2 (2007) 13–21.