One parameter quadratic C1-spline collocation method for solving first order ordinary initial value problems

Applied Mathematical Modelling 23 (1999) 153±160

One parameter quadratic C 1-spline collocation method for solving ®rst order ordinary initial value problems S. Sallam *, M. Naim Anwar Kuwait University, Department of Mathematics and Computer Science, P.O. Box 5969, Safat 13060, Kuwait Received 4 December 1997; accepted 17 July 1998

Abstract The convergence and stability analysis of a ``variable'' quadratic C 1 -spline collocation method for solving the initial value problem y 0 …x† ˆ f …x; y†;

y…0† ˆ y0 ;

x 2 ‰0; bŠ

will be considered. Letting the interior (non-nodal) collocation point xk‡b ˆ xk ‡ bh be dependent on some parameter b 2 …0; 1Š, it will be shown that the proposed method is strongly unstable if b < 12 and it turns out that the method is a continuous extension of the well-known mid-point and trapezoidal methods, if b ˆ 12 and b ˆ 1, respectively. Moreover, a wider region of absolute stability is achieved if b ! 1ÿ . Error bounds in the uniform norm for ks…i† ÿ y …i† k; i ˆ 0; 1 if y 2 C 3 ‰0; bŠ; together with illustrative examples will also be presented. Ó 1999 Elsevier Science Inc. All rights reserved. AMS classi®cation: 65 L 05; 65 D 05 Keywords: Initial value problems; Quadratic spline; Collocation methods; Absolute stability; A-stability; Sti€-equations

1. Introduction In this paper we present an analysis of a particular one parameter collocation method, based on quadratic C 1 -splines, to ®nd approximate solution to the initial value problem: y 0 …x† ˆ f …x; y†;

y…0† ˆ y0 ; x 2 ‰0; bŠ:

…1†

Letting the interior (non-nodal) collocation point xkÿ1‡b ˆ xkÿ1 ‡ bh, be dependent on some parameter b 2 …0; 1Š, we well establish the so-called b-method where an interesting relation is revealed between the b-method and the well-known explicit mid-point rule, if b ˆ 12, which has no interval of absolute stability and is of order two, hence (see [1, 2]) the method may be regarded as a continuous extension of the mid-point rule for solving Eq. (1) in the sense that it reproduces the values generated by the mid-point rule at the mesh points. Moreover, the method provides a continuous approximation to the ®rst derivative as well. Similar results are obtained if b ˆ 1, and *

Corresponding author. Tel.: +965 481 1188/5604; fax: +965 481 7201; e-mail: [email protected]

0307-904X/99/$ ± see front matter Ó 1999 Elsevier Science Inc. All rights reserved. PII: S 0 3 0 7 - 9 0 4 X ( 9 8 ) 1 0 0 7 3 - 2

154

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

the b-method is a continuous extension of the trapezoidal rule (which is A-stable) and is of order two. Global methods for solving ®rst order ordinary initial value problems have been investigated in 1967 by Loscalzo (see [4]). Since then several authors have presented similar procedures based on the use of splines (see, for example [2,5] and the references cited therein). Recently, a class of generalized trapezoidal formulas for solving Eq. (1) has been introduced by Chawla [1]. Another linear one-step method which can be of value in the context of sti€ness is the so-called Theta method (see Lambert [3]). For a given positive integer N the interval [0, b] is partitioned into N equal subintervals Ik ˆ ‰xkÿ1 ; xk Š; k ˆ 1…1†N , with stepsize h ˆ b=N . Let p2 denote the collection of all real polynomials of degree at most 2, and 1 SN;2 ˆ fs…x†: s 2 C 1 ‰0; bŠ; s 2 p2 ; for x 2 Ik ; k ˆ 1…1†N g: …`†

If g: ‰0; bŠ ! R is a real-valued function and b 2 R; …0 < b 6 1† is a real parameter, then gk‡b stands for g…`† …xk ‡ bh† …k ˆ 0…1†N ÿ 1†; ` P 0. Now, given the real numbers s0k‡b …k ˆ 0…1†N ÿ 1†, …1† s0 and s00 , then the unique quadratic spline s 2 SN;2 de®ned in ‰xk ; xk‡1 Š will be s…x† ˆ sk ‡ hs0k A…t† ‡ hs0k‡b B…t†;

…2†

where A…x† ˆ x ÿ

x2 ; 2b

B…x† ˆ

x2 ; 2b …1†

and x ˆ xk ‡ th; t 2 ‰0; 1Š, with a similar expression for s…x† in ‰xkÿ1 ; xk Š. Since s 2 SN ;2 , then we have   1 h sk ˆ skÿ1 ‡ h 1 ÿ s0kÿ1 ‡ s0kÿ1‡b ; 2b 2b   …3† 1 1 0 0 0 ; k ˆ 1…1†N ; sk ˆ 1 ÿ s ‡ s b kÿ1 b kÿ1‡b and hence s is uniquely determined in [0, b]. Based on the quadratic spline scheme discussed above, we are now able to develop a global collocation method to ®nd an approximate solution s…x† to the exact solution y…x† of the initial value problem (1). Let f 2 …C 2 ‰0; bŠ  R) and satisfy the Lipschitz condition jf …x; y† ÿ f …x; y†j 6 Ljy ÿ yj; in ‰0; bŠ  R. The approximate solution s…x† will be constructed as follows: On ‰xkÿ1 ; xk Š; k ˆ 1…1†N skÿ1‡b ˆ skÿ1 ‡ 12hbs0kÿ1 ‡ 12hbs0kÿ1‡b ;

…4a†

where s0kÿ1‡b ˆ f …xkÿ1‡b ; skÿ1‡b †;    1 1 0 sk ˆ skÿ1 ‡ h 1 ÿ s ‡ f …xkÿ1‡b ; skÿ1‡b † ; 2b kÿ1 2b   1 0 1 0 sk ˆ 1 ÿ skÿ1 ‡ f …xkÿ1‡b ; skÿ1‡b †; k ˆ 1…1†N ; b b

…4b† …4c†

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

155

and skÿ1‡b , for k P 1, using the initial conditions s0 ˆ y0 ; s00 ˆ f …0; y0 †, can be found via Eq. (4a) by using a quasi-Newton method or the contraction mapping principle. It can be easily veri®ed, using the contraction mapping principle, that the quadratic C 1 -spline approximation to Eq. (1), is uniquely determined if Lhb=2 < 1. Remark 1. If b ˆ 1, formula (4b) reduces to the well-known classical arithmetic mean trapezoidal formula and hence the one-parameter b-method described by Eqs. (4a)±(4c) may be regarded as a generalization to the classical trapezaidal rule. Now, by virtue of Eqs. (4a)±(4c) we shall prove that the method is im-practical if b < 12 according to the following. Lemma 1. The quadratic C 1 -spline method de®ned by Eqs. (4a)±(4c) is strongly unstable, and hence divergent, if b < 12. Proof. Denote by sk the vector sk ˆ …sk ; hs0k †T . Then the application of the method to the test equation y 0 ˆ 0, yields   1 sk ˆ skÿ1 ‡ 1 ÿ hs0kÿ1 ; 2b   1 0 0 sk ˆ 1 ÿ …5† s ; b kÿ1 or in matrix notation (5) will be sk ˆ Askÿ1 and the coecient matrix A has spectral raduis …1 ÿ b†=b > 1, if b < 12. 2. Convengence and rate of convergence In this section, we proceed to obtain a priori error estimates, in the uniform norm, for the proposed b-method. It will be shown that the method is of O…h2 † if f 2 …C 2 ‰0; bŠ  R†. An estimation for the local truncation error (LTE) of the b-method can easily shown to be (by direct application of Taylor's expansion about xkÿ1‡b )     1 b3 1 b2 3 1 h ‡ O…h4 † ˆ …2 ÿ 3b†h3 ‡ O…h4 †: Tkÿ1‡b ˆ …1 ÿ b†3 ‡ ÿ 1 ÿ 6 6 2b 2 12 …r†

…r†

…r†

…r†

Setting ek ˆ sk ÿ yk , where yk ˆ y …r† …xk †, we have Lemma 2. je0kÿ1‡b j 6 Ljekÿ1‡b j;

k ˆ 1…1†N :

…6†

Proof. Immediate consequence of Lipschitz continuity of f . Remark 2. If b ˆ 12, the method which corresponds to the embedded mid-point rule, has LTE h3 …3† y …ak †, whereas if b ˆ 1, the method which correspond to the embedded trapezoidal rule, has 24 LTE

156

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

ÿ

h3 …3† y …bk †; ak ; bk 2 Ik : 12

Now, if y 2 C 3 ‰0; bŠ, then a bound for the global error ekÿ1‡b of the b-method can be given by   2 9 2 h jekÿ1‡b j 6 1 ÿ b ‡ 6b …7† ky …3† k1 ; 6 2 and consequently   2 9 h je0kÿ1‡b j 6 L 1 ÿ b ‡ 6b2 ky …3† k1 : 6 2

…8†

Since (using Eq. (4c))   1 0 1 bh2 …3† h2 y …ak † ÿ y …3† …bk †; e0k ˆ 1 ÿ ekÿ1 ‡ ekÿ1‡b ‡ 2 2 b b which is a ®rst order di€erence equation, with e00 ˆ 0, hence kÿj k  X 1 0 ek ˆ 1ÿ cj b jˆ1 and  kÿj k  X 1 1ÿ ; b

je0k j 6 c where

jˆ1

c ˆ max jcj j; j

  2 L 9 h2 2 h c6 ky …3† k1 ‡ …1 ‡ b† ky …3† k1 : 1 ÿ b ‡ 6b 6 2 b 2

< 1†, we have For b > 12 i.e., …1ÿb b ( ) L…1 ÿ 92 b ‡ 6b2 † ‡ 3b…1 ‡ b† 2 …3† 0 h ky k1 ; jek j 6 6…2b ÿ 1†

b 2 …12; 1†:

In particular if b ˆ 12 (the mid-point rule), we have je0kÿ1 j 6 L 2

h2 …3† ky k1 24

and je0k j 6 …L ‡ 9†

h2 …3† ky k1 12

and   1 h h3 h3 he0kÿ1 ‡ e0kÿ1‡b ‡ by …3† …ak † ÿ y …3† …bk † ek ˆ ekÿ1 ‡ 1 ÿ 4 6 2b 2b or ek ˆ ekÿ1

  1 h h3 h3 ˆ 1ÿ he0kÿ1 ‡ e0kÿ1‡b ‡ by …3† ÿ y …3† …bk †; 4 6 2b 2b

…9†

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

i.e.,

k X ek ˆ dj

)

157

jek j 6 kd;

jˆ1

thus

    2 1 9 h 2 ky …3† k1 jek j 6 xk L 1 ÿ b ‡ 6b ‡ 3b…1 ‡ b† 6 eb 2    2 L 1 h 3 3 2 ky …3† k1 ‡ …1 ÿ b† ‡ b ‡ 3b 1 ÿ 2 2b 2b 3 h ‡ …2 ‡ 3b† ky …3† k1 : 12

…10†

Theorem 1. Let y 2 C 3 ‰0; bŠ, then for all x 2 ‰a; bŠ, we have js…j† …x† ÿ y …j† …x†j < kj …b†h2 ;

j ˆ 0; 1;

where kj …b† denote generic constants independent of h, but dependent on b and upon the order of the derivative. Proof. On ‰xkÿ1 ; xk Š, we have e0 …x† ˆ s0 …x† ÿ y 0 …x† ˆ s0 …x† ÿ u0 …x† ‡ u0 …x† ÿ y 0 …x†, where u0 …x† is the linear interpolant of y 0 …x† at xkÿ1 and xkÿ1‡b ; i.e.,   t t 0 y 0 …xkÿ1 † ‡ y 0 …xkÿ1‡b †; u …x† ˆ 1 ÿ b b hence   t 0 t 0 0 ekÿ1 ‡ e0kÿ1‡b ; s …x† ÿ u …x† ˆ 1 ÿ b b and js0 …x† ÿ u0 …x†j 6 je0kÿ1 j ‡ je0kÿ1‡b j; or

( 0

0

js …x† ÿ u …x†j 6

) 2bL…1 ÿ 92 b ‡ 6b2 † ‡ 3b…1 ‡ b† 2 …3† h ky k1 : 6…2b ÿ 1†

Using the fact that u0 …x† is the linear interpolant of y 0 …x† at xkÿ1 ; xkÿ1‡b , then b2 h2 …3† ju …x† ÿ y …x†j 6 ky k1 : 8 Hence ( ) 2 2 9 16bL…1 ÿ b ‡ 6b † ‡ 24b…1 ‡ b† ‡ 6b …2b ÿ 1† 2 h2 ky …3† k1 ; je0 …x†j 6 48…2b ÿ 1†     9 b2 2 2 6 8L 1 ÿ b ‡ 6b ‡ 3…2b ‡ 3b ‡ 4† h2 ky …3† k1 : 24…2b ÿ 1† 2 On ‰xkÿ1 ; xk Š; Zx e…x† ˆ e0 …x† dx ‡ ekÿ1 ; 0

0

xkÿ1

can be determined using Eqs. (10) and (11). This completes the proof of Theorem 1.

…11†

158

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

3. Stability analysis Applying the method to the test equation y 0 ˆ ky;

k 2 C;

Rek < 0…kh ˆ z†

we get b zb skÿ1‡b ˆ skÿ1 ‡ h s0kÿ1 ‡ skÿ1‡b ;  2  2 1 z s0kÿ1 ‡ skÿ1‡b ; sk ˆ skÿ1 ‡ h 1 ÿ 2b 2b   1 z hs0k ˆ 1 ÿ hs0kÿ1 ‡ skÿ1‡b ; b b consequentely

   2b ÿ b2 z ‡ z 1 z skÿ1 ‡ 1 ÿ ‡ hs0kÿ1 ; sk ˆ b…2 ÿ bz† 2b 2…2 ÿ bz†    2z 1 z hs0k ˆ skÿ1 ‡ 1 ÿ ‡ hs0kÿ1 : b…2 ÿ bh† b 2 ÿ bz

…12†

Or in matrix notation, Eq. (12) will be sk ˆ Mskÿ1 ; where

0 M ˆ@

2bÿb2 z‡z b…2ÿbz†

2bÿb2 z‡bzÿ1 b…2ÿbz†

2h b…2ÿbz†

2bÿb2 z‡2bzÿ2 b…2ÿbz†

1 A:

If b ˆ 1 (which leads to the classical trapezoidal rule), then M reduces to ! 2 2ÿz 2z 2ÿz

1 2ÿz z 2ÿz

with eigenvalues l1 ˆ 0 and l2 ˆ …1 ‡ 2z †=…1 ÿ 2z † and hence the method is A-stable. De®nition 1. The b-method de®ned by Eqs. (4a)±(4c) is said to be absolutely stable in a region R of the complex plane if, for all z 2 R, the eigenvalues l1;2 of M, satisfy jl1;2 j < 1: The intervals of absolute stability were deduced numerically for di€erent values of b 2 ‰12 ; 1Š, as shown in Table 1. The region of absolute stability enlarges to occupy the left-half of the complex plane, as b ! 1ÿ . 4. Numerical examples The purpose of this section is to present two examples to illustrate computationally the results established in Theorem 1. In these examples the following algorithm is used: · Set s0 ˆ y0 . · Use Eq. (4a) to compute skÿ1‡b ; k ˆ 1…1†N.

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

159

Table 1 b

Interval of absolute stability

0.50 0.60 0.75 0.90 1.00

No interval (vanishing interval) ()0.833, 0) ()2.666, 0) ()8.88, 0) ()1, 0)

· · · See

Use Eqs. (4b) and (4c) to compute sk and s0k , respectively. Use Eq. (2) to compute s…x† in each subinterval ‰xkÿ1 ; xk Š; k ˆ 1…1†N : s0 …x† is computed from the derivative of Eq. (2). for example Table 2 and Table 3.

Example 1. First, we consider the sti€ equation y 0 ˆ ÿ10…10 ÿ x†y;

y…0† ˆ 1; x 2 ‰0; 1Š

whose exact solution is   x  y…x† ˆ exp ÿ 10x 10 ÿ : 2 Example 2. In this example we consider the problem y 0 ˆ ÿ50…y ÿ x3 † ‡ 3x2 ;

y…0† ˆ 0; x 2 ‰0; 1Š;

3

with exact solution y…x† ˆ x :

Table 2

Absolute maximum error for s(x) and its continuous derivative for Example 1 Error bound For ks ÿ f k1 ks0 ÿ f 0 k1

b ˆ 0:5 3.6 8.3 1.2 4.0 1.7 5.2

´ ´ ´ ´ ´ ´

b ˆ 0:75 17

10 1016 1016 1020 1020 1019

1.7 2.4 4.1 1.8 4.7 1.2

´ ´ ´ ´ ´ ´

ÿ3

10 10ÿ4 10ÿ5 100 10ÿ1 10ÿ1

b ˆ 0:8 2.0 2.9 3.9 1.5 3.9 9.9

´ ´ ´ ´ ´ ´

bˆ1 ÿ3

10 10ÿ4 10ÿ5 100 10ÿ1 10ÿ2

5.1 1.2 3.1 1.4 4.2 1.1

´ ´ ´ ´ ´ ´

h ÿ3

10 10ÿ3 10ÿ4 100 10ÿ1 10ÿ1

0.004 0.002 0.001 0.004 0.002 0.001

Table 3

Absolute maximum error for s(x) and its continuous derivative for Example 2 Error bound For

b ˆ 0:5

ks ÿ f k1

7.2 1.0 1.2 7.6 2.2 4.9

ks0 ÿ f 0 k1

´ ´ ´ ´ ´ ´

102 102 101 105 105 104

b ˆ 0:75 3.3 8.4 3.0 1.4 3.1 7.6

´ ´ ´ ´ ´ ´

10ÿ8 10ÿ9 10ÿ9 10ÿ5 10ÿ6 10ÿ7

b ˆ 0:8 1.5 5.3 1.4 1.0 2.5 6.1

´ ´ ´ ´ ´ ´

10ÿ8 10ÿ9 10ÿ9 10ÿ5 10ÿ6 10ÿ7

bˆ1 1.7 4.4 1.1 1.1 2.9 7.6

´ ´ ´ ´ ´ ´

h 10ÿ7 10ÿ8 10ÿ8 10ÿ5 10ÿ6 10ÿ7

0.004 0.002 0.001 0.004 0.002 0.001

160

S. Sallam, M. Naim Anwar / Appl. Math. Modelling 23 (1999) 153±160

5. Concluding remarks In summary, we established a ``variable'' C 1 -quadratic spline procedure, that depends on a parameter b 2 …0; 1Š; for solving ®rst order ordinary initial value problems. The b-method is a continuous extension of the explicit midpoint and the trapezoidal rules, if b ˆ 12 and 1, respectively. Moreover, the method is shown to be of order two and numerical results indicate that b ˆ 0.8 may be the optimal value. Also a wider region of absolute stability is achieved, as b ! 1ÿ : References [1] M.M. Chawla, M.A. Al Zanaidi, D.J. Evans, A class of generalized trapezaidal formulas for the numerical integration of y 0 ˆ f …x; y† , Intrn. J. Computer Math. 67 (1996) 131±142. [2] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Di€erential Equations I, Nonsti€ Problems, Springer, New York±Berlin±Hiedelberg, 1987. [3] J.D. Lambert, Numerical Methods for Ordinary Di€erential Systems, Wiley, Chichester, New York, 1991. [4] F.R. Loscalzo, An introduction to the application of spline functions to initial value problems, in: T.N.E. Greville (Ed.), Theory and Applications of Spline Functions, Academic Press, New York, 1969, pp. 37±64. [5] S.P. Norsett, Splines and collocation for ordinary initial value problems, in: S.P. Singh et al. (Eds.), Approximation Theory and Spline Functions, Reidel, Dordrecht, 1984, pp. 397±4178.