Pandey SpringerPlus (2016)5:326 DOI 10.1186/s40064-016-1969-z

Open Access

RESEARCH

Solving third‑order boundary value problems with quartic splines P. K. Pandey* *Correspondence: [email protected] Department of Mathematics, Dyal Singh College, University of Delhi, Lodhi Road, New Delhi 110003, India

Abstract  In this article, we present a novel second order numerical method for solving third order boundary value problems using the quartic polynomial splines. We establish the convergence of the method. We present numerical experiments to demonstrate the efficiency of the method and validity of our second order method, which shows that present method gives better results. Keywords:  Boundary-value problems, Finite-difference methods , Obstacle problems, Quartic polynomial splines Mathematics Subject Classification:  65L10, 65L12, 65L20

Background In this article we consider a quartic splines method for the numerical solution of the third order boundary value problems given as u′′′ (x) = f (x, u),

a ≤ x ≤ b,

(1)

subject to the boundary conditions

u(a) = α,

u′ (a) = β

and u′ (b) = γ

where α, β and γ are real constant. In environments and in most other areas of natural and applied sciences, the differential equations that govern the behavior of model systems are well-known. For instance, to describe the evolution of physical phenomena in fluctuating environments governed by third order differential equation (Ahmad et  al. 2012). The study of aero elasticity, sandwich beam analysis and beam deflection theory, electromagnetic waves, theory of thin film flow and incompressible flows and regularization of the Cauchy problem for one-dimensional hyperbolic conservation laws (Bressan 2000) are some other model systems in natural and applied sciences where the third order boundary value problems arise. The theoretical concepts of existence, uniqueness and convergence of the solution and some specific solution of problem (1) can be found in the literature (Howes 1982; Agarwal 1986; Gupta and Lakshmikantham 1991; Gregus 1987; Murty and Rao 1992; Feckan 1994; Henderson and Prasad 2001; Li 2010). The specific assumption to further ensure existence and uniqueness of the solution to problem (1) will not be considered. Thus the © 2016 Pandey. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Pandey SpringerPlus (2016)5:326

Page 2 of 10

existence and uniqueness of the solution to problem (1) is assumed. Further we assume that problem (1) is well pose. The emphasis in this article will be on the development of an efficient numerical method to deal with approximate numerical solution of the third order boundary value problem. The quality of a numerical method depends on the accuracy of the method to a great extent. Some efficient and accurate numerical methods for solving higher order boundary value problems are available in literature. Some researchers have studied and solved in particular third order boundary value problems with different boundary conditions using different methods for instance some literary work in Finite Difference Method (AlSaid 2001), Quintic Splines (Khan and Aziz 2003), Non polynomial spline method (Islam et  al. 2005, 2007; Srivastava and Kumar 2012), Quartic B-splines (Gao and Chi 2006), Haar wavelets method (Fazal-i-Haq and Ali 2011), Collocation quantic spline (Noor and Khalifa 1994), Reproducing Kernel Method (Li and Wu 2012) and references therein can be found. With advent of computers it gained important to develop more accurate numerical methods to solve higher order boundary value problems. Hence, the purpose of this article is to develop an efficient numerical method for solution of third order boundary value problems (1). We present our work in this article as follows. In the next section we derive a finite difference method. In section “Convergence analysis”, we discuss convergence of the proposed method under appropriate condition. The application of the proposed method on the test problems and illustrative numerical results so produced to show the efficiency in section “Numerical results”. Discussion and conclusion on the performance of the proposed method present in section “Conclusion”.

The difference method We define N finite numbers of nodal points of the domain [a, b], in which the solution of the problem (1) is desired, as a ≤ x0 < x1 < x2 < · · · < xN = b using uniform step length h such that xi = a + i.h, i = 0, 1, 2, . . . , N . Also, we let u(x) be the exact solution of (1) and we denote the numerical approximation of u(x) at node x = xi as ui. Let us denote fi as the approximation of the theoretical value of the source function f(x, u(x)) at node x = xi , i = 0, 1, 2, . . . , N . Thus the boundary value problem (1) at node x = xi may be written as u′′′ i = fi ,

(2)

a ≤ xi ≤ b,

subject to the boundary conditions

u0 = α,

u′0 = β

and

u′N = γ

We wish to determine the numerical approximation of the theoretical solution u(x) of the problem (1) at the nodal point xi , i = 1, 2, . . . , N . Define xi− 1 = xi − 21 h, i = 1, 2, . . . , N nodes in [a, b]. Let si− 1 be an approximation to 2 2 ui− 1 = u(xi− 1 ) solution of the problem (1) at these nodes obtained by the quartic spline 2 2 Si (x) passing through the points (xi− 1 , si− 1 ) and (xi+ 1 , si+ 1 ). For each ith segment, we 2 2 2 2 write polynomial quartic spline Si (x) in the form

Pandey SpringerPlus (2016)5:326

Page 3 of 10

 4  3  2 Si (x) = ci0 x − xi− 1 + ci1 x − xi− 1 + ci2 x − xi− 1 2 2 2   + ci3 x − xi− 1 + ci4

(3)

2

where ci0 , ci1 , ci2 , ci3 , ci4 are real finite constants. Then the quartic spline defined by

s(x) = Si (x),

i = 1, 2, . . . , N − 1 and

s(x) ∈ C 3 [a, b].

(4)

To determine constants ci0, ci1, ci2, ci3, ci4 we assume that Si (x) satisfies problem (1) with boundary conditions at xi− 1 and xi+ 1 . Following the idea in Fazal-i-Haq and Ali (2011) 2 2 we let

      ′ Si xi− 1 = si− 1 , Si xi+ 1 = si+ 1 , Si′ xi− 1 = si− 1, 2 2 2 2 2 2     ′′′ ′′ ′′ S xi− 1 = si− 1 and S xi− 1 = Ti− 1 , i = 1, 2, . . . , N − 1. 2

2

2

(5)

2

Thus we will obtain

  1  1 h Ti+ 1 − Ti− 1 , ci1 = Ti− 1 , ci2 = − Ti+ 1 + 3Ti− 1 , 2 2 2 2 2 24h 6 24  1 ci3 = s 1 − si− 1 , ci4 = si− 1 , i = 1, 2, . . . , N − 1. 2 2 h i+ 2

ci0 =

(6)

Using method of undetermined coefficients and Taylor’s series expansion, we discretize problem (2) at these nodes in [a, b],

3h3 T 1 + ti , i = 1 2 2 8 i− 2  h3  165Ti− 1 − 45Ti+ 1 + ti , i = 2 − 15si− 3 + 10si− 1 − 3si+ 1 = −8si−2 − 2 2 2 2 2 48   3 h Ti− 3 + Ti− 1 + ti , 3 ≤ i ≤ N − 1 si− 5 − 3si− 3 + 3si− 1 − si+ 1 = − 2 2 2 2 2 2 2   3 h −25Ti− 3 + 21Ti− 1 + ti , i = N si− 5 − 3si− 3 + 2si− 1 = hsi′ + 2 2 2 2 2 48 ′ 9si− 1 − si+ 1 = 8si−1 + 3hsi−1 −

(7)

where ti , i = 1, 2, . . . , N is truncation error. In discretization we have used boundary conditions in a natural way. After neglecting the ti in (7), at nodal points xi− 1 , i = 1, 2, . . . , N , we will obtain the 2 N × N linear or nonlinear system of equations depends on the source function f(x, u) in unknown si− 1. We have to solve a system of equations by an appropriate method. We 2 have applied either Gauss Seidel or Newton–Raphson iterative method to solve above system of Eq. (7) respectively for linear and nonlinear system of equations by following algorithms,

  S(k+1) = (L + d)−1 · −U · Sk + b ,

Pandey SpringerPlus (2016)5:326

Page 4 of 10

in which terms are defined as follow, 

S(k+1)

(k+1)

s1



  2  (k+1)    s3   2    =  ···      ···     (k+1)  s 1 N− 2

0



,

N ×1

0 10 −3 1 ··· ···

−3 0

··· ···

··· ···

3 −3 ··· ···

3 ··· ··· 1

0

0

−1 0

     U =  ···   ···  

9  −15   1   L+d =  ···   ···   

−1 0 ··· ···

−1 ··· ···

··· ···

0

··· ··· −3 1

··· ··· 3 −3



      ···   ···    2

N ×N



      ···   ···   −1  0

,

b = (bi )N ×1

N ×N

where

 3h3  8 Ti− 12 ,  8α + 3hβ − � �   3  5h  −8α − 11T 1 − 3T 1 , i− 2� i+ 2 � 16 bi = h3  + T , − T 3  2 i− 12 i−  �2 �    h3 hγ + 48 −25Ti− 3 + 21Ti− 1 , 2

and

2

  Fi s 1 , s 3 , . . . , sN − 1 ) = 0, 2

2

i=1 i=2 3≤i ≤N −1 i=N

i = 1, 2, . . . , N .

2

Let s and F are the all vector of values si− 1 and Fi respectively. 2





F(s + δs) = F(s) + Jδs + O (δs)2 . Neglecting the term O((δs)2 ) and assuming s + δs is root of F, then

sm+1 = sm + δs, where J−1 =

and

∂(f1 ,f2 ,...,fN ) ∂(s 1 ,s 3 ,...,sN − 1 ). 2

2

2

δs = −J−1 F,

m = 0, 1, 2 . . . .

Also some time we consider    δs = sm+1 − sm.

We compute numerical value of uN by using following second order approximation, 1 ′ uN = sN − 1 + hsN . (8) 2 2

Convergence analysis We will consider following linear test equation for convergence analysis of the proposed method (7). u′′′ (x) = f (x),

a ≤ x ≤ b.

subject to the boundary conditions u0 = α, posed method (7) in the matrix form as

(9) u′0 = β and u′N = γ . We can write the pro-

Pandey SpringerPlus (2016)5:326

Page 5 of 10

(10)

Du = a + t where



9  −15   1   D=  ···  ···   0

0

−1 10 −3 1 ··· ···

−3 3 −3 ··· ···

−1 3 ··· ··· 1

−1 ··· ··· −3 1

··· ··· 3 −3



     ···   ···   −1  2 N ×N

and u = (ui− 1 ), a = (ai ), and t = (ti ) are N-dimensional column vectors defined as, 2

 3h3   8α + 3hβ − 8 Ti− 21 ,    � �  3   11T 1 − 3T 1 ,  −8α − 5h i− 2 i+ 2 16 ai = � � 3  h   − 2 Ti− 3 + Ti− 1 ,  2 2   � �    hγ + h3 −25T 3 + 21T 1 , i− i− 48 2

ti =

 27h5 (5)  u 1, − 1920   i− 2     7h5 (5)  − u ,

i=2

1920 i− 2

i=N

8

i− 12

 o(h6 ),         31h5 u(5)1 ,

2

i=1 i=2 3≤i ≤N −1 i=N

i=1

3≤i ≤N −1

Let us solve test problem (9) by proposed method (7) after neglecting the terms ti. We will obtain a system of linear equations in si− 1. Solving the system of equations so 2 obtained by an iterative method, we get an approximate solution. Let us write system of equation in matrix form,

(11)

Ds = a

where s = (si− 1 ) is N-dimensional column vector of approximate solution of system of 2 equations obtained from (7). Let us define

ei− 1 = ui− 1 − si− 1 2

2

2

(12)

where si− 1 is an approximate value of ui− 1 , i = 1, 2, . . . , N . Thus from (10) and (11) we 2 2 can write an error equation

De = t

(13)

where e = (ei− 1 ), i = 1, 2, . . . , N is N-dimensional column vector. Let K = (kij ) be the 2 explicit inverse of nonsymmetric Toeplitz matrix D and defined as Jain et  al. (1987), Varga (2000), Horn and Johnson (1990), [24],

Pandey SpringerPlus (2016)5:326

Page 6 of 10

   

(2i−1)(4N −2i+1) , 24N 2 (2i − 1) c1 ,

 

(N +1−j)(2j+1)2 , 40N (8j−(2j−5)2 )



N +1−j 8N ,

c2 =

 

(N +1−j)(4N (j−1)+1)−8(j−1) , 24N

j=2

(N +1−j)(4N (j−1)+1)−8(j−1) , 8N

j+1≤i