Applied Mathematics and Computation 154 (2004) 175–182 www.elsevier.com/locate/amc
Parametric cubic spline solution of two point boundary value problems Arshad Khan Department of Applied Mathematics, Faculty of Engineering and Technology, Aligarh Muslim University, Aligarh 202-002, India
Abstract In this paper, we use parametric cubic spline function to develop a numerical method, which is fourth order for a specific choice of the parameter, for computing smooth approximations to the solution for second order boundary value problems. Some numerical evidence is also included to demonstrate the superiority of our method. 2003 Elsevier Inc. All rights reserved. Keywords: Parametric cubic spline function; Finite difference method; Boundary value problems; NumerovÕs method
1. Introduction We consider the two-point boundary value problem y 00 ðxÞ ¼ f ðxÞyðxÞ þ gðxÞ; yðaÞ ¼ a0 ; yðbÞ ¼ a1 ;
a 6 x 6 b;
ð1:1Þ
where f ðxÞ and gðxÞ are continuous functions on ½a; b and a; b; a0 ; a1 are arbitrary real finite constants. Such problems arise in the theory which describes the deflection of plates and a number of other scientific applications. In general it is difficult to obtain the analytical solution of (1.1) for arbitrary choices of f ðxÞ and gðxÞ. We usually employ some numerical method for obtaining an approximate solution of the problem (1.1). The standard numerical methods
E-mail addresses:
[email protected],
[email protected] (A. Khan). 0096-3003/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(03)00701-X
176
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
for the numerical treatment of (1.1) consist of finite difference methods discussed by many authors, see, for example [7,11]. The subject of obtaining spline solutions for the initial as well as boundary value problems is discussed by Ahlberg et al. [1]. Since then Albasiny and Hoskins [2], Bickley [8], Fyfe [10] and Sakai and Usmani [15] have used the cubic spline for obtaining approximations and Al-Said [4] has used cubic spline for the solution of (1.1). Al-Said [5] has used quadratic spline for obtaining smooth approximations for the solution and its first derivative, of second order obstacle problems and of boundary value problems of the form (1.1). The use of higher order spline functions and collocation methods with splines as basis functions for solving various second order boundary value problems were demonstrated by different authors, see, for example [1,6,12–16,18] and references therein. In this paper, we have derived a uniformly convergent uniform mesh difference scheme using parametric cubic spline for the solution of (1.1). The main idea is to use the condition of continuity as a discretization equation for (1.1). In comparison with the finite difference methods, spline solution has its own advantages. For example, once the solution has been computed the information required for spline interpolation between mesh points is available. This is particularly important when the solution of boundary value problem is required at various locations in the interval ½a; b. Analysis of the method shows that it is second order convergent for arbitrary a; b such that a þ b ¼ 1=2 and fourth order convergent for a ¼ 1=12, b ¼ 5=12. Our second order method outperforms other second order methods and fourth order method reduces to well-known NumerovÕs method and fourth order quartic spline method of Usmani [18]. For the choice of a ¼ 1=6, b ¼ 1=3 the parametric cubic spline reduces to ordinary cubic spline and our method reduces to the well-known Bickley scheme [8]. In Section 2, we develop the numerical technique for solving (1.1). Section 3 is devoted to the error analysis. In the last section, some numerical evidence is included to show the practical applicability and superiority of our method.
2. Parametric cubic spline method We develop a smooth approximate solution of (1.1) using parametric cubic spline function. For this purpose, we discretize the interval ½a; b using equally spaced knots xi ¼ a þ ih, i ¼ 0; 1; 2; . . . ; N , x0 ¼ a, xN ¼ b and h ¼ ðb aÞ=N ;
ð2:1Þ
where N is a positive integer. A function Sðx; sÞ of class C 2 ½a; b which interpolates yðxÞ at the mesh points xi , depends on a parameter s, reduces to cubic spline in ½a; b as s ! 0 is termed as parametric cubic spline function. The spline functions Sðx; sÞ ¼ SðxÞ satisfying in ½xi ; xiþ1 , the differential equation
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
177
S 00 ðxÞ þ sSðxÞ ¼ ðS 00 ðxi Þ þ sSðxi ÞÞðxiþ1 xÞ=h þ ðS 00 ðxiþ1 Þ þ sSðxiþ1 ÞÞðx xi Þ=h;
ð2:2Þ
where Sðxi Þ ¼ yi and s > 0 is termed as cubic spline in compression. Solving (2.2) and determining the constants of integration from the interpolatory conditions Sðxi Þ ¼ yi and Sðxiþ1 Þ ¼ yiþ1 we get, after writing x ¼ hs1=2 h2 xðx xi Þ xðxiþ1 xÞ S 00 ðxiþ1 Þ sin SðxÞ ¼ 2 þ S 00 ðxi Þ sin h h x sin x 2 2 h ðx xi Þ 00 x þ 2 S ðxiþ1 Þ þ 2 Sðxiþ1 Þ h x h 2 ðxiþ1 xÞ 00 x þ S ðxi Þ þ 2 Sðxi Þ : ð2:3Þ h h Differentiating (2.3) and using continuity conditions which lead to the tridiagonal system h2 ðaMi1 þ 2bMi þ aMiþ1 Þ ¼ yiþ1 2yi þ yi1 ; 2
ð2:4Þ 2
00
where a ¼ ðx cosec x 1Þ=x , b ¼ ð1 x cot xÞ=x , Mi ¼ S ðxi Þ, i ¼ 1ð1Þ N 1. The condition of continuity given by (2.4) ensures the continuity of the first order derivatives of the spline Sðx; sÞ at interior nodes. We may write (1.1) in the form of Mj ¼ fj yj þ gj and substituting in Eq. (2.4), we get the following tridiagonal system which gives the approximations y1 ; y2 ; . . . ; yN 1 of the solution yðxÞ at x1 ; x2 ; . . . ; xN 1 . ð1 þ ah2 fi1 Þyi1 þ ð2 þ 2bh2 fi Þyi þ ð1 þ ah2 fiþ1 Þyiþ1 ¼ h2 ðagi1 þ 2bgi þ agiþ1 Þ
with y0 ¼ a0 ; yN ¼ a1 ; i ¼ 1ð1ÞN 1: ð2:5Þ
Remark 2.1. For a ¼ 1=6, b ¼ 1=3, the cubic spline in compression reduces to the ordinary cubic spline and the method (2.5) reduces to the well-known Bickley Scheme [8]. Remark 2.2. For a ¼ 1=12, b ¼ 5=12, the relation (2.5) is the same as the fourth order method of Usmani [18] based on quartic spline.
3. Convergence analysis In this section, we investigate the error analysis of the spline method discussed in Section 2. We can write our method in matrix–vector form as
178
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
AY ¼ R þ T ;
ð3:1Þ
AY ¼ R;
ð3:2Þ
AE ¼ T ;
ð3:3Þ
where E ¼Y Y
and
A ¼ A0 þ Q;
ð3:4Þ
Q ¼ h2 BF , F ¼ diagðfi Þ, i ¼ 1; 2; . . . ; N and A0 ¼ ðaij Þ is the tridiagonal matrix defined by 8 3; i ¼ j ¼ 1; N 1; > > < 2; i ¼ j ¼ 2; 3; . . . ; N 2; ð3:5Þ ai;j ¼ 1; ji jj ¼ 1; > > : 0; otherwise: The tridiagonal matrix B is given by
ð3:6Þ
For the vector R, we have 8 < ð1 ah2 f0 Þa0 h2 ag0 ; ri ¼ h2 ðagi1 þ 2bgi þ agiþ1 Þ; : ð1 ah2 fN Þa1 h2 agN ;
i ¼ 1; 2 6 i 6 N 2; i ¼ N 1:
ð3:7Þ
Note that the ith equation of the linear system is ðah2 fi1 1Þyi1 þ ð2 þ 2bh2 fi Þyi þ ðah2 fiþ1 1Þyiþ1 ¼ h2 ðagi1 þ 2bgi þ agiþ1 Þ þ ti ; where ti , i ¼ 1; 2; . . . ; N 1 are the local truncation errors given by h4 ð4Þ h6 ð6Þ ui þ ð1 30aÞ u : 12 360 i
ð3:8Þ
xi1 < ni < xi ; i ¼ 1ð1ÞN 1
ð3:9Þ
ti ¼ f1 2ða þ bÞgh2 u00i þ ð1 12aÞ For a ¼ 1=12, b ¼ 5=12 ti ¼ ð1=240Þh6 y ð6Þ ðni Þ þ Oðh7 Þ; is the truncation error.
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
179
Thus ktk ¼
h4 M4 ; 84
ktk ¼
h6 M6 ; 240
M4 ¼ max jy ð4Þ ðxÞj
ð3:10Þ
x
and M6 ¼ max jy ð6Þ ðxÞj; x
ð3:11Þ
where k k represent the 1-norm in matrix–vector. The explicit inverse of A0 , namely A1 0 ¼ ðgij Þ, is defined by, see Usmani [17]. 8 2ðN jÞ þ 1 > > ; i 6 j; < ð2i 1Þ 4N ð3:12Þ gij ¼ > > : ð2j 1Þ 2ðN iÞ þ 1 ; i P j; 4N 1 which indicates that A1 0 > 0, that is, all the entries of A0 are positive. In fact we have
kA1 0 k6
N2 þ 1 : 8
ð3:13Þ
The equality in (3.13) holds only if N is odd. Using (2.1) the inequality (3.13) can be written in the form kA1 0 k6
ðb aÞ2 þ h2 : 8h2
ð3:14Þ
The following result gives the sufficient condition for which the system (3.2) has a unique solution. Theorem 3.1. The discrete boundary value problem (3.2) has a unique solution if Kjf ðxÞj < 1;
2
where K ¼ ½ðb aÞ þ h2 =8:
ð3:15Þ
Proof. The proof follows by observing that kA1 0 Qk < 1 and the following statement. If W is an N N matrix and kW k < 1, then ðI þ W Þ is non-singular, see [9]. We now derive a bound on kek. From Eq. (3.3) we have 1 1 1 E ¼ A1 T ¼ ðA0 þ QÞ T ¼ ðI þ A1 0 QÞ A0 T and it follows that kek 6
kA1 0 kktk 1 kA1 0 kkQk
provided that kA1 0 kkQk < 1.
ð3:16Þ
180
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
Now since kQk < h2 jf ðxÞj, then (3.10), (3.14)–(3.16) give kek 6
KM4 h2 ffi Oðh2 Þ 84½1 Kjf ðxÞj
ð3:17Þ
and (3.11), (3.14)–(3.16) give kek 6
KM6 h4 ffi Oðh4 Þ: 240½1 Kjf ðxÞj
ð3:18Þ
These inequalities show that (3.2) is a second order convergent method for arbitrary a; b such that a þ b ¼ 1=2 and fourth order convergent method for a ¼ 1=12, b ¼ 5=12.
4. Numerical results and discussion We consider two numerical examples illustrating the comparative performance of our present method with some earlier numerical methods for solving (1.1). All the computations were carried out using double precision arithmetic in order to keep the rounding errors negligible as compared to the discretization errors. Problem 1 2 1 y ; 2 x x yð2Þ ¼ 0; yð3Þ ¼ 0 y 00 ¼
with the exact solution yðxÞ ¼ ½5x2 þ 19x 36=x=38. Problem 2 y 00 ¼ 100y; yð0Þ ¼ yð1Þ ¼ 1 with the exact solution yðxÞ ¼ coshð10x 5Þ= cosh 5. Problems 1 and 2 were solved using the parametric cubic spline method described in Section 2 with a number of h values. For the sake of comparison with other methods, we tabulate the results by some well-known methods. The observed maximum errors in absolute values are given in Tables 1–4. From Table 1 we can notice that if the step size h is reduced by a factor of 1/2, then the maximum errors are approximately reduced by a factor of 1/16. Thus the numerical results confirm that our present method is fourth order convergent for a ¼ 1=12, b ¼ 5=12. For all other values of a and b, the method is second
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
181
Table 1 Observed maximum errors for a ¼ 1=12, b ¼ 5=12 (fourth order method) h
Problem 1
Problem 2
1/8 1/16 1/32
1.74 · 107 1.09 · 108 6.85 · 1010
1.74 · 103 1.12 · 104 7.29 · 106
Table 2 Observed maximum errors for Problem 1 Our method
h ¼ 1=4
h ¼ 1=8
h ¼ 1=16
a ¼ 1=18, b ¼ 4=9 a ¼ 1=10, b ¼ 2=5 a ¼ 1=14, b ¼ 3=7 a ¼ 1=16, b ¼ 7=16 Alsaid [3] Alsaid [5] Sakai and Usmani [15] Albasiny and Hoskins [2] Cent-diff
5.14 · 105 3.50 · 105 2.05 · 105 3.79 · 105 5.49 · 105 1.60 · 104 7.93 · 105 1.65 · 104 2.79 · 104
1.36 · 105 8.46 · 106 5.74 · 106 1.01 · 105 1.87 · 105 2.66 · 105 2.06 · 105 4.17 · 105 5.42 · 105
3.46 · 106 2.09 · 106 1.47 · 106 2.59 · 106 5.07 · 106 5.58 · 106 5.20 · 106 1.04 · 105 1.19 · 105
Table 3 Observed maximum errors for Problem 1 Our method a ¼ 1=10, b ¼ 2=5 a ¼ 1=14, b ¼ 3=7 a ¼ 1=16, b ¼ 7=16 Alsaid [3]
h ¼ 1=6
h ¼ 1=12 5
h ¼ 1=24
6
1.49 · 10 9.78 · 106 1.75 · 105 3.08 · 105
9.29 · 107 6.60 · 107 1.15 · 106 2.28 · 106
3.74 · 10 2.61 · 106 4.60 · 106 8.84 · 106
Table 4 Observed maximum errors for Problem 2 Our method a ¼ 1=10, b ¼ 2=5 a ¼ 1=14, b ¼ 3=7 a ¼ 1=18, b ¼ 4=9 Alsaid [3] Sakai and Usmani [15] Khalifa and Eilbeck [12]
h ¼ 1=16 3
1.28 · 10 7.22 · 104 1.83 · 103 2.27 · 103 3.06 · 103 –
h ¼ 1=32 4
3.07 · 10 2.06 · 104 4.91 · 104 6.84 · 104 7.58 · 104 –
h ¼ 1=20 4
8.17 · 10 5.00 · 104 1.22 · 103 1.57 · 103 – 1.8 · 103
h ¼ 1=40 1.95 · 104 1.34 · 104 3.16 · 104 4.53 · 104 – 4.7 · 104
order convergent as shown in Tables 2–4. It may be noted that from Tables 2–4 that our present method produces numerical results that are better than those obtained by the other methods. From (3.8) and (3.9), we notice that the parametric cubic spline method has smaller local truncation error than the
182
A. Khan / Appl. Math. Comput. 154 (2004) 175–182
cubic and quadratic spline method of [2,3,12,13]. This is in agreement with our numerical results given in Tables 2–4. From the above discussions and experiments we can conclude that second- and fourth order parametric cubic spline method produce better results than most of the other methods.
References [1] J.H. Ahlberg, E.N. Nilson, J.L. Walsh, The Theory of Splines and their Applications, Academic Press, New York, 1967. [2] E.L. Albasiny, W.D. Hoskins, Cubic spline solution of two point boundary value problem, Comput. J. 12 (1969) 151–153. [3] E.A. Al-Said, Cubic spline method for solving two point boundary value problems, Korean J. Comput. Appl. Math. 5 (1998) 759–770. [4] E.A. Al-Said, Spline solution for system of second order boundary-value problems, Int. J. Comput. Math. 62 (1996) 143–154. [5] E.A. Al-Said, Quadratic spline solution of two point boundary value problems, J. Nat. Geom. 12 (1997) 125–134. [6] E.A. Al-Said, M.A. Noor, Al-Shejari, Numerical solution for system of second order boundary value problems, Korean J. Comput. Appl. Math. 5 (1998). [7] U.M. Ascher, R.M. Metheij, R.D. Russell, Numerical solution of boundary value problems for ordinary differential equations, SIAM, Philadelphia, PA, 1995. [8] W.G. Bickley, Piecewise cubic interpolation and two-point boundary-value problems, Comput. J. 11 (1968) 206–208. [9] C. Froberg, Numerical Mathematics, Theory of Computer Applications, Benjamin/Cummings, Reading, MA, 1985. [10] D.J. Fyfe, The use of cubic splines in the solution of two-point boundary-value problems, Comput. J. 12 (1969) 188–192. [11] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley, New York, 1961. [12] A.K. Khalifa, J.C. Eilbeck, Collocation with quadratic and cubic splines, IMA J. Numer. Anal. 2 (1982) 111–121. [13] G. Mullenheim, Solving two-point boundary-value problems with spline function, IMA J. Numer. Anal. 11 (1982) 503–518. [14] M.A. Noor, A.K. Khalifa, Cubic splines collocation methods for unilateral problems, Int. J. Eng. Sci. 25 (1987) 1527–1530. [15] M. Sakai, R.A. Usmani, Quadratic spline and two point boundary value problems, Publ. RIMS, Kyoto Univ. 19 (1983) 7–13. [16] M. Sakai, R.A. Usmani, On consistency relations for cubic spline on spline and their asymptotic error estimates, J. Approx. Theory 55 (1985) 195–200. [17] R.A. Usmani, Bounds for the solution of a second order differential equation with mixed boundary conditions, J. Eng. Math. 9 (1975) 159–164. [18] R.A. Usmani, M. Sakai, A connection between quartic spline and Numerov solution of a boundary value problems, Int. J. Comput. Math. 26 (1989) 263–273.