On the convergence of finite difference methods for a class of two-point boundary value problems with periodic boundary conditions

li N O g ~ - ItOUAND

On the C o n v e r g e n c e of Finite Difference M e t h o d s for a Class of T w o - P o i n t B o u n d a r y Value P r o b l e m s w i t h Periodic B o u n d a r y Conditions C. P. K a t t i a n d S. N. B a b o o

School of Computer and Systems Sciences Jawaharlal Nehru University New Delhi, 110067 India and S. S i v a l o g a n a t h a n

Department of Applied Mathematics University of Waterloo Waterloo, Ontario, N2L 3G1 Canada

Transmitted by Melvin Scott

ABSTRACT In this paper we present a detailed convergence analysis to establish second-order convergence for a three-point two parameter family of finite difference methods (based on one evaluation of f ) for the class of two-point boundary value problems (TPBVPs)

y" + f ( x , y ) = O , y(0) = y(1),

0~
y'(O) = y ' ( 1 )

under suitable conditions on 3 f / a y and for y ~ C~[O, 1]. Some numerical examples axe provided to illustrate the results.

APPLIED MATHEMATICS AND COMPUTATION 75:287-302 (1996) © Elsevier Science Inc., 1996 655 Avenue of the Americas, New York, NY 1 0 0 1 0

0096-3003/96/$15.00 SSD10096-3003(95)00144-7

288 1.

C.P. KATTI ET AL. INTRODUCTION

We consider the class of two-point boundary value problems with periodic boundary conditions: y" + f ( x , y) = O,

0 ~< x ~< 1,

y(O) = y(1),

y'(O) = y'(1)

(1)

under the assumption that f has continuous partial derivatives with respect to y and that either

u,

af

(2)

or

0
~f

(3)

holds, where

sup

~f

0 < z < l o~Y

and

u, =

inf

~f

-0~
-oo < y < o o .

(4)

It is known from [11] that [0, 1] is then an oscillation free interval and (1) has a unique solution y(x) ~ 610, 1]. Accurate and fast numerical solutions of two-point boundary value problems with periodic boundary conditions are necessary in many important application areas in physics and engineering. High-order finite difference schemes and their convergence analysis for two-point boundary value problems with natural and mixed boundary conditions have been given in [2-7, 12, 13, 15, 16]. Although the finite difference scheme (based on one evaluation of f ) for the solution of (1) is well known, its convergence analysis does not seem to have been considered in the literature so far. In the following, we present a convergence analysis to establish O(h 2) convergence of a three-point two parameter family of finite difference methods for the solution of (1). For particular values of the parameters, the scheme reduces to the typical finite difference scheme for the solution of (1).

289

Convergence Methods for Periodic Problems

2.

F I N I T E D I F F E R E N C E M E T H O D A N D ITS C O N V E R G E N C E

For a p o s i t i v e integer N > t 4 , let h = 1 / N , x k = kh, k = 0 , 1 , 2 , . . . , N and set Yk = Y(Xk), fk = f(xk, Yk). Discretization of (1) leads to the standard three-point scheme: -Yk-1

k = 1,2,...,

+ 2Yk - Yk+l - h2fk = t k ( h )

N.

(5)

Where Y0 = YN and tk(h) = O(h4). It is assumed t h a t y ~ C4[0, 1]. For the discretization of the b o u n d a r y condition Y'0 = Y'N, we obtain the following equation --Yl -- YN-1 -[- 2 y N -- h 2 [ ( 2

--

~ ) f N "}- ~ f N - 1 "{- ~ f l + ( Y2 - IJ.) f ( Xo, YN)] = tN

(6)

where

tN(h) = h ~ [ ( ' - - ~ ) Y ( ~ ) ( X o )

+(~--~)Y(~)(~N)]

+ 0 ( t ' ~)

(7)

and t~, tt are arbitrary parameters. For p~ = tt -- 0, the equations (5)-(6) reduce to the simply, typical finite difference scheme for the solution of (1). W e m a y note t h a t for t~ = 8 = ~ the local truncation error is uniformly O(h4), whereas for /~ = 8 = 0, /~ = 8 = ~1 the error in tN(h) is O(h 3) if y(3)(0) ¢ y(3)(1). It m a y be noted that, computationally finite difference scheme (5)-(6) is 1 slightly superior when /~ = 8 = ~ as seen in Section 3. T o discuss our method we first write the system (5) and (6) in matrix form. Let D = (d~j)i, j = 1, N denote the quasitridiagonal m a t r i x with i = 1,2,...,

dii = 1,

i = 1,2,...,

N-

N,

1,

d~,~± 1 = - 1 ,

dl, N -~ dN, 1 =

--1

and let

Y = ( y,, ~ . . . . . yN)~', T = (tl(h),

a(Y) ...,

= ( g~, g~ .... , g~)~,

t~(~))

~

290

C. P. KATTI ET AL.

where g*= fi,

i= l,2,...,

N-1,

g N = (½ -- 8 ) f N + SfN_l + # f l + (½ -- p t ) f ( Xo, ym).

(8)

Thus our discretizations (5), and (6) can then be written in matrix form as DY-

h 2 a ( Y ) = T.

(9)

Our method providing approximations l~ = ( Yl, Y2,-.., ~tN)T for Y for (1) consists of solving the N × N system DY-

h2a( Y ) = O.

(10)

The nonlinear system (10) can be solved by the use of Newton's method and an adaptation of Gauss elimination for linear systems. In the following we establish O(h2)-convergence of our method for all f such that the inequality (2) or (3) is satisfied. Let E - ( % e 2 , . . . , eN) T = Y - ]?. We may write fk - ]k = f ( xk, Yk) - f ( xk, Yk) = ukek,

k= 1,..., N

(11)

and f ( x0, YN) -- f ( Xo, ~IN) = ~tgeg

(12)

where

~f

~f

Yk < 7/~ < Yk,

YN < TIN < ~]N"

Using (8), (11), and (12) from (9) and (10) we obtain the error equation: ( D-

h 2 M ) E = T.

(13)

291

Convergence Methods for Periodic Problems

Where M is the N X N quasi-diagonal matrix given by M = (

mij ) i = j,

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

Izua,

i = N,

j=l,

6UN_I,

i = N,

j=N-1,

(½--6)UN+(½--tZ)UN,

i=N,

j=-N,

0

otherwise.

-ui,

~-

1,

In the following we let

[o,½l.

(14)

Case(i). u , ~< 6 f / 6 y < u* < 0 . Let M* be the N × N quasi-diagonal matrix given by

M* = ( m ' j ) =

u*, /~u*, 6u*, [~1 - ( 6 + ~ ) ) u *

i = j, i = N, i = N,

j=l,

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

i=N,

j=

1,

j=N-1, N,

otherwise. Since u k ~< u* and /*, 6 satisfy condition (14), we have D - h2 M >I D - h2 M *.

Let

ho = m i n

[1 {11}] ~

i / ~ - ' !/~

"

15,

Then for h < h 0 both D - h 2M and D - h 2M* are irreducibly diagonally dominant with off-diagonal entries negative. It then follows (see e.g. Varga [17]) that D - h=M and D - h2M * are both monotone and (O-

h2M) -1 ,< ( D -

h2V*) -1.

(16)

292

C. P. KATTI ET AL.

We next establish the following result:

LEMMA 2.1. F o r 0 <. x <<. 1, - ~ < y < ~ let u , <<. 6 f / 6 y ~< u* < 0, /z, 6 ~ [0, ±]2 and h 1 = min(h0, ¼) where h o is given by (15). L e t B = ( D h2M*) -1 = (b~j)i, j = 1, N. T h e n f o r all h < h 1 C m a x bi N < - i<~i<~N h

(17)

where c is a positive c o n s t a n t which depends only on u*, tz and 6.

PROOF.

Let the matrix

A = D - h2 M * = ( a~j)i, j = 1, N be written

as al N

A =

All

aN-1, N aN1 "'" aN, N - 1

aNN

All is symmetric, tridiagonal, irreducible, diagonally d o m i n a n t with nonpositive off-diagonal elements. It follows t h a t A n is monotone [17]. Following similar arguments for inverting a tridiagonal m a t r i x as given in Henrici [10] we obtain -1

N

=

s i n h ( i O ) , sinh( N - j ) 0 s i n h ( 0 ) , sinh(NO)

'

i < j

(18)

where

sinh(O/2) = (h -¢-:--U/2).

(19)

We next proceed to obtain biN, i = 1 . . . . , N. Following the a r g u m e n t s for inverting a m a t r i x as given in [18], we obtain bNN = 1 / F

Convergence Methods for Periodic Problems

293

where

r = a N N - (aN~, a N 2 , . . . , aN, N - , ) A ; , I ( a ~ N , . . . ,

ag-,,N) T

and

bzN= --(1/F) A711(alN,...,ag_,,g) T,

i---1,2,...,N-1.

(20)

Substituting for A~-I1 from (18) into (20) we obtain cosh( N - 2 i) 0/2 b,g =

J

,

i = 1,2,...,

N,

(21)

where J = [2 + h2u*(/z + 8) - h2u*]cosh(NO~2)

-[2

+ h2u*(/z + 8 ) ] c o s h ( N -

20/2).

(22)

Substituting for h2u * from (18) into (22), after some simplification, we obtain J = 4 sinh(0/2) [sinh(( N - 1)0/2) + 2 p s i n h 2 ( O / 2 ) s i n h ( ( N - 1)0/2) +cosh((g-

2) 0/2)sinh( 0/2)]

(23)

where p = 1 - ( ~ + 8). Since sinh and cosh are monotonically increasing functions on [0, oo] it follows from (23) that J >/4 sinh( 0/2)sinh(( N - 1) 0/2).

(24)

294

C.P. KATTI ET AL.

From (19), we now have, 1( -

=

=

-

sinh-

2

~o 1

-

1

~/1 + t 2

dt

h i %-- u*

>1 h2u *

11/2

+-7and hence

y/•--U* Co- ~* >/NO t>

1/2-

(25)

h~u*

1

4

Use of (25) in (24) leads to

J >~ 2h~ f-L-- u* sinh

4-:- u* (1 - hi)

(26)

for all h < h 1. Substituting for J from (26) into (21) and maximizing over i we obtain (17).

THEOREM 2.1. L e t f, t~, ~ satisfy the hypothesis o f L e m m a 2.1. Also, let A = D - h2M * = ( a O N y = l , B = A -1 = (b0Nj=I andho, h 1 be as defined in L e m m a 2.1. Then f o r all h < h 1 our method defined by (10) is second order convergent and ]l Ell~ <~ K l h 2

where K 1 is a positive constant which depends only on u*, tt and B.

(27)

Convergence Methods/or Periodic Problems

295

PROOF. Let st denote the sum of the elements in the jth row in the matrix D - h 2M*. Then N

N

N

E biysj= E b o E j=l

j=l

N

ajk= E

k=l

N

N

E bijajk= E 6i~= 1.

k=l j=l

(28 /

k=l

Since sj >/ - h2 u* for 1 ~< j ~ N from (30) we have N

N

- h 2 u * ~ b0 <~ ~ bijs j = 1 j=l

j=l

and hence N-1

E

j=l

N

1

b~j ~< E bij ~< _ h2 u-------7"

j=l

(29)

Writing the error equation (13) componentwise and using (16) we obtain N-1

leil <<. ~

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

bijltjl + biyltg[,

(30)

j=l

Thus, with the help of (5), (7), (17), and (29) from (30) we obtain (27). Case(ii). 0 < u , ~< 6 f / 6 y < u* < Ir 2. We establish the following result:

N LEMMA 2.2. Let P = D - h2 M = ( pij)i,j=l, 0 <~ ~, ~ <~ ~1 as before and H = [12/qr2(1 - u*/qr2)] 1/2. Then P is non-singular for h ~ H. Also if Q = p-1 = (qo)i, j= 1, then

max

I<~
I qiNI

]~2 <~ - -

h

(31)

and N- 1

max E I q J <

l~i<~N j = l

ka

(32)

296

C.P. KATTI ET AL.

for sufficiently small h where k 2 and k 3 are positive constants which depend only on u*, u , , /z and &

PROOF.

Let the matrix P be written as

P1N

P1 a p =

P N a,- N ...

,

where P l a --- D l a - h2Mll, On and M n being submatrices of order N - 1 consisting of first ( N - 1) rows and columns of D and M respectively. Let P u = Dn - h2u*Iand Pu = Dn - h2u. I. The fact that Pla, t511 and P u are invertible for all h ~< H follows from [9]. Denoting P i l l = ( Pij)i,j=a, " N-1 -, N-a P~a-l a _ - - ( Pij)i,j=a and p~a = ( P- , ij)i,N -j =1 a it is obvious that p~a and p~a are symmetric and

P , ij < [9,j < Pij,

i, j = 1 , 2 , . . . , N -

1.

(33)

As in [8] we have

15*j = sin( N - j) ~ba sin( i~b1) sin( N~bl)sin ~b1 '

i ~< j,

sin( ~ba/2 ) = ( h ~ i - / 2 ) (34)

and sin( N - j) ¢2 sin( i~b2) ~9, ij

sin( N~b2)sin 4)2

i<~j,

sin(~b2/2 ) = ( h u ~ . / 2 ) . (35)

297

Convergence Methods for Periodic Problems

Let Q = (q#)Nj= 1 be defined by U-1

( qij)i, j=l = Pu I + Pul[ P l N , ' ' ' , PN-1, N]T~'-'[ PN1,''', PN, N-1]Pu 1,

(36) [ qlN,''', qN-1,N]T= _p~l[ PlN,''', PN-1, N] T~-I,

(37)

[qN,,''-, qN, N-l] = --~-'[ PN1,''', PU, N-,]P5 ~,

(38)

and qNN =- 'r- 1,

.r= PUN- [ P r o , ' ' ' , PN, N - , ] P u ' [ P1N, "'', PN-1, N]"

(39) To show that Q exists, it is enough to show that ~"4= 0. Since P = D - h2M = (Pij) from (39) we have T~--- 2 -- ( P l l "~ Pl, N-1 "Jr"PN-I,1 -[- PN-1, N-1) -

h2[(1/2 - 6 ) u N + (1/2 - /z)5 N + ~t£Ul(p, 1 -I" Pl, N - , )

"JC'SUN-I(PN-I,1-]- BN-1, N-1)] ~< 2 -

- (~,1, + ~,1, u-1 + ~, u-,,~ + ~, u-l, u-,) h2u.[(1 - 6 - / z ) + / z ( 1 f . u + i5.1, N_,)

+8(~,u_1,, + ~,u-l,u-1)] (using (33))

(40)

--= - 2 h ui/-~-, [ s i n ( ( N = 1)~b2/2) ] cos( N¢~2/2 ) ]

- A2u, [(1 - 8 - ~) + ( ~ + 8) (using (35)).

cos( N - 2) ~b2/2

cos(NCJ2) (41)

298

C.P. KATTI ET AL.

Since, from [8] it follows that, ~ < N~b1 < 7r and ~ ~< N¢2 < 7r for all h ~< H, both the quantities in the parenthesis of (41) are bounded by positive constants. It then follows easily that ~- < 0 for h ~ H. It can then be seen that (see e.g. [12]) PQ

=

QP

= I.

Where I is the identity matrix of order N. Hence P is nonsingular having inverse Q. To prove the next part of the theorem, from (40), we note that

ITI ~ [(PSI1 "1-P$1, N-1 "l-PSN-I,1-1- PSN-I,N-1) - 2]

(42)

Using (34) and (35) in (42) and after some simplification we obtain, for sufficiently small h, ] r l > h u ~ , [sin( u ~ , / 2 ) ] .

(43)

Since, 1

1

I

IqiN[ = T(Ril + ~3iN) ~ ~T~(])~I "~ ~3i$,N-1)

<

cos(( N - 2 i)~b~/2)

]~'ICOS(N~bl//2)

< [r[cos(fll/2)

(31) follows by the use of (43).

'

fll =

sup (N4h). O<<.h<.H

(44)

299

Convergence Methods for Periodic Problems

To establish (32), we note that for 1 ~< i, j ~< N small h N- 1

N~:~1

~b1

l < imax ~ < N - 1 jE = l Lqijh 4 sin(N~b,)

sin(~bl)

cos(( N - 2 i) ¢1//2)

COS(N~l/2 ) X

N2

~

1 and sufficiently

1 + ITI •1

N(N-

1)

2 sin(N~b~)

[2+ u*(tz+6)} "16

Since

1 < sin(N~bl)

< 1+

sl,

1<

sin~b 1 < 1 + s 2

where ~1, ~2 are positive constants depending only on u*, (32) follows by the use of (43).

THEOREM 2.2. Let 0 < u . < u* < 1r 2 and ~, ~ ~ [0, 1/2]. As before let P = D - h2M = (Pij)i,j=l, Q -= (D -- h2M) -1 =-(qij). Then, for suffiN ciently small h our method defined by (10) is second-order convergent and IIEI1~ ~< k, h2

(45)

where k 4 is a positive constant which depends on u*, u , , I~ and ~ only.

PROOF. The result follows easily by the use of Lemma 2.2 and analogous arguments as used in proving Theorem 3.1. 3.

NUMERICAL EXAMPLES

We illustrate the O( h~)-convergence of the finite difference method with the following examples.

300

C.P. KATTI ET AL.

EXAMPLE 3.1. We consider the linear two-point b o u n d a r y value problem studied by Ciarlet et al. [9]:

(

y" - y - s i n ( 2 1 r x ) ( - 1 - 4 w

2) x 3 -

(

8

- 4Ir cos(27rx) 3 x 2 - - x + 3 y(0) = y ( 1 ) ,

4

-x 2 + 3

+6x-

-

:_)

= O,

(46)

(47)

y'(0) = y'(1)

with exact solution

y(x) =

( x 3 - -4x 2 + 3

sin(21rx).

(48)

For this example 8 f / f l y = - 1 and y(3)(0):/: y(3)(1). For N = 2 k, k = 1 2, 3 , . . . , 9 the corresponding values of ]l EI]~ are calculated for t~ = 8 = and tL = 8 = 0. The results are summarized in Table 1. 1 It m a y be noted t h a t the values of II EIJ~ for /z = 8 = ~ are slightly superior to those of t¢ = 8 = 0.

EXAMPLE 3.2.

We consider the solution for the equation

y,, _ y _ ~y3 = 0.8sin(0.27 × l l . 5 2 w x ) ,

w = 0.92845

(49)

with periodic b o u n d a r y conditions y(0) -- y ( 1 ) ,

(50)

y'(0) = y'(1).

TABLE 1 NORM OF THE ERROR FOR VARIOUS MESH REFINEMENTS AND PARAMETER VALUES

N iz=8=~ ]]EII:¢

8

16

32

64

128

256

4.7(- 3)

9.6(- 4)

2.0(- 4)

4.4(- 4)

1.0(- 5)

2.4(- 6)

6.3(-3)

2.2(-3)

6.1(-4)

1.6(-4)

4.0(-5)

1.0(-5)

1

1~ = 8 = 0

HElL

Convergence Methods for Periodic Problems

301

TABLE 2 NORM OF THE ERROR FOR VARIOUS MESH REFINEMENTS AND PARAMETER VALUES

N

8

16

32

64

128

256

/ x = t t --- ~1 IIEll~ ~=~=0

4.7( - 3)

1.3(- 3)

3.2(- 4)

8.0(- 5)

1.9(- 5)

3.9(- 6)

]lEll~

6.3(-3)

1.6(-3)

3.9(-4)

9.6(-5)

2.3(-5)

4.6(-6)

Here 3 f / 3 y = - 1 - ly2 ~< _ 1 and y3(0) ~: y(3)(1). The solution of the nonlinear problem (49) with associated periodic b o u n d a r y condition (56) was c o m p u t e d and the error norms (llEll~) for N = 2 k, k = 2,3 . . . . ,9 are presented in Table 2.

EXAMPLE 3.3. As a last example we consider the nonlinear two-point b o u n d a r y value problem studied by Ciarlet et al. [9]: y,, = y + y3 + esin(2~rx)[47r2 cos2(2,rrx) _ 4zr 2 sin(2"rr x)

--

e 2sin(2"n'x)

0 < x< 1 y(0) = y(1),

y'(0) = y'(1)

--

1],

(51) (52)

with exact solution y(X)

=

e sin(2~rx).

(53)

For this example 6 f / 3 y = - 1 - 3 y2 ~ _ 1 and y(3)(0) = y(3)(1). For N = 2 k k = 2, 3 , . . . , 9 the corresponding values of IIEll~ are shown in Table 3.

TABLE 3 NORM OF THE ERROR FOR VARIOUS MESH REFINEMENTS AND PARAMETER VALUES

N /z=8=~

8

16

32

64

128

256

2.4(-1)

4.1(-2)

7.6(-3)

1.6(-3)

3.5(-4)

8.7(-5)

1.0(- 1)

2.4(- 2)

5.8(- 3)

1.5(- 3)

3.6(- 4)

9.1(- 5)

l

IIEIl~ / z - (~ = 0

IIElla

302

C.P. KATTI ET AL.

REFERENCES 1 O.R. Asfar and A. M. Hussein. Numerical solution of linear two-point boundary value problems via the fundamental matrix method, Int. J. Numerical Methods in Engineering 28:1205-1216. 2 A. Boutayeb and E. H. Twizell. Finite difference methods for twelfth-order boundary value problems, J C A M 35:133-138 (1991). 3 J . R . Cash and A. Singhal. High-Order methods for the numerical solution of two-point boundary value problems B I T 22:184-199 (1982). 4 M.M. Chawla. A sixth-order tridiagonal finite difference method for nonlinear two-point boundary value problems, B I T 17:128-133 (1977). 5 M.M. Chawla. An eight-order tridiagonal finite difference method for nonlinear two-point boundary value problems, B I T 17:281-285 (1977). 6 M. M. Chawla. High accuracy tridiagonal finite difference approximation for nonlinear two-point boundary value problems, J. IMA 22:203-209 (1978). 7 M. M. Chawla. A sixth-order tridiagonal finite difference method for general two-point boundary value problems with nonlinear boundary conditions, J. IMA 24:35-42 (1979). 8 M.M. Chawla and C. P. Katti. A finite difference method for a class of singular two-point boundary value problems, IMAJNA 4:457-466 (1984). 9 P.G. Ciarlet, M. H. Schultz, and R. S. Varga. Numerical methods of high-order accuracy for nonlinear boundary value problems IV (Periodic b.c.'s), Numerishe Mathematik 12:266-279 (1968). 10 P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley & Sons, New York, 1962. 11 L . F . Shampine. Boundary value problems for ordinary differential equations, SIAM J. Num. Ana. 2:219-242 (1968). 12 R . P . Tewarson. A seventh-order numerical method for solving boundary value nonlinear ordinary differential equations, Int. J. Numerical Methods in Engineering 18:464-469 (1980). 13 R.P. Tewarson. A seventh-order numerical method for solving boundary value nonlinear ODE's, Int. J. Num. Methods in Engineering 18:1313-1319 (1980). 14 R.P. Tewarson and Vin Zhang. Solution of two-point boundary value problem using splines, Int. J. Numerical Methods in Engineering 23:707-710 (1986). 15 E. H. Twizell. Numerical methods for sixth-order boundary value problems, International Series for Numerical Mathematics 86:495-506 (1988). 16 E.H. Twizell and A. Boutayeb. Numerical methods for the solution of special and general sixth-order BVP's with applications to Ben£rd Layer eigen-problem, Proceedings of Royal Society London 431:433-450 (1990). 17 R. S. Varga, Matrix Iterative Analysis, Prentice Hall, Englewood Cliffs, N J, 1962. 18 F. Ayres, Schaum' s Outline of Theory and Problems of Matrices, McGraw-Hill, Singapore, 1982.