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.