Cubic spline polynomial for non-linear singular two-point boundary value problems

Cubic spline polynomial for non-linear singular two-point boundary value problems

Applied Mathematics and Computation 189 (2007) 2017–2022 www.elsevier.com/locate/amc Cubic spline polynomial for non-linear singular two-point bounda...

134KB Sizes 1 Downloads 55 Views

Applied Mathematics and Computation 189 (2007) 2017–2022 www.elsevier.com/locate/amc

Cubic spline polynomial for non-linear singular two-point boundary value problems A.S.V. Ravi Kanth

*

Department of Mathematics, Vellore Institute of Technology, Vellore 632 014, Tamil Nadu, India

Abstract We present a cubic spline polynomial for the solution of non-linear singular two-point boundary value problems. The quesilinearization technique is used to reduce the non-linear problem to a sequence of linear problems. The resulting sets of differential equations are modified at the singular point and are treated by using cubic spline for finding the numerical solution. The method is tested on three physical model problems from the literature. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Two-point boundary value problems; Cubic spline; Quesilinearization; Singular point

1. Introduction We consider a class of non-linear singular two-point boundary value problem a y 00 ðxÞ þ y 0 ðxÞ ¼ f ðx; yÞ; x y 0 ð0Þ ¼ 0;

ð1Þ ð2Þ

yð1Þ ¼ b;

ð3Þ of oy

where a P 1 and b is a constant. We assume that f(x, y) is continuous function, and exists, is continuous and P 0. Due to the singularity at x ¼ 0 on the left side of the differential equation (1), direct numerical techniques face convergence difficulties. Attempts by many researchers for the removal of singularity are based on using the series expansion procedures in the neighborhood ð0; dÞ of singularity (d is vicinity of the singularity) and then solve the regular boundary value problem in the interval ðd; 1Þ using any numerical method. A threepoint finite difference scheme has been proposed by Russell and Shampine [6] with Newton’s iteration procedure for solving non-linear singular problems. A shooting algorithm based on a Taylor series method is developed by Rentrop [5]. The application of fourth order finite difference method to these problems has been discussed by Chawla et al. [2].

of oy

*

Present address: Department of Mathematics, Eritrea Institute of Technology, P.O. Box 106 99, Asmara, Eritrea. E-mail address: [email protected]

0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.01.002

2018

A.S.V. Ravi Kanth / Applied Mathematics and Computation 189 (2007) 2017–2022

In this paper, we discuss a direct method based on cubic spline approximation [1] for the solution of nonlinear singular two-point boundary value problems. An advantage of the method is that the coefficient matrix of the system is of the system of Hessenberg form. First we use the quesilinearization technique to reduce the given non-linear problem to a sequence of linear problems. The resulting sets of differential equations are modified at the singular point and are treated by using cubic spline for finding the numerical solution. The numerical method is tested for its efficiency by considering three physical model problems from the literature. 2. Method of solution The quesilinearization technique, has been used to reduce the given non-linear problem (1)–(3) to a sequence of linear problems. We choose a reasonable initial approximation for the function yðxÞ in f ðx; yÞ, call it as y(0)(x) and expand f ðx; yÞ around the function y(0)(x), we obtain    ð1Þ   ð0Þ   ð1Þ  of ð0Þ ¼ f x; y þ y y f x; y þ  ð6Þ oy ðx;y ð0Þ Þ or in general we can write for k ¼ 0; 1; 2; . . . (k = iteration index)         of f x; y ðkþ1Þ ¼ f x; y ðkÞ þ y ðkþ1Þ  y ðkÞ þ  oy ðx;y ðkÞ Þ

ð7Þ

Eq. (1) can be a ðkþ1Þ y xx ðxÞ þ y xðkþ1Þ ðxÞ þ aðkÞ ðxÞy ðkþ1Þ ðxÞ ¼ bðkÞ ðxÞ; x where aðkÞ ðxÞ ¼ 

  of ; oy ðx;y ðkÞ Þ

k ¼ 0; 1; 2; . . . ;

bðkÞ ðxÞ ¼ f ðx; y ðkÞ Þ  y ðkÞ



of oy

ð8Þ

 ðx;y ðkÞ Þ

subject to the boundary conditions y xðKþ1Þ ð0Þ ¼ 0; y

ðkþ1Þ

ð9Þ

ð1Þ ¼ b:

ð10Þ

To solve the set of linear singular boundary value problems given by (8) subject to the boundary conditions (9) and (10). We first modify Eq. (8) at the singular point x ¼ 0 and then apply the cubic spline method. Following [4] by L’Hospital rule we transform the given problem (8)–(10) into ðkþ1Þ y xx ðxÞ þ pðkÞ ðxÞy xðkþ1Þ ðxÞ þ qðkÞ ðxÞy ðkþ1Þ ðxÞ ¼ rðkÞ ðxÞ;

y xðKþ1Þ ð0Þ ðkþ1Þ

y

k ¼ 0; 1; 2; . . . ;

ð11Þ

¼ 0;

ð12Þ

ð1Þ ¼ b;

ð13Þ

where

( ðkÞ

p ðxÞ ¼

0; ðx ¼ 0Þ a ; x

ðx 6¼ 0Þ

( ;

ðkÞ

q ðxÞ ¼

aðkÞ ð0Þ ; aþ1 ðkÞ

ðx ¼ 0Þ

a ðxÞ; ðx 6¼ 0Þ

and ( ðkÞ

r ðxÞ ¼

bðkÞ ð0Þ ; aþ1 ðkÞ

ðx ¼ 0Þ

:

b ðxÞ; ðx 6¼ 0Þ

The range of the independent variable x is [0, 1]. We choose points 0 ¼ x0 ; x1 ; . . . ; xn ¼ 1. Starting at x0 a function y ðkþ1Þ ðxÞ (k = iteration index) is represented in ½x0 ; x1  by 2

3

S ðkþ1Þ ðxÞ ¼ a þ bðx  x0 Þ þ cðx  x0 Þ þ d 0 ðx  x0 Þ ;

k ¼ 0; 1; . . .

A.S.V. Ravi Kanth / Applied Mathematics and Computation 189 (2007) 2017–2022

2019

3

Proceeding to the next interval, ½x1 ; x2  we add a term d 1 ðx  x1 Þ ; proceeding into the next interval ½x2 ; x3 , we 3 add another term d 2 ðx  x2 Þ , and so on until at last we reach xn . So in general y ðkþ1Þ ðxÞ can be represented in the form S ðkþ1Þ ðxÞ ¼ a þ bðx  x0 Þ þ cðx  x0 Þ2 þ d 0 ðx  x0 Þ3 þ d 1 ðx  x1 Þ3 þ d 2 ðx  x2 Þ3 þ    þ d n1 ðx  xn1 Þ3 ;

x 2 ½x0 ; xn ;

ð14Þ

where a; b; c; d 0 ; d 1 ; d 2 ; . . . ; d n1 are unknowns to be found out. From this we get representation of the derivatives 2

2

2

S xðkþ1Þ ðxÞ ¼ b þ 2cðx  x0 Þ þ 3d 0 ðx  x0 Þ þ 3d 1 ðx  x1 Þ þ 3d 2 ðx  x2 Þ 2

þ    þ 3d n1 ðx  xn1 Þ ; ðkþ1Þ S xx ðxÞ

k ¼ 0; 1; . . . ;

ð15Þ

¼ 2c þ 6d 0 ðx  x0 Þ þ 6d 1 ðx  x1 Þ þ 6d 2 ðx  x2 Þ þ    þ 6d n1 ðx  xn1 Þ;

k ¼ 0; 1; . . .

ð16Þ

The cubic spline S ðkþ1Þ ðxÞ interpolates the function y ðkþ1Þ ðxÞ at the knots xi ¼ x0 þ ih ði ¼ 0; 1; . . . ; nÞ which is continuous together with the first and second derivatives in the interval [0, 1] and satisfies S ðkþ1Þ ðxi Þ ¼ y ðkþ1Þ ðxi Þ. If S ðkþ1Þ ðxÞ is a cubic polynomial on ½xi1 ; xi  then in general S ðkþ1Þ ðxÞ can be written as 2

3

3

3

S ðkþ1Þ ðxÞ ¼ a þ bðx  x0 Þ þ cðx  x0 Þ þ d 0 ðx  x0 Þ þ d 1 ðx  x1 Þ þ d 2 ðx  x2 Þ 3

þ    þ d i1 ðx  xi1 Þ ;

k ¼ 0; 1; . . .

ð17Þ

The spline approximation at the ith node for Eq. (11) is ðkÞ

ðkÞ

ðkÞ

ðkþ1Þ S xx ðxi Þ þ P i S ðkþ1Þ ðxi Þ þ Qi S ðkþ1Þ ðxi Þ ¼ Ri ; x

i ¼ 0; 1; . . . ; n  1; k ¼ 0; 1; . . .

ð18Þ

Substituting S ðkþ1Þ ðxÞ in the above equation, we get h i h i ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ aQi þ b Qi ðxi  x0 Þ þ P i þ c Qi ðxi  x0 Þ2 þ 2P i ðxi  x0 Þ þ 2 þ

i1 X

h i ðkÞ ðkÞ ðkÞ 3 2 d j Qi ðxi  xj Þ þ 3P i ðxi  xj Þ þ 6ðxi  xj Þ ¼ Ri ;

0 6 x 6 1;

ð19Þ

j¼0 ðkÞ

ðkÞ

ðkÞ

where P i ¼ P ðkÞ ðxi Þ, Qi ¼ QðkÞ ðxi Þ, and Ri ¼ RðkÞ ðxi Þ. We can see that the above system consists of ðn þ 3Þ coefficients to be evaluated. The boundary conditions (12) and (13) together with (19) provides us with sufficient number of equations for determination of these coefficients. Substituting S ðkþ1Þ ðxÞ in the boundary conditions (12) and (13), we obtain b¼0

ð20Þ

and 2

a þ bðxn  x0 Þ þ cðxn  x0 Þ þ

n1 X

3

d m ðxn  xm Þ ¼ b:

ð21Þ

m¼0

Using Eqs. (19)–(21), we get the system as Az ¼ B; ð22Þ  T ðkÞ ðkÞ T ðkÞ where B ¼ b; Rn ; Rn1 ; . . . ; R0 ; 0 , z ¼ ðd n1 ; d n2 ; . . . ; c; b; aÞ and A is an upper Heisenberg matrix given by 3 2 z1;nþ3 z1;1 z1;2    z1;nþ2 6 z2;1 z2;2    z2;nþ2 z2;nþ3 7 7 6 7 6 6 0 z3;2    z3;nþ3 z3;nþ3 7 7 6 ; 6 0 0 z4;3  z4;nþ3 7 7 6 7 6 .. 7 .. .. .. 6 .. 4 . . 5 . . . 0



0

znþ3;nþ2

znþ3;nþ3

2020

A.S.V. Ravi Kanth / Applied Mathematics and Computation 189 (2007) 2017–2022 3

2

where z1;m ¼ ðxn  xj Þ for j ¼ 0; 1; 2; . . . ; n  1 and m ¼ 1; 2; . . . ; n such that m  j ¼ 1 and z1;nþ1 ¼ ðxn  x0 Þ , z1;nþ2 ¼ ðxn  x0 Þ; z1;nþ3 ¼ 1. 3 2 ðkÞ zi;j ¼ QðkÞ i ¼ 2; 3; . . . ; n þ 2: j ¼ 1; 2; . . . ; n: m ¼ n þ 2  i and m ðxm  xl Þ þ 3P m ðxm  xl Þ þ 6ðxm  xl Þ, l = n j such that if m  l 6 0 then zi;j ¼ 0: ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ zm;nþ1 ¼ Qi ðxi  x0 Þ2 þ 2P i ðxi  x0 Þ þ 2, zm;nþ2 ¼ Qi ðxi  x0 Þ þ P i , and zm;nþ3 ¼ Qi for i ¼ n; n  1; . . . ; 0 and m ¼ 2; 3; . . . ; n þ 2 such that m þ i ¼ n þ 2. znþ3;j ¼ 0, for 1 6 j 6 n þ 2, znþ3;nþ2 ¼ 1 and znþ3;nþ3 ¼ 0: The matrix equation (22) consists of upper Heisenberg matrix A can be solved efficiently by Gauss elimination method, using back substitution. 3. Computational results and discussion In this section, we have described a computational method for solving non-linear singular two-point boundary value problems. It is a practical method and can easily be implemented on a computer to solve such problems. We have implemented the present method on three physical model examples: (i) Astronomy; (ii) Chemistry; (iii) Thermal explosions. The applicability of the results shows that the present method approximates the solution very well. Example 1. First, we consider the problem arising in Astronomy; the equilibrium of isothermal gas sphere can be described by 2 y 00 ðxÞ þ y 0 ðxÞ þ y 5 ðxÞ ¼ 0 x pffiffi with boundary conditions y 0 ð0Þ ¼ 0 and yð1Þ ¼ 23. This problem has earlier been considered by Russell and 1 . The computational Shampine [6], Hoog and Weiss [3] and Rentrop [5] and its exact solution is yðxÞ ¼ pffiffiffiffiffiffiffiffiffiffi 1þx2 =3 results are presented in Tables 1 and 2. Example 2. We consider the problem arising in Chemistry; the formulation of heat and mass transfer within porous catalyst particle leads to   cbð1yðxÞÞ 2 0 1þbð1yðxÞÞ 2 00 y ðxÞ þ y ðxÞ  / yðxÞe ¼0 x with boundary conditions y 0 ð0Þ ¼ 0 and yð1Þ ¼ 1. The computational results are presented in Table 3. Example 3. As a last example, we solve 1 y 00 ðxÞ þ y 0 ðxÞ þ eyðxÞ ¼ 0 x has earlier p been with boundary conditions y 0 ð0Þ ¼ 0 and yð1Þ ¼ 0. This  problem  ffiffiffi discussed by Russell and Shampine [6]. The exact solutions are yðxÞ ¼ 2 log BxBþ1 where B ¼ 3  2 2. The computational results 2 þ1 are presented in Tables 4 and 5. Table 1 Computational results for Example 1 x

Spline solution, h ¼ 1=50

Spline solution, h ¼ 1=100

Exact solution

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.99999776585428 0.99833533048663 0.99339731881833 0.98532764416099 0.97435344735938 0.96076805778639 0.94491067534081 0.92714531931739 0.90784126427522 0.88735654900981 0.86602540378444

0.99999944168089 0.99833694863610 0.99339878024765 0.98532886939410 0.97435438938421 0.96076870639181 0.94491105559580 0.92714548535674 0.90784129026725 0.88735651929041 0.86602540378444

1.00000000000000 0.99833748845958 0.99339926779878 0.98532927816429 0.97435470369244 0.96076892283052 0.94491118252306 0.92714554082311 0.90784129900320 0.88735650941611 0.86602540378443

A.S.V. Ravi Kanth / Applied Mathematics and Computation 189 (2007) 2017–2022

2021

Table 2 Computational results for Example 1 x

Spline solution

Exact solution

Finite difference

Collocation

Patch bases

0 0.125 0.250 0.375 0.500 0.625 0.75 0.875 1

0.99999863666149 0.99740466902369 0.98974221822877 0.97735472815483 0.96076839462410 0.94063390078531 0.91766286544420 0.89256962335054 0.86602540378444

1.000000000000000 0.997405961908059 0.989743318610787 0.977355554850441 0.960768922830522 0.940634162003544 0.917662935482247 0.892569603966722 0.866025403784438

0.99982 0.99731 0.98967 0.97730 0.96073 0.94060 0.91765 0.89256 0.86603

0.99994 0.99735 0.98969 0.97731 0.96074 0.94061 0.91765 0.89256 0.86603

0.99992 0.99737 0.98971 0.97733 0.96075 0.94062 0.91765 0.89257 0.86603

Table 3 Computational results for Example 2 x

Spline solution, h ¼ 1=50; a ¼ b ¼ c ¼ 1

Spline solution, h ¼ 1=50; a ¼ 6; b ¼ 0:05; c ¼ 20

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.83675985838429 0.83836491959750 0.84318428772800 0.85123026453741 0.86252249405263 0.87708663616863 0.89495243141739 0.91615107927542 0.94071183074544 0.96865767679753 1.00000000000000

0.00186369813659 0.00216698024782 0.00325730758512 0.00580362490765 0.01143983393732 0.02393441219970 0.05182774206917 0.11388182383972 0.24819925031990 0.51929683445694 1.00000000000000

Table 4 Computational results for Example 3 x

Spline solution, h ¼ 1=20

Spline solution, h ¼ 1=50

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.31666276936043 0.31323481598371 0.30298581783045 0.28601992615051 0.26250678586427 0.23267602814188 0.19681005729614 0.15523560923546 0.10831458301100 0.05643463769780 0

0.31668933643689 0.31326089850075 0.30301069507864 0.28604289758920 0.26252723764750 0.23269346654184 0.19682412855195 0.15524610889540 0.10832145570468 0.05643796864326 0

The computational results for Example 1 for different values of h are presented in Table 1. These solutions corresponds to the iteration index k ¼ 5 by taking y 0 ¼ 0 as the initial approximation. The iterations were ðkþ1Þ ðkÞ stopped when the absolute error criterion jy i  y i j 6 107 for all i is appeared. As it is apparent from Table 1, the computational results converges to the exact solution. The comparison of the computed solution of this example for h ¼ 641 with other methods is given in Table 2. It can be observed that our solutions compare well with the solutions obtained by other methods. It can be also be noticed that the solution near the singular point is not affected much due to the removal of singularity at x ¼ 0. The computational results for Example 2 are given for b ¼ c ¼ / ¼ 1 and b ¼ 0:05, c ¼ 20, / ¼ 6 in Table 3. These solutions corresponds to the iteration index k ¼ 5 by taking y 0 ¼ 0 as the initial approximation. The iterations were stopped

2022

A.S.V. Ravi Kanth / Applied Mathematics and Computation 189 (2007) 2017–2022

Table 5 Computational results for Example 3 x

Spline solution

Smaller solution

Finite differences

Patch bases

0 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000

0.31669129814469 0.31133689338372 0.29535914687153 0.26900997941952 0.23269475929685 0.18695217670046 0.13242901163954 0.06985212883356 0

0.31669 0.31134 0.29536 0.26901 0.23270 0.18695 0.13243 0.069853 0

0.31643 0.31135 0.29530 0.26897 0.23267 0.18693 0.13242 0.069847 0

0.31672 0.31135 0.29535 0.26902 0.23270 0.18696 0.13243 0.069854 0

ðkþ1Þ

ðkÞ

when the absolute error criterion jy i  y i j 6 107 for all i. The computational results for Example 3 for different values of h are given in Table 4. These results corresponds to the iteration index k ¼ 5 by taking y 0 ¼ 0 as the initial approximation. The iterations were stopped when the absolute error criterion ðkþ1Þ ðkÞ jy i  y i j 6 107 for all i. It can be observed pffiffiffi that the results agree well with the smaller solution, which is exact solution with the value of B ¼ 3  2 2. The comparison of the computed solution of this example for h ¼ 641 with other methods is given in Table 5. As it is evident from the computational results, the method gives O(h4) accuracy. The results obtained using this method are better than using the usual finite difference method with the same number of knots. In comparison with the finite difference methods, spline solution has its own advantages, once the solution has been computed; the information required for spline interpolation between mesh points is available. Whereas the finite difference methods only solves the problems at the chosen knots. References [1] W.G. Bickley, Piecewise cubic interpolation and two point boundary problems, Comput. J. 11 (1968) 206–208. [2] M.M. Chawla, R. Subramanina, H. Sathi, A fourth order method for a singular two point boundary value problems, BIT 28 (1988) 88– 97. [3] F.R. Hoog, R. Weiss, Difference methods for boundary value problems with singularity of the first kind, Math. Research Center MRC 1536, University of Wisconsin-Madison, 1975. [4] A.S.V. Ravi Kanth, Y.N. Reddy, Cubic spline for a class of singular two-point boundary value problems, Appl. Math. Comput. 170 (2005) 733–740. [5] P. Rentrop, A Taylor series method for the numerical solution of two point boundary value problems, Numer. Math. 31 (1979) 359– 375. [6] R.D. Russell, L.F. Shampine, Numerical methods for singular boundary value problems, SIAM J. Numer. Anal. 12 (1975) 13–36.