THE UNIVERSITY OF BRITISH COLUMBIA FACULTY OF GRADUATE STUDIES PROGRAMME OF THE FINAL ORAL EXAMINATION FOR THE DEGREE OF DOCTOR OF PHILOSOPHY Of RIAZ AHMAD USMANI 3J.Sc, The Ali g a r h Muslim University., India, 1956 '4i.Sc. The Al i g a r h Muslim University, India, 1957 THURSDAY, A p r i l 13, 1967 at 3°-30 p.m. In Room 102 Mathematics Building External Examiner: Dr. Ralph G. Stanton Department of Mathematics University of Manitoba. Research Supervisor: C. Froese. COMMITTEE IN CHARGE Chairman:: I. McT. Cowan C. Clark A . H. Cayford E. Leimanis' C. Froese J. Brenner M. MacMillan NUMERICAL SOLUTION OF BOUNDARY VALUE PROBLEMS FOR ORDINARY..DIFFERENTIAL EQUATIONS Abstract In the numerical solution of the two point "boundary value problem, y" = f ( x , y ) , y(a) = y &,y(b) = y^- 0 0 <a <b < • • .... {if the usual method i s to approximate the problem by a . f i n i t e ' difference analogue of the form .Vi yn+1 = h 2 \ P ± y , , n + ± > y 0 = y a ' yN+l = y b n = 0,1,...,N-1 . . . ( 2 ) with k = 2, and the truncation error T.E. = O(h^) or 0 ( h 6 ) , where h i s the step-size. Varga (1962) has obtained error bounds for the former when the problem ( l ) i s l i n e a r and of class M . In t h i s t h e s i s , more accurate f i n i t e difference methods are considered. These can be obtained i n e s s e n t i a l l y two d i f f e r e n t ways, either by increasing the value k i n difference equations ( 2 ) , or by introducing higher order derivatives. Several methods of both types have been derived. Also, i t i s shown how the i n i t i a l value problem y' = <i>(x,y) can be formulated as a two point boundary value problem and solved using.the l a t t e r approach. Error bounds have been derived f or a l l of these methods for l i n e a r problems of class M . In p a r t i c u l a r , more accurate bounds have been derived than those obtained by Varga (1962) and Aziz and Hubbard ( 1 9 6 4 ) . Some error estimates are suggested for the case where •2— < 0 , but these are not accurate bounds, *f e s p e c i a l l y when -2— i s not a constant. In the case of non-linear d i f f e r e n t i a l equations, s u f f i c i e n t conditions are derived for the convergence of the solution of the system of equations (2) by a generalized Newton's method. Some numerical r e s u l t s are included and the observed errors compared with t h e o r e t i c a l error bounds. -GRADUATE STUDIES Real Variable Theory S. Nash Point Set Topology S. Cleveland Theory of Functions R. Cleveland Functional Analysis C. W. Clark C e l e s t i a l Mechanics < E. Leimanis Numerical Analysis C. Froese. PUBLICATIONS "Certain I n i t i a l Value Problems Solved by Boundary Value Techniques," J.ACM 1 3 , 287-295 ( 1 9 6 6 ) . "Accurate Finite-Difference Approximat-ions of a Boundary Value Problem," Proceedings of the May 1966 Conference of the Computer Society of Canada. NUMERICAL SOLUTION OF BOUNDARY VALUE PROBLEMS IN ORDINARY DIFFERENTIAL EQUATIONS Riaz Ahmad Usmani B.Sc. (Hons.) , 1956; M.Sc, 1957, The Aligarh Muslim University, ALIGARH, India. A thesis submitted in partial fulfilment of the requirements for the degree of Doctor of Philosophy in the Department of Mathematics We accept this thesis as conforming to the required standard The University of British Columbia March, 1967 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an advanced degree a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x -t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by the Head o f my Depar tmen t o r by h i s r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n -c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Depar tmen t o f The U n i v e r s i t y o f B r i t i s h Co lumb ia Vancouve r 8, Canada Date ABSTRACT In the numerical solution of the two point boundary value problem, y" = f(x,y), y(a) = y &, y(b) = y b > -» <a <b < «, , ...(1). the usual method is to approximate the problem by a finite difference analogue of the form k 2 k $ a i yn+i = h I 3 i yn+i' y0 = V yN+l = V n = •"(2) 1=0 i=0 4 6 with k = 2, and the truncation error T.E. = 0(h ) or 0(h ), where h is the step-size. Varga (1962) has obtained error bounds for the former when the problem (1) is linear and of class M, i.e. f(x,y) is continuous, bounded and f > 0, but his method can be extended immediately to include y - 3 the latter case as well. In this thesis, more accurate finite difference methods are considered. These can be obtained in essentially two different ways, either by increasing the value of k in difference equations (A.E.) (2), or by introducing higher order derivatives. Methods with k = 4 are con-sidered as a one parameter family. When k > 2, additional A.E.'s are required for points near the boundaries. Such A.E.'s for k = 4 have been derived so that the overall resulting error in these methods turns out to be O(h^). These methods are shown to be"more accurate than a similar one proposed by Bramble and Hubbard (1964).: Difference' methods involving higher order derivatives of the type k k k • k I a. y . - h* I 3.y'*.+h 3 l y. y"' +.'..+ h C Y 9 , y ( ^ . . . ( 3 ) with k = 2, t = 3 and 4 are considered when (1) is linear. Now addi-tional equations are required to eliminate y' from the expressions for the higher derivatives obtained from (1) through over differentiation. Such equations have been determined so that the overall method has T.E. 8 10 = 0(h ) and 0(h ) for t = 3 and 4, respectively. Also, i t is shown how the i n i t i a l value problem y' = <j)(x,y) can be formulated as a two point boundary value problem and solved using these equations. Error bounds have been derived for a l l of these methods for linear problems of class M. In particular, more accurate bounds have been derived than those obtained by Varga (1962) and Aziz and Hubbard 9 f (1964). Some error estimates are suggested for the case where — < 0, 3 f but these are not accurate bounds, especially when is not a constant. In the case of non-linear differential equations, sufficient conditions are derived for the convergence of the solution of the system of equations (2) by a generalized Newton's method. Some numerical results are included and the observed errors compared with theoretical error bounds. The ratio of the. maximum abso-lute error and the theoretical bound was usually greater than 0.1, for the methods based on (3). Various methods are compared by treating the error as a function of cost, assuming the major cost is either (i) in the evaluations of the functions involved or (ii) in the solution of the system of equations. Under both assumptions, the method based on (3) with t = 4 gave the lowest error for a given cost over a wide range. It also compared favourably with well-known predictor-corrector methods for initial.value problems. iv TABLE OF CONTENTS CHAPTER PAGE I INTRODUCTION 1 II BOUNDARY VALUE PROBLEM OF CLASS M 9 III DIFFERENCE EQUATIONS OF HIGHER ORDERS 29 IV- A CLASS OF HIGH ACCURACY DIFFERENCE FORMULAS OF LOW ORDERS 49 V LINEAR BOUNDARY VALUE PROBLEMS WHERE 8 < 0. . . . 61 VI NON-LINEAR DIFFERENTIAL EQUATIONS . . . . 69 VII BOUNDARY VALUE TECHNIQUES AS APPLIED TO INITIAL VALUE PROBLEMS 79 VIIB EXPERIMENTAL RESULTS AND" CONCLUDING REMARKS . . . . . . 90 BIBLIOGRAPHY.' 117 LIST OF TABLES TABLE PAGE I EXPERIMENTS WITH PROBLEM (2.18). . . . . . . . . . ..•.-•'27 II DISTRIBUTION OF ROOTS OF G(s) = 0 for 2<s<3. . . . . . 34 III TABULATION OF THE FUNCTION -103 G(s) . . . 35 IV EXPERIMENTS WITH PROBLEM (5.10) 68 V . COEFFICIENTS OF THE A.E. (7.6) 82 VI OPTIMUM CHOICE OF c IN M- 3 ( c ) . . . . . . . 91 VII EXPERIMENTS WITH PROBLEM (A) WITH h = 1/8. . . . . . . . 92 VIII EXPERIMENTS WITH PROBLEM (D) WITH h = 1/8 103 IX EXPERIMENTS WITH PROBLEM (D) WITH h = 1/2 USING M-5 • AND M-6 104 X EXPERIMENTS WITH PROBLEM (B) . . 105 XI EXPERIMENTS WITH PROBLEM (E) 106 ; XII EXPERIMENTS WITH PROBLEM (8.2) 107 XIII EXPERIMENTS WITH NON-LINEAR PROBLEM (8.3), h = 0.1, u = 3/2 . . . . HO XIV EXPERIMENTS WITH PROBLEM (8.4) USING M-1, M-4 and M-5 ...-..•-. . - 111 . XV EXPERIMENTS WITH PROBLEM (8.4) USING M*-T AND M*-5 . . 113 XVI COMPARISON OF M-1, M-4 AND M-5 WITH STANDARD INITIAL, VALUE TECHNIQUES 116 vi LIST OF FIGURES FIGURE PAGE 1 Lo8;LO l E r r o r l versus log 2 (l/2h) for Problem (A) . . . . 94 2 L o 2io l E r r o r l versus log 2 Cost, Problem (A), Assumption 8.1(i) 95' 3 Log^Q |Error| versus log 2 Cost, Problem (A), Assumption 8.1(ii) . 96 4 L o g i o l E r r o r l versus log^ (l/2h), Problem (A), Single Precision Arithmetic 97 5 Log^ Q |Error| versus log 2 Cost, Problem (A), Assumption 8.1(i), Single Precision Arithmetic. . . . . 97 6 Log^ |Error| versus log 2 Cost, Problem (A), Assumption 8.1(i), Single Precision Arithmetic 98 7 L°Sio I E r r o r l versus log 2 (l/2h), Problem (D) 99 8 Lo§io l E r r o r l versus log 2 Cost, Problem (D), Assumption 8.1(i) . . . . . . . . 100 . 9 °^S^ o l E r r o r l versus log 2 Cost, Problem (D), Assumption 8.1 ( i i ) . 101 ACKNOWLEDGEMENT The author wishes to thank his advisor, Professor Charlotte Froese for guidance and helpful suggestions given during the prepara-tion of this thesis and Professor A.H. Cayford for reading and comment ing on the manuscript in the final form. The author also expresses his appreciation to the National Research Council of Canada for the financial assistance through its summer grants and the University of British Columbia Computing Centre for making available to him the facilities for the numerical experi-ments on their IBM 7040 computer. v i i i To my wife •i CHAPTER I INTRODUCTION 1.1 GENERAL CONSIDERATIONS The general solution of an nth order differential equation (D.E.) y ( n )(x) = F (x, y(x), y'(x),..., y ( n _ 1 ) ( x ) ) , for a real function y(x), normally depends on n parameters. These para-meters are determined in an i n i t i a l value problem by prescribing the values y ( m )(x) = A m (m = 0,1,..., n-1) at a fixed point x = a. If the conditions are based on more than one point, then the problem is called a boundary value problem. The boundary condi-tions usually have the form V. [y(x), y*(x),..., y ( n _ 1 ) ( x ) ] = 0 (i = 0,1,..., n-1) at the boundary points x = a and x = b which may include + ». Here y-fm^ (x) denotes the value of the mth derivative of y(x) at the point x. The func-tions F and V\ are in general nonlinear. In the present thesis, we shall deal with the special class of boundary value problems for which the D.E. is of second order (linear or nonlinear). 1.2 BOUNDARY VALUE' PROBLEM OF CLASS M A boundary value problem is said to be of class M, see [12], Chapter 7, i f i t is of the form y"(x> = f(x,y), y(a) = y &, y(b) = y b ...(1.1) where (i) - o o < a < b < ° ° (ii) y a and y^ are arbitrary constants ( i i i ) f(x,y) is a continuous function, of two variables with ^^gy continuous, bounded and non-negative in the strip S defined by a _< x <_ b and - 0 0 < y < °° . Although in what follows, we will also consider boundary value 3 f (x ) problems in which the above condition viz. — £ — > 0 , will be violated 9y -in S, we will concentrate mostly on problems of class M. The existence of a unique solution of the boundary value problem of class M is assured (Henrici [12, pp. 347]') . 1.3 NUMERICAL SOLUTION OF BOUNDARY VALUE PROBLEMS Numerically there are two methods for solving a boundary value . problem. The "shooting" technique is that of solving the i n i t i a l value problem y" = f(x,y), y(a) = y , y'(a) = c using some appropriate method. A correction procedure is then used to choose new values of y >'(a) such that y(b) converges to the desired value y^ This technique will not be considered in this thesis, even though there are problems where i t is a feasible procedure. The second method is the numerical approximation of the D.E.. (1.1) at a finite set of grid points, x n = a + nh (n = 0,1,..., N+1) where x Q = a, x ^ = b and h = (b-a)/(N+l) . Here N is an appropriate positive integer. A scheme is then designed for the determination, of the numbers y n , which i t is hoped, will closely approx imate the values y(x ) of the true solution of (1.1). A convenient way J n to obtain such a scheme is to demand that the y n satisfy a difference equation (A.E.) of the form I a. y , . = h 2 I p. y" (n = 0,1 x=0 x=0 > • • • j N-k+1) ...(1.2) Here y^ = £U±,y±) , a k + 0 and | c t Q | + | B Q |' * 0. We normalize (1.2) by choosing = 1. The positive integer k is called the order of the A.E. This method always leads to a system of N-k+2 sim-ultaneous equations involving y., , y 9 ,..., y , possibly in a nonlinear fashion, tional equations are required. The purpose of this thesis is to formulate more accurate proced-ures than the known numerical techniques for the solution of the boundary value problems and to obtain accurate error bounds to assist one in the selection of the step size h. It is aimed to show how such error bounds may be obtained for certain finite difference analogues approximating a class of problems of the form (1.1). In the case of linear boundary value problems, this purpose has been achieved in two w a y s , (a) .By increasing the order k of the A.E. from k = 2 to k = 4. In Chapter II we approximate the linear D.E. by a well known A.E. of order k = 2 of the 4 form (1.2) for which the resulting error is 0(h ). But in Chapter III, a general one parameter family of A.E.'s of order k = 4 is developed (most of the known A.E.'s of order k = 4 become special cases of this A.E.) and i t is proved that the resulting error is O(h^). (b) By using A.E.'s involving higher derivatives. The second important aspect is to modify A.E. (1.2) to the following form and N+l being determined by the boundary conditions. For k>2, addi-ii i n+i : i - 0 X k t>3 • (1.3) 4. and to use i t in approximating the linear D.E. (1.1). These ideas are -expanded in Chapter IV in connection with linear two point boundary value problems of class M. The linear boundary problems not of class M are dealt with in Chapter V using the same ideas. In Chapter VI, we consider nonlinear bound-ary value problems of class M. The iterative procedure used is merely a generalization of Newton's method for solving a system of nonlinear algebraic equations. Finally in Chapter VII, we show that the boundary value techniques studied in Chapters II and IV can be suitably applied to solve a class of i n i t i a l value problems of the' form y' = f(x) y.+ g(x), y(a) = y . C*. The last chapter is devoted to numerical results and conclusions. 1.4 Before we start a systematic discussion of linear boundary value problems in the next chapter, we would like to introduce some terms employed frequently in the thesis. We assume that the true solution y(x) of the problem (1.1) satis-fies a A.E. of the form k k T a. y(x ..) = h I 3. y"(x ,.) + T , ., ...(1.4) . L n l J n+i >_ i ' n+i n+.k ' i=0 i=0 where T n +^ is the truncation error (T.E.). Here y"(x^) = f(x^, y(x^)), and i f we introduce the discretization error £n = yiXr? ~ yn ' we obtain the error equation by subtracting (1.2) from (1.4) namely k 2 k y a. e , . = h g. a,. e + T . , ' ...(1.5) . L n i n+i .L_ i n+i n+i . n+k i=0 i=0 5. where a. e. = f(x., y(x.)) - f(x., v.), so that a. will be some value of x i x x x y x x \ ?^ by "mean value theorem' in two variables. 8y In actual practice, in any machine calculations, equation (1.2) will not be solvable exactly because of round-off and other errors which may be introduced in the evaluation of f(x,y). The round-off error depends primarily on the number of digits being carried out during computation but partly on the way in which the calculation is organized. However in this thesis, we want to concentrate on the accuracy of the difference approxi-. mation itself and therefore a l l such errors will be neglected. Difference Operator We associate with the A.E. (1.2) the operator k k . L[y(x); h] = I a. y(x +. ih) - h z £ g y"(x + ih) ...(1.6) 1=0 i=0 This operator acts on any function y(x) which has a second order derivative, but we assume that y(x) has continuous derivatives of sufficiently high orders. We may then expand (1.6) using Taylor's formula in powers of h and obtain, on collecting terms containing like powers 'of h, L[y(x); h] = C Q y(x) +.C± h y'(x) + C 2 h 2 y"(x) + ...,_ ... (1.7) where k : 0 i=0 1 k . ' • C = I i a i=0. , k • k C2 = 2^ I 1 " i " I 3 i ...(1.8) Z 1 ' .1=0 i=0 1 and C q =TT J Q ± q " i " 7Tb)T j 0 ^ \ <<-3.«.-> ' Now we define the degree of the A.E. as the unique integer p such that Cq = 0 (q=0,l,...,p+l); C p + 2 * 0. Then by definition we have L[y(x); h] = C p + 2 h P + 2 y ( p + 2 ) (x) + 0 ( h P + 3 ) , ...(1.9) We also note [using (1.3)] that the T.E. is given by However, for a number of difference operators, i t turns out that we can write L[y(x); h] = C p + 2 h P + 2 y ( p + 2 ) (O ...(1.11) where £ is a suitable number in the interval (x, x+kh). We shall refer to (1.11) as the generalized mean value theorem. This theorem does not hold for a l l operators, but every operator of degree p can be represented in the form k i L[y(x); h] = h P + 2 / G(s) y ( p + 2 ) (x + hs) ds , ...(1.12) 0 by using Taylor's formula with the integral representation of the remainder in the expansion of (1.6). The kernel G(s) is given as follows: G^ " \ (k(; f i y . - gk V-1);" » S £ K - 1 - V and for s e [u,u +1] (u = 0,1,..., k-2) G(s) = , [a k ( k - s ) P + 1 + a k_ 1 ( k - l - s ) P + 1 + ... + a j + 1(u+l-s) P + 1] " T ^ T T l [ ek ( k - s ) P ~ 1 + e k - i <k-i-s) p" 1 + . . . + B y + 1(y+i-s) p- 1] . An operator L[y(x); h] is called definite i f the kernel G(s) does not change sign in [0,k]. For a definite operator with y^ P + 2^(x) continuous, the gen-eralized mean value theorem holds since by the "second law of mean" (1.12) can be written as (1.11). However, i f G(s) does change sign, we s t i l l have |L [y(x); h]| < h P + 2 G M k where G = / |G(s)| ds , 0 p+2 (1.13) and M _ = max. P x<C<x+kh y ( p + 2> CO I. Associated Polynomials Sometimes i t is convenient to associate with the A.E. (1.2), the polynomials . k a(x) = £ a. x with a, = 1. 1=0 1 k (1.14) B(x) = I 3. x 1 i=0 It should be noted that i f the A.E. (1.2) has a positive degree p^l, then we have CQ = C x = C 2 = 0 . . ...(1.15) In terms of the polynomials a(x) and 8(x), the above conditions may be expressed in the form a(l) =0 a'(1) = 0 ...(1.16) _a"(l) = 2 8(1) The first two conditions of (1.16) imply that x = 1 is a root of a(x) = 0 of at least multiplicity 2. Thus we can write a(x) in the form a(x) = (x - l ) 2 y(x) ...(1.17) where y O O is a monic polynomial of degree k-2. Vector Norm Let v be a vector with components v^ ,' V2»..., v n . The vector norm ||v|| employed is defined by 8. | | v | | = max. | v ± | . ... ' . . . ( 1 . 1 8 ) i T h i s a l s o i n d u c e s a n o r m o n t h e m a t r i c e s a s f o l l o w s : ||A|I = max. ||Av|| w i t h | | v | | = 1 , w h e r e A = ( a — ) . I t t u r n s o u t t h a t t h e m a t r i x n o r m | |A| | s u b o r d i n a t e t o t h e v e c t o r n o r m d e f i n e d b y ( 1 . 1 8 ) i s g i v e n b y n | |A| | = max. I |a | . . . . ( 1 . 1 9 ) . i j = l 1 3 F o r a d e t a i l e d d i s c u s s i o n o f t h e n o r m s o n v e c t o r s a n d m a t r i c e s , t h e r e a d e r i s r e f e r r e d t o F a d d e e v a [ 7 ] . CHAPTER II BOUNDARY VALUE PROBLEMS OF CLASS M 2.1 We consider the boundary value problem of class M viz. y" = f(x,y), y(a) = y&t y(b) = y b ...(2.1) with 3 f ^ y y ^ — 0 a n d replace i t by a A.E. of the form (1.2). The coeffici-ents a's and $'s are usually chosen in such a way that the A.E. (1.2) turns out to be of as high degree as possible. The truncation error of the A.E. is thus T.E. - C p + 2 h P + 2 y ( p + 2 ) + 0 ( h P + 3 ) , C p + 2 * 0 . ...(2.2) The resulting system of algebraic equations for y n is implicit in any case* however another difficulty may arise. As in any system of equations, we need as many equations for the determination of the unknowns as there are unknowns. Since y^ and y^^ are determined by the boundary conditions, the unknowns in our case will be y , y2,..., y^ . If the order k of the A.E. is greater than 2, the new unknowns such as y_^ or y^+£ a r e introduced for which there are no equations. This difficulty is- overcome by suitably mod-ifying the A.E. near the boundaries; however i t does not arise at a l l i f k=2, the smallest possible value of k. In view of relations (1.16) i f k=2 and the A.E. has a positive degree p, then i t is necessarily of the form ^ " 2 y n + l + ^ 2 - h ' ( e0 y"n + h + 32 ^n-Kp ' •' ( 2 ' 3 ) with 6 Q + &1 + B 2 = 1 . Now we shall apply representation (1.12) in order to prove the generalized mean value theorem (1.11) for the operator associated with the_—-10. A.E. (2.3). Expanding around x, using Taylor's formula with integral remainder, we get L[y(x); h] = y(x) - 2 y(x+h) + y(x+2h) - [S y"(x) + &± y"(x+h) + 32y"(x+2h)] = hP + 2 / G(s) y 0 where G(s) = p+1 (p+2) (2-s) (x+hs) ds p-1 (p+D: "2 (p-i): , 1 _< s < 2 (2-s) p+1 (1-s) p+1 (P+D : (P+D : (2-s) P-I (l-s) p-1 'i (P-D: P I ( P -D : x , , o < S < I For p=l, 3Q + 3 1 +• 8 2 = 1 , G(s) = \ ((2-s) 2. - 2 S 2) , 1 < s < 2 \ (2 BQ - s 2) , 0 < s < 1 It i s easily shown that G(s) changes sign i f (D o < e 2 < \ or i f o < e Q < -2-( i i ) both 3Q,, B 2 >. \ ( i i i ) 3Q , 3 2 < 0 , and G(s) is of constant sign i f 3 2 >. ^ , 3Q £ 0 o r ^ - 2' 2^ - 0 For p=2, we have by (1.8) 3Q +^ -,^ + 32 = 1, 3 X + 23 2 = 1 The above conditions together imply that 3Q =32 so that G(s) = YT (2-s) 3 - 3 2 (2-s) , 1 <. s < 2 G(2-s) , 0 <. s _< 1 i.e. the function G(s) is symmetric about the line s=l. 11 Also G(s) 1 0 , s 2 = 0 For p=3 , by (1.8), we have + 3]_ + S 2 = 1 , + 23 2 = 1 , J l + 4 S2 - £ ' The solution of this system is given by 3Q = 3 2 = 1/12 , 3 2 = 10/12 , and the corresponding coefficient vanishes. Thus there is no operator of degree p=3 associated with the A.E. (2.3). However = -1/240, hence the operator associated with the A.E. (2.3) with 3Q=32=:1/12,3-^ =10/12 turns out to be of degree p=4 and the corresponding kernel G(s)<0 in [0,2]. The question of p being greater than 4 does not arise at a l l , since a l l the S's are determined already with p=4. Thus the operator associated with the A.E. (2.3) is definite, provided p=l, 2 or 4. In case'p=l, either 3Q£.0 or 3^1/2, 32£0 and i f p=2, then either 32=0 or i f 3 2^ 0, then 3 2^l/6. 2.2 A.E. OF ORDER 2 ' If we make use of A.E. (2.3) in approximating the D.E. (2.1) set-ting n=0,l,..., N-1 in (2.3), we get the following system of equations Jy + h 2 D f(y) = d , (2.4) where J = 2 -1 -1 2 -1 -1 2 -1 -1 2 " -1 1 50. 2 51 h l0 h >0 p l p2 3 3, 0 1 y = 'N f(y) = f ( x l ' y l ) _ f (x2,y2-) f(x 3,y 3) and d = y a-3 Qh f(a,y a) 0 0 0 y b-3 2h 2f(b,y b) Here we are concerned with the case where the given D.E. (2.1) is linear i.e. f(x,y) is of the form f(x,y) = g(x) y + s(x) , g(x) > 0 in S. Defining a new diagonal matrix G G = and the vector! s = N where g^ = g(x^) and s^ '= s(x^), we may write the vector f(y) introduced in (2.4) as f (y) = G y + s . The algebraic system (2.4) now reduces to the system of linear equations Ay - b ...(2. where A = J + h 2 D G and b = d - h 2 D s . The numerical solution of the system (2.5) is relatively easy owing to the fact that the matrix A shares with the matrices J and D the property of having non-zero elements only on the main diagonal and the two diagonals adjacent to i t . Matrices of this kind are called tridiagonal. . In fact, i f A = ( a ^ j ) j w e h a v e b y v i r t u e o f ( 2 . 5 ) a = 2 + h 2 3 1 g n ( n = 1, 2,..., N) n , n a n d a n , n - l = " 1 + h B 0 § n - 1 ( n = 2> 3 " - - > N ) a n , n + l " -1 + h ' B 2 § n + l ( n = 1 ' 2,..., N-1) a = 0 f o r I n - ml >1 . L i n e a r s y s t e m s i n v o l v i n g n o n - s i n g u l a r t r i d i a g o n a l m a t r i c e s a r e m o s t e a s i l y s o l v e d b y 1 a m o d i f i c a t i o n o f t h e G a u s s i a n e l i m i n a t i o n a l g o r i t h m . We l e a v e t h e d e t a i l s o f t h e a l g o r i t h m f o r b r e v i t y [ 1 2 , p p . 3 7 5 ] , We m i g h t n o t e t h a t t h e w h o l e p r o c e s s o f s o l v i n g t h e s y s t e m ( 2 . 5 ) r e q u i r e s a p p r o x i m a t e l y 3N a d d i t i o n s , 3N m u l t i p l i c a t i o n s a n d 2N d i v i s i o n s . 2.3 THE D I S C R E T I Z A T I O N ERROR I N L I N E A R PROBLEM -Now we s h a l l e s t a b l i s h b o u n d s f o r t h e d i s c r e t i z a t i o n e r r o r . We n o t e t h a t t h e t r u e s o l u t i o n y ( x ) o f t h e l i n e a r p r o b l e m s a t i s f i e s t h e s y s t e m A y ( x ) = b + T ( y ) . . . ( 2 . 6 ) w h e r e ' , 5 r ( p + 2 ) r(y) - c p + 2 h p+2 y ( p + 2 ) ( X l + 9 l h ) ( x 2 + 0 2 h ) H e r e IQ^lji x , i t b e i n g a s s u m e d t h a t t h e o p e r a t o r a s s o c i a t e d w i t h t h e A . E . ( 2 . 3 ) i s d e f i n i t e s o t h a t t h e g e n e r a l i z e d m ean v a l u e t h e o r e m h o l d s . Now s u b t r a c t i n g ( 2 . 5 ) f r o m ( 2 . 6 ) , we g e t t h e e r r o r e q u a t i o n Ae = T ( y ) ' . . . ( 2 . 7 ) w h e r e e d e n o t e s t h e e r r o r v e c t o r . T h e m a t r i x A h a s t h e f o l l o w i n g p r o p e r t i e s p r o v i d e d 3^, >. 0 a n d n 8^ §^1 (i = 0, 1, 2), where g = max. g(x) a<x<b (i) A is real, irreducible matrix (ii) A is monotone matrix i . e . A _ 1 > 0 (2.8) and 0 < A - 1 < j f 1 J ^ is known, in fact i f J ^ = (j ), then ' Jmn 'mn .n(N+l-m) N+1 m(N+l-n) n < m N+1 , n _> m , See Rutherford [23]. The truth of the above properties of the matrix A follows from the following well known theorems: Theorem 2.1 A tridiagonal matrix A = (a.^) is irreducible i f (i) a . ^ 0 (ii) a i > i + 1 ^ 0 (i = 2,3,..., n) (i = 1,2,..., n-1) Theorem 2.2 Let the matrix A = (a..) be irreducible and satisfy the conditions (i) a.. < 0 n (ii) I j = 1 Then A is monotone. Theorem 2.3 a. . i ^ j , (i , j = 1,2,..., n) > 0, i = 1,2,...,n > 0, for at least one i Let the matrices A and B be monotone, and assume that A-B >_ 0, then -1 - I B - A > 0 . • For the proof of the above theorems, see Henrici [12, pp. 360]. Now we consider the error equation (2.7) and rewrite i t in the form e = A - 1 T(y) . • From the above equation, we deduce — |e | < |C /_| h P + 2 M . , f j ' m' — ' p+2' p+2 L_ Jmn , n=l using (2.8 (ii)) , < ,!C • I hP+2 M x m(N+l-m) - ' V 2 1 p+2 2 (x -a) (b-x ) < I h" ( V a ) .(b-xm) , (m = 1,2,... , N) where M - |C p + 2| M p + 2 l Also . • ||e|| =' max. |e | = M h P (b-a)2/8 l<m<N We summarize these results in the following theorem. Theorem 2.4 Let y(x) be the true solution of the linear boundary value prob-lem of class M and y be its discrete approximation with an error e . Let n r n the error equation be Ae = T(y) , (see (2.7)). where A is a monotone matrix such that A - J > 0 and ||r(y)|| <_. | C p + 21 h P + 2 M p + 2 , 'then 16. |e±| <Mh P (x±-a) (b-x±)/2 , (i = 1,2,...,N) and | |e| |- ' < M h P (b-a)2/8 , where M = |C p + 2| M p + 2 , or equivalently ||e|| = 0(hP) . A special case of Theorem 2.4 was proved by Varga [25], see Theorem 6.2, pp. 165, in which the approximation of the linear D.E. is based on the A.E. " yn + 2 y n + l - yn +2 + h ' y " n + l = 0 ' • • • <2--10a> with p = 2 and M = | r M, . Thus 12 4 I l e l I £ h 2 M 4 (b-a)2/96 . A more accurate approximation is based on the A.E. " yn + 2 y n + l - yn +2 + ^ 2 / 1 2> <y"n + 1 0 y " n + l + ^'n^ = 0 '' ' ( 2 ' 1 0 b ) with.p = 4 and M = TTTTT M, and thus 240 6 l l e l l 1 h 4 M 6 (b-a)2/l,920 . ,..(2.11) The difference operators associated with the above A.E.'s have been proved to be definite and hence the generalized mean value theorem (1.11) holds,(Section 2.1). 2.4 MORE PRECISE ERROR BOUNDS The error bounds obtained in Theorem 2.4 depends on the fact that 0 £ A _ 1 £ J _ 1 . and thus on r: T ,„, „ N , n AT N(N+2)/8 , N even IA"1!! < UJ"1! (N+l)2/8 , N odd < (N+l)2/8 = (b-a)2/8 h 2 . 17. However, in view of Theorem 2.3 and properties of the matrix A , given by (2.8), j | A "'"I I monotonically increases as g(x) decreases and attains i t s maximum value (b-a) /8 h when g(x) = 0. Setting A^ = matrix A in whi ch each element g. has been replaced by g„, = max. g(x) and A = matrix A in a<_x<b which each g. has been replaced by g = min. g(x), then from the theory of m a<x<b montone matrices, i t follows that I A " 1 ! ! < \ \A^\\< M J " 1 ! (2.12) , - l i equality holds i f f g(x) = 0 . In order to obtain | |A M || which would enable us to get a better error bound than already obtained, we consider -2 + h 2 3 L G M - l + h 2 B 2 g m - l + h 2 & o g m 2 + h 2 - l + h 2 8 2 g m A = m -1 + h 2 g g 2 + h 2 0 g -1 + h 2 3 9 g m 0 m 1 m z. m -1 + h '0 *m 2 + h h Sm with 8 = 8 and 8 + 8-, + 8 9 = 1 and this w i l l be the case for p = 2 or A. Then m . u -1 -1 v -1 -1 y -1 = P N (u) Q , where X = 1 and y = (2 + h 2 8 X 8 m ) / ( l - h 2 3 Q g m) Therefore 18. -1 -1 where Q We now proceed to determine exact expressions for P^(y) and Let D denote the determinant of the matrix P (p) , then i t is easy to show using Laplace's method of expanding determinants that sat-isfies a A.E. with constant coefficients viz. (2.13a). D = y D - D (n = 1,2,...) n n-l n- z. ... .(2.13a) D _ 1 = 0, D Q = 1 - ; Please note that for a given problem and for a fixed h, u is constant. On solving the above A.E. with i n i t i a l conditions = 0, DQ = 1, we obtain D = n "sinh (n+1) 9 sinh 9 (n+1) sin (n+1) 9 sin 9 , i f u = 2 cosh 9 > 2 , i f y = 2 , i f u = 2 cos 9 < 2 (2.13b) We first determine the inverse of the matrix P (u) f° r smaller value of say n = 1,2,3 and 4. This can readily be done using the Gauss elimination method or the adjoint method for finding the inverse of a matrix i f n is not very large. The p "'"(v) f° r n = 1.2,3 and 4 is given below: p 2 1 ( ^ " - k V i DoDo DoDo D1 D0 19. and i _ D, P^(y) i_ D. D0 D2 V i Do Do D o D i D1 D1 D1 D0 _DoDo D i D o D2 D0 V 3 D 0 D 2 D o D i V o D0D2- D1 D2 D1 D1 D2°0 D o D i °1D1 D2°l V o V o D2°0 D2 D0 D3 D0 The form of the above inverse matrices suggests that in general P > ) = ( P i j) ...(2.13c) where p. . = p. . i J 31 D. D ./D 3 < ± 2 - 1 n-i n ' J — We will prove the above result (2.13c) using induction argument below. We assume that the above result is true for the matrix-P -, (y) , n-1 i.e. i f P-1., (y) = (p..), then p.. = p.. = D' D . . /D n , j < i . Let n-1 i j 13 J i j-1 n - l - i n-1 ' J — the sequence of the matrices P.. (y) , P„(y),...P (y) be written as follows: Px(y) -1 P2(y) 0" P1(y) = [y], P2(y) = , P3(v) = -1 -1 • • H. _0 •' -1 y_ Pn(y) = Pn-1^> 0 0...0 -1 0 - I y In general, we partition the matrix Pn(p) into 4 sub-matrices which we write as Pn(y) - '11 21 412 22 H e r e a . , = P 11 n - 1 ( p ) , a 2 1 = [0 0 . . . 0 - 1 ] , ( l x n - 1 ) row m a t r i x , a ^ , a ( n - l x l ) c o l u m n m a t r i x a n d e q u a l s t h e t r a n s p o s e o f a n d f i n a l l y a 2 2 I f we w r i t e 11 21 12 !22 p a r t i t o n e d e x a c t l y i n t h e same f o r m a s P (p) i s , t h e n we h a v e '11 '21 12 22 11 '21 12 l22 = I n w h e r e I i s t h e u n i t m a t r i x o f o r d e r n . On m u l t i p l y i n g t h e 2 m a t r i c e s o n t h e l e f t , we g e t f o u r m a t r i x e q u a t i o n s n a m e l y 11 312 + a12 6 22 = 0 21 hi + a22 8 2 2 = 1 11 hi + a12 . B 21 = n - 1 21 hi + a 22 hi = 0 . S o l v i n g t h e s e m a t r i x e q u a t i o n s f o r 8^ ^^ , B ^ 2 , 8 2 ^ a n d 8 2 2 > we h a v e B 2 2 = D 8 2 1 = -D ( a 2 1 8 1 2 = -D ( a " J a 1 2 ) 8 n = + (all o 1 2 ) D ( a ^ . o ^ ) = al]_ + ( a 1 2 a 2 1 ) / D w h e r e D = ( a 2 2 - a^) On s u b s t i t u t i n g = ^ n ^ 2 . ^ ^ a n ^ t * i e v a ± u e s o r t h e ' m a t r i c e s ' a ^ 2 > a2i a n c ^ a 2 2 i n t h e a b o v e r e l a t i o n s , we g e t , 322 = D ,/D n - 1 n 521 ' (Do> V ' D n - 2 ) / D n 5-^ 2 b e i n g e q u a l t o t h e t r a n s p o s e o f 82^ 21. and B l l = P n - l ( ^ + 6 1 2 e21 / D = _ < Y l D n - l - i / D n - l > + D D ~ ^ - 1 Di-1> n n-1 J i D D + D 1 /T, n n-l - i l - l . = F~ ( D j - i " D ~ ) ' 3 — x n J n-1 = ( D D . / D ) , j < i j-1 n-i n J — ' using the identity (Dn V l - i + Di-1 ) / Dn-1 = °n-i (i - 0,1 ,2 . . . , n) ...(2.13d) Proof of the Identity We will prove the identity for u ^ 2 using (2.13a). Consider the A.E. (2.13a) in the form °r = y Dr-1 " °r-2 ( r = i. 2' 3'---) D - 1 = ° ' D 0 = 1 We now consider the above A.E. for r = n and r = n-i, and eliminate y between these two relations and write the result in the form ( * ) D D . , - D . D . = D . D . 0 - D „ D . . . n n-i-1 n-1 n-i n-1 n-^ i-2 n-2 n-i-1 We then consider the above A.E. for r = n-1 and r = n-i-1 and prove ( * * ) D D . 0 - D 0 D . T = D 0 D - .' - D _ D . 0 . n-1 n-i-2 n-2 n-i-1 n-2 n-3-i n-3 n-i-2 Next we consider the A.E. for r = n-2 and r = n-i-2;,'... ; and finally r = i+1 and r = 1. Combining the results of the form (*) and (**), we get ( * * * ) D D . N - D . D . = Constant = X (say), (r = n, n-1,..., i) . r r-i-1 r-1 r - i J ' In order to evaluate X, we put r = i in (***) and thus X - - D i _ ^ • Hence -D D . . - D . D . = X = - D . ' (i = 0,1,2,..., n) , n n-i-1 n-1 n-i i - l 22. or (D D . . - D. _)/D = D . , n n-i-1 i - l n-1 n-i and this completes the proof of identity (2.13d). Now substituting 3-j_i' ^ 2.2' ^ 21 & n C* 2^2 W £ ^ i n a- 1-ly obtain ^ =(Pij> where p . . = p.. = (D. , D ./D) , j < i . . . r i j r j i j-l n-i n ' J — and this completes the proof of (2.13c). For the determination of | |P ^(y)|| , we will require two "more identities which we give below: n-1 / JQ °j = T^ 2 C°i " °i-l " 1 ] ' ( y * 2 ) * ...(2.13e) Proof of Identity We'write (2.13a) for n = 1,2,..., i and then summing up a l l the i equations columnwisewe get ' D l = * D 0 " °-l : D2 =' y Dx - DQ D3 = ^ °2 - °1 • D. = y D. . - D. „ I i - l i-2 i - l i - l i - l -D.+D +( I D ) = y I D - ( I D - D ) ° 1 j=0 2 j=0 3 j=0 3 i - l (y-2) I D. = (D. - D._1 - 1) i - l h 2 ~ V-2 L " i " i - l or I D. = 73- [D. - D._n - 1] , (y/2) . 3=0 23. D . D . - D . D . , = D n-i i i - l n - i - 1 n .. (2.13f) The above identity can also be proved using the A.E. (2.13a) and following the technique employed for proving (2.13d) We now state the theorem giving | |P "'"(u) | | below: Theorem 2.5 IP" 1^)!! < n — '[!•- sech (£±i ©)]/(p-2). , i f u > 2 (n+1) /8, i f p = 2 n sin &. sin (n+1) 9 , i f P < .2 equality holds only i f n is odd positive integer and p >_ 2 Proof For p > 2 n i n y | P..i = ^ - t y D . . D . + y • D . D . . ] . ' • ' u i 1 D • i 3 - 1 n-i . .., n-3 i - l 3 = 1 J . n 3 = 1 J 3 = 1 + 1 - [ D . y D . . + D . . y D n-i 3-I i - l . h,.. n -3 ' n 3 = 1 . J 3=1+1 J D _ . D [ D , - D - 1] + T H T I T [ D (p-2) D i i - l n (p-2) D L n-i n-i-1 n using (2.13e) , 1] , D . D . - D . . D . , D . + D . _ 1 r l n-i i - l n-i-1 n-i i - l p-2 1 V-2 D n n D _ . + D . _ [1 - — ] , using-(2.-13f) n Since P^^Xp) is symmetric and also is such that i n j=l 13 V l - i . j . » ' n therefore, row sums increase with i , hence £ p.. is maximum for 3=1 1 J 1 = n/2 , n even — — , n odd thus, P-2 (1 -y-2 1 ± ^ 2 D + D ,„ n/2 n/2 - 1 even 2 D (1 - n-1/2 D ), n odd ,n+l [1 - sech ( — 9)] , y> 2 , (equality holds for n odd p o s i t i v e integer).-The case for y = 2 has been given in.'the beginning of the Section 2.4. For y < 2, D = sin(n+1)6/sinQ } j=l n a n d I IP,- J = I D. D , 3-1 n - i D n n + I j=i+l D . D. /D n-j i - l n = I sin(,j9) s i n ( n - i + l ) 9 sinG sin(n+l)9 n + i j=i+l sin(n-j+l)9 s i n ( i 9 ) sin9 sin(n+l)9 < ± [ i + ( n - i ) ] • |sin9 sin(n+1)9| . .< n/|sin9 sin(n+1)9| (an expression independent of i ) . < n/|sin9 sin(n+l)9|, y .< 2 . This completes the proof of the theorem 2.5. The above r e s u l t s enable us( to write down the norm of another matrix P ( y ) , c l o s e l y r e l a t e d to the matrix P ( y ) . The norm of the matrix n n P ^ ( y ) defined below w i l l be required l a t e r i n Section 3.3. Thus I|P 1 ( y ) ' ' n 25. P n(y) -U 1 1 y ' l = [-P n(-u)] ..(2.14) Let D~ denote the determinant of the matrix P (y) n n v minant then leads' to the following A.E. 5 n = * V l " 6n-2 » 5 - l = 0 > 5 0 = 1 from which immediately follows (see (2.13a))that Expansion of the deter-Dn(u). = D n(|y|) The inverse of the matrix P (y) i s known e x p l i c i t l y . I f P (y) = (p..) , n n • i j p.. = p.. = ( - 1 ) 1 + 3 D. . D ./D , j < i * i j * j i j-1 n - i n ' J — then see Rutherford (1952), [23]. The elements of the matrix Pn^"(y) can also be deduced using (2.13c). Now since D (y) = D (|y ), and D i s an even or odd function of y according as n n n ' 1 n ° i s an even or odd p o s i t i v e integer, we notice that the elements of the i n -verse matrices P ^~(u) and P "*"(y) are the same i n magnitude, but d i f f e r i n n n sign only. Hence ..(2.15) In our present discussion y = (2 + h 2 B 1 g m ) / (1 - h 2 B 0 g m) > 2 . 2 This follows from the fact that 3 > 0 , g(x) > 0 and h 3 Q g m < 1 . The-2 l a t t e r condition v i z . h 3Q g m < T," i s an immediate consequence of the assumption h^ $. g M < 1 ( i = 0,i;2) 26 which we made i n proving the i r r e d u c i b i l i t y of the matrix A (see 2 . 8 ( i ) ) . Thus f i n a l l y , f o r y> 2 I I A T 1 1 I £ H N ( y ) / h 2 g m . ...(2.16) where N + 1 Hjjdi) = [1 ~ sech ( — 9] , 2 cosh 9 = y (N r e f e r s to the order of the associated matrix) , - [1 -:sech ( ^ l o g £ V + ' f * - * ) ] . Thus again using the e>ror equation (2.7) we get for the l i n e a r boundary value problem of class M the more accurate error bound | |e| | < M h P V ^ m ' ...(2.17) Now i n view of the above, we restate the Theorem 2.4 as follows: Theorem 2.6 p+2 Let y(x) e C be the true s o l u t i o n of the l i n e a r boundary value problem of class M (given by (2.1) for f(x,y) = g(x)y + s(x)) and l e t y^ be i t s d i s c r e t e approximation based on A.E. (2.3). Furthermore, l e t 2 min. g(x) = g > 0 and h be such that h 8. g < 1 ( i = 0,1,2), then a<x<p | |e| | < M h P V y ) / g m ° r • ||e|| = 0(h p) . P a r t i c u l a r Cases of Theorem 2.6 (i) Using p = 2, B Q = B 2 = 0, B = 1, M = M /12, we get ||e|| < h 2 M 4 ^Cv)/12 g m , y = 2 + h 2 ^ ...(2.18) ( i i ) Using p = 4, B Q = B £ = 1/12, $ ± = 10/12 and M = Mg/240, we get ' , , , , 4 2+(10/12)h 2 g Me|| 1 h 4 ^ H N ^ / ^ Q g m ^ = 1 _ ( 1 / 1 2 ) h 2 g m m •••( 2- 1 9> The error bound given by (2.19) i s better than the one given by (2.11) as i s displayed by the following Table I. For experimentation, we chose the problem y" - y = x 2 - 2, y(0) = 0, y(b) = (2 sinh b/sinh 1) - b 2 ...(2 2 with y(x) = (2 sinh x/sinh 1) - x for a serie s of values of b and h res p e c t i v e l y , (Table I ) . TABLE I EXPERIMENTS. WITH PROBLEM (2.20) b h Error bound E 1 Error bound E l ^ E 2 based on (2.11) based on (2.19) 1/4 0.4069 E-5** 0.3684 E-5 1.10 1/8 0.2543 E-6 0.2303 E-6 . 1.10 1/16 0.1589 E-7 0.1439 E-7 1.10 1/32 0.9934 E-9 0.8995 E-9 1.10 1/4 0.6510 E-4 0.2390 E-4 2.72 1/8 0.4069 E-5 0.1494 E-5 • 2.72 1/16 0.2543 E-6 0.9336 E-7 2.72 1/32 0.1589 E-7 0.5835 E-8 2.72 1/4 0.1042 E-2 0.3253 E-4 32.02 1/8 0.6510 E-4 0.2033 E-5 32.02 1/16 0.4069 E-5 0.1271 E-6 32.02 1/32 0.2543 E-6 0.7942 E-8 32.02 ** We write 0.4069 E-5 for 0.4069 X..10 . We w i l l adopt t h i s convention throughout t h i s t h e s i s . 28. An 0(h ) three point f i n i t e d i f f e r e n c e analogue was also formu-lated by Aziz and Hubbard (1964), [2], In fac t the problem considered by them was L y = -y" + g(x) y = s(x) 7(0)= y ( l ) = 0 . . (2.21) 4 6 where g(x) _> 0, g(x) , s(x) are of class C on [0,1] and y(x) e C . The main r e s u l t of t h e i r paper i s that the d i s c r e t i z a t i o n error s a t i s f i e s the i n e q u a l i t y 2 i -48 gM -1 7 h 4 M, ,4 §. + k_ 720 6 c(4) (x) M (2.22) The error bound given by (2.19) i s better-than the one given by e i t h e r (2.11) or (2.22). For instance i n the problem y" - y = x e X , y(0) = y ( l ) = 0 , the error bounds given by (2.11)r, (2.19) and (2.22) are 0.51 E-5, 0.46 E-5 and 26.0 E-5 r e s p e c t i v e l y with h = 0.1. In the l a s t chapter, we w i l l compare the accuracy of the method proposed by Aziz and Hubbard with other known techniques and with those developed i n the present t h e s i s . The numerical technique i n which the two point l i n e a r boundary value problem i s replaced by A.E. (2.10 ) with the error bound as given by (2.19) w i l l henceforth be designated f o r reference purposes as M-1. S i m i l a r l y , the method proposed by Aziz and Hubbard w i l l be r e f e r r e d to as M-AH. CHAPTER I I I DIFFERENCE EQUATIONS OF HIGHER ORDERS 3.1 In the l a s t chapter, we considered the use of A.E. (1.2) with k = 2 for the numerical s o l u t i o n of l i n e a r boundary value problem of 2 class M, We proved that the r e s u l t i n g error i n one case i s 0(h ), while 4 i n the other case i t i s 0(h ). The two d i f f e r e n t error bounds each of 4 0 ( h ) are given by (2.11) and (2.19) r e s p e c t i v e l y . Now i f we want to consider f i n i t e d ifference analogues i n which the r e s u l t i n g error i s 0(ht), t>4, then one way to do t h i s i s to consider (A.E.) (1.2) with k>2. Other, methods f o r obtaining accurate A.E.'s without increasing the order k w i l l be considered i n the next chapter. We have already seen i n Chapter I (1.17) that a A.E. of p o s i -t i v e degree has an. associated polynomial k i 2 a(x) = 1 a. x = (x-1) Y(x) , ...(3.1) i=0 1 where y(x) i s a monic polynomial of degree k-2. In fac t any choice of the polynomial Y(x) r e a d i l y determines a's. We then determine the c o e f f i c -ients B's from a set of l i n e a r equations C = 0 ( i = 2,3,. . . , k+2) , ...(3.2) see (1.8). The above system of l i n e a r equations s a t i s f i e d by B's can be wri t t e n i n matrix form as follows: 30. 1 1 1 2 0 . k-1 k 2 2 . (k-l) Z k Z (k-l) k k k r 1 1.2 1 2.3 V - 2 L 1 <\-(k+1)(k+2) .k+2 a a. ...(3.3) where a l l summations in (3.3) extend from i=0 to i=k. The matrix assoc-iated with the above system is a Vandermonde matrix, hence i t is non-sing-ular and the coefficients 3's are therefore uniquely determined. The determination ofB^'s for a fixed value of k can in general be effected in two ways. We can either invert this Vandermonde matrix, Parker [21], or we can use Gauss elimination method. In next section,, we will consider A.E. (1.2) for k = 4 and develop a one parameter family of A.E.'s to approximate the linear boundary value problem.of class M. 3.2 A.E. OF ORDER k = 4 2 In view of (3.1), a(x) = (x-1) y(x) where y(x) is now a monic polynomial of degree 2. Hence 2 2 a (x) = (x-1) (x + cx + d) = x 4 + (c-2) x 3 + (d-2c+l) x 2 + (c-2d) x + d , where c and d are real arbitrary constants. Thus 0 31. Any choice of c and d determines a's.and then we determine 3's by solv i n g the system (3.3) for k = 4, and T.E. using (1.10). An easy c a l c u l a t i o n shows that ( - c + 19 d - l)/240 (12 c + 102 d + 2)/120 (97 c + 7 d + 7)/120 ...(3.5) (12 c + 2 d + 102)/120 ( - c - , d + 19)/240 0 and T.E. d - l , 7 (7) , , . 240 h y , d * 1 31c-190 ,8 (8) . . , 190 607O0-h Y ' d = 1 ' ° *~31 ...(3.6) -79 ,10 (10) _ 190 n y , d=l, c = _585,900 " ' » " -» - 31 However, f or obtaining symmetric A.E.'s, we w i l l require that d=l. The 8 T.E. of the corresponding A.E. w i l l be 0(h ) and we might expect that the r e s u l t i n g error i n the D.E. w i l l be 0(h ). at most. Thus for d=l, the c o e f f i c i e n t s a's and B'S are given as. follows: a o = a 4 = 1 a l = a 3 = c - 2 a 2 = -2 c + 2 e o = e v = (-c + 18)/240 3 1 = S3 " (12 c + 104)/120 h - (97 c + 14)/120 that 3's are non-negative for ..(3.7) -14/97 ,< c < 18. Thus we a r r i v e at a one parameter family of symmetric 32. A.E.' s namely, " ( y n + W + ( 2 " C ) ( y n + l + y n + 3 } + ( 2 c " 2 ) W + h I 3 1 x=0 y " n + . = 0 ...(3.8) where c i s an a r b i t r a r y parameter and 3's are given by (3.7). The T.E. associated with (3.8) i s given by (3.6). Each equation (3.8) may be considered as being centred on x n+2 and represents an equation for y . • But because two a d d i t i o n a l points are n+z involved on e i t h e r s i d e , i t cannot be used as an equation for e i t h e r y^ or y^ and therefore two a d d i t i o n a l equations must be introduced. The A.E.'s used near the boundaries are given by (3.9) below. The reason for t h i s p a r t i c u l a r choice of A.E.'s (3.9) i s that the matrix J(c) defined below by (3.11) becomes r e a l , symmetric, i r r e d u c i b l e and monotone for any c>2 . „ 3 (i ) -c y Q + (2c-l) y± + (2-c) Y 2 " Y 3 + h . I B. y * ^ = i=0 2 3 ( i i ) -y* T_o + ( 2- c) y M _ i + ( 2 c - D y M - c + h I • i=0 0 = 0 ...(3.9) where B Q = c / 1 2 , g = ( 1 0 c + l ) / 1 2 , &2 = ( c + 1 0 ) / 1 2 and B 3 = 1 / 1 2 . The T.E. associated with (3.9) i s given by "-(c+1) .6 (6) , . 2 4 0 h y . . c ^ 1 T.E. = - 1 240 h7 c = -1 for A.E. (3.9 i ) 1 h 7 y ( 7 ) , c = -1 for A.E. ( 3 . 9 i i ) 240 (3.10) The usual d e t a i l s of the d e r i v a t i o n of equation (3.9) have been omitted for b r e v i t y , We notice that.the c o e f f i c i e n t s 3's are p o s i t i v e f or any; value of c>0 . Now we would l i k e to f i n d those values of the parameter c 3 3 . for which the A.E.'s (3.8) and (3.9) are such that the difference opera-tors associated with them become definite and the generalized mean value theorem (1.11) holds. In case of A.E. (3.8), we have L[y(x); h] = [y(x) + y(x+4h)] + (c-2) [y(x+ h) +. y(x+3h)] + (2-2c) y(x+2h) - h 2 FfJ^ p [y"(x) + y" (x+4h) ] + [y"(x+h) + y"(x+3h)] + ^gi4- y"(x+2h)} 30 = h G(s) y ( 8 ) (x+hs) ds where G(s) = "(4-s) _ -c+18 : (4-s) 7 ! 240 ' 5 I 7 ,„ x7 3<s<4 (4-s) + (c_2) (3-s) _ -c+18 . _ 3c+260 . (3-s)" 2 7 : + ( c Z ) i : 240 5 : 30 5 : - -and G(s) = G(4-s), 0£s<2 . i.e. the function G(s) is symmetric about the line s=2. Now for 3<s<4 provided c <_ 86/7 . Furthermore, G(2) = (96c - 448)/(5 x 8 !) <. 0 , provided c <^ 14/3 , and G(3) = (7c - 86)/(5 x 8 !•).<_ 0 ,• c <_ 86/7 . Numerical evidence indicates that G(s) has no zeros in the interval (2,3). 14 The distribution of the zeros for a series of values of c in [2, — ] is shown in Table II. 34. TABLE II DISTRIBUTION OF ZEROS OF G(s) FOR 2<s<3 e The i n t e r v a l i n which a l l the r e a l zeros- of G(s) l i e 2 (1, 1.5) (3, 3.5) (12.5, 13). 2.5 (1,_1.5)_- (3, 3.5) (10, 10.5) 3 (1, 1.5) (3, 3.5) ( 9, 9.5) 3.5 (1, 1.5) (3, 3.5) ( 8.5, 9) 4 (1, 1-5) (3, 3.5) ( 8, 8.5) 4.5 (1, 1.5) (3, 3.5) ( 7.5, 8) 4.67 (1, 1.5) (3, 3.5) ( 7, 7.5) The function -10 G(s) , 2<s<3, i s tabulated below i n Table m f o r the same sequence of values of c. This table also shows that the function G(s) i s of constant sign i n [2,3] . TABLE I I I TABULATION OF THE FUNCTION -10 3 G(s) C 2 2.2 2.4 2.6 2.8 3.0 2 1.270 1.263 1.191 .987 .674 .357 2.5 1.032 1.051 1.038 .898 .631 .339 3 .793 .838 • .885 .808 .588 .322 3.5 .556 .626 .733 .719 .545 . . 3 0 5 -4 .317 .413 .580 .630 .502 .287 4.5 .079 .201 " .428 .541 .459 .270 4.67 0 .130 .377 .511 .444 .264 Thus i t follows that G(s) <_ 0 f o r 2<_s<3 . Combining the r e s u l t s proved above we get, G(s) < 0 for s e[0,4] , provided 2 £ c £ 14/3 , so that the generalized mean value theorem holds for these values of c. In a s i m i l a r way, i t can be proved that the operator associated with the A.E. (3.9) i s also d e f i n i t e and s a t i s f i e s G(s) <_ 0 for c>2 . Now we make use of the A.E. (3.8) at the i n t e r i o r points and the A.E.'s (3.9) at the boundaries and consequently the l i n e a r boundary value problem of class M i s replaced by A(c) y = b ...(3.11) where A(c) = J(c) + h 2 D(c) G and b = d- h 2 D(c) G . Here 36. J(c) 2 c - l 2-c -1 2-c. • 2c-2 2-c -1 2-c 2c-2 -1 2-c -1 2c-2 2-c 2-c 2 c - l ..(3.12) D(c) = 3 1 3 2 e3 3 2 3 3 l0 h h 4 1 50 d = c y a - h 2 3 Q f (a, y a)" Y a " h ' 30 f ( a ' y a } 0 0 y b - h 3 4 f (b, y b ) c y b - h 2 3 Q f (b, y b ) 2 the matrix G and vector s have been given before, Section 2.2. F i n a l l y , the e r r o r equation i s obtained i n the.form A(c) e = r(y) , . . .(3.13) where now, T(y) = 240 y u r - 31C-190 ,8 . (8) , 60,480 h y ( ? 2 } - 31c-190 ,8 (8) . 60,480 h / ( W c+1 240 i 6 (6) . h y ( y f o r any c s a t i s f y i n g 2<c £ 14/3 , and | |r(y) j | < h 2 K. . Here K = max. [ h 4 ( c+l) M 6 / 2 4 Q , h 6 |31c-190| Mg/60, 480] ...(3.14) 3.3 PROPERTIES OF MATRICES J(c) and A(c); We s h a l l prove that the matrices J(c) and A(c) are i r r e d u c i b l e and monotone. The geometrical concept of i r r e d u c i b i l i t y by means of graph theory i s quite u s e f u l , see Varga [25], Chapter I. We w i l l state here the main theorem which w i l l f i n d a p p l i c a t i o n • i n our a n a l y s i s . Let A= ( a ^ j ) be any (n*n) complex matrix. Consider any n d i s -t i n c t points P^, P2>..., P n i n the plane, which we s h a l l c a l l nodes. For every non-zero entry a.. of the matrix, we connect node P. to P. by means J . i j 1 2 of a path P^Pj directed from P^ to P_. , as shown i n the adjoining f i g u r e . When i = j and a _ ^ 0, we get what i s c a l l e d a loop. A directed graph i s strongly connected i f , for any ordered p a i r of nodes and P , there e x i s t s a path >P m n m =i r-1 r o P ± (i=j) P. P , P P connecting P. to P. i 3 Theorem 3.1 An (nxn) complex matrix. A = ( a_jj) 1 S i r r e d u c i b l e i f f i t s directed graph i s strongly connected. For proof see Varga [25], Theorem 1.6. From the above theorem we immediately deduce the following c o r o l l a r y . Corollary 3.2 I f A = (a..) i s a f i v e band matrix such that a± $ 0 ' ( i = 2,3,...,n) a i i+1 ^ 0 ( i = 1,2,...,n-l) . , then A i s i r r e d u c i b l e . Proof: In the proof of the theorem, we only require the directed subgraph G(A) of the matrix A = ^ a i j ^ namely G(A) : 1 "2 "3 ' 1 "n-2 n-1 n C l e a r l y the subgraph G(A) of the matrix A i s strongly connected, hence G(A), i t s directed graph, i s strongly connected as w e l l . Therefore X i s an i r r e d u c i b l e matrix. As a consequence., of t h i s c o r o l l a r y we deduce that the matrix 39. A(c) defined by (3.11) i s i r r e d u c i b l e i f "(i) (2-c) + h 2 3 g. i 0 ( i = 2, N-1) 2 1 ...(3.15) ( i i ) (2-c) + h & g ± ^ 0 (j = l , 3 ) , ( i = j,j+l,...,N+j-3) . These conditions w i l l be s a t i s f i e d f o r c>2 whenever i - 2 -h 6 CT < c-2 2 °M y ...(3.16) h 6j g M < c-2, (j = 1,3) . This i s the case when.either i ^ " = 3^ ' = 3^ = 0 or else f o r s u f f i c i e n t l y small values of h. . Having proved A(c) to be i r r e d u c i b l e , i t can be proved using Theorem 2.2 to be montone as w e l l f o r s u f f i c i e n t l y small values of h sat-i s f y i n g (3.16), provided the parameter c s a t i s f i e s the i n e q u a l i t y . . 2 < c < 18 ...(3.17) .In an analogous manner we can also prove that the matrix J(c) defined by (3.12) i s i r r e d u c i b l e and montone. The f a c t that the matrices A(c) and J(c) are both i r r e d u c i b l e and montone w i l l be c r u c i a l f o r the error estimates i n our l a t e r development. We further have A(c) >. J(c) i f c s a t i s f i e s the i n e q u a l i t y (3.17), hence from Theorem 2.3 i t follows that 0 < [ A ( c ) ] " 1 < [ J ( c ) ] " 1 ...(3.18) provided h and c s a t i s f y (3.16) and (3.17) r e s p e c t i v e l y . Now we proceed to determine | |.J(c) ^| |. The matrix J(c) can be factored as a product of two t r i d i a g o n a l matrices J and P^(c) and hence . j | J ( c ) - 1 | | = | |[J P N ( c ) ] _ 1 | | ...(3.19) ' . i l l P N 1 ( C ) ' N • l l J _ 1 H < I I P'^Cc) j I . (N+l) 2/8 40. where J i s defined by (2.4). The matrix P^(c) was introduced i n Section 2.4 and shown to have an inverse such that II?N 1 ( C )II <-H N(c)/(c-2) . We f i n a l l y obtain | | J - 1 ( c ) j | ^ H N ( c ) (b-a) 2/[8(c-2)h 2] . ...(3.20) Now going back to the error equation (3.13), we get | | e| | <. | | A _ 1 ( c ) | | | |r(y)| | ' ...(3.21) < H N(c) (b- a ) 2 K/8(c-2) , where K i s introduced•in (3.14). We summarize these r e s u l t s i n the follow-ing theorem. Theorem 3.3. 8 Let y(x) e C be the true s o l u t i o n of the l i n e a r boundary value problem of class M and l e t y^ be i t s d i s c r e t e approximation with an error e^ defined by the equation (3.13). Let A(c) be a monotone matrix such. that A(c) >_ J(c) , and l e t | JT(y) j J <_ h 2 K, then | |e| j i s given by (3.21) . 8 The numerical method i n which we use A.E. (3.8) with T.E. = 0(h ) at the i n t e r i o r points and A.E.'s (3.9) with T.E. 0(h ) near the boundaries w i l l henceforth be referred to as M-2(c). The error estimate w i l l be based on (3.21). 3.4 AN 0(h 6) FINITE-DIFFERENCE ANALOGUE Since the approximations near the boundaries i n the l a s t s ection were not as accurate as at the i n t e r i o r points, the o v e r a l l accuracy was 0(h 4) rather than O(h^) as desired, we therefore introduce more accurate A.E.'s instead of (3.9) near the boundaries.. We explain b r i e f l y how A.E. (3.9) i s modified to obtain new A.E.'s. The T.E. associated with (3.9) i s 41. 6 8 0 ( h ) and we want to develop new A.E. whose T.E. w i l l be 0(h ). This can be done by considering a A.E. of order k = 3 connecting y^, y^, y^, y^ and then Y N _ 2 > V N_I» Y N A N D Y N+l r e s P e c t i v e l y through the usual analysis i n the form 9 3 * (i) -c y Q + (2c-l) y x + (2-c) 7 2 ~ Y 3 + ^ I 3. y^ + ^ + 7 ^ ' ) = 0 \ 3 ( i i ) - y N _ 2 + (2-c) y N _ 1 + (2c-l) y N - c y ^ + h J y ^ _ 2 + . i=0 + h 3 (y2y^+Y3V^+1) = 0 , ...(3.22) where 3 Q = (58c - 57)/360 B = (91c + 16)/120 3 2 = (10c + 115)/120 3 3 = (-c + 24)/360 Y Q=-Y 3 = (2c - 3)/60 Y 1=-Y 2 = (c - 4)/20 and the T.E. associated with (3.22) i s given by ~71c-349 8 (8) , 349 302,400 y ' ° f 71 T.E. = -1,051 .9 (9) 349 , ... 4,294,080 h y , c = — f o r (3.22i) ...(3.23) 1,051 , 9 • (9) 349 , ..... _4,294,080 h y ' C=-Jl f° r (3' 22ll) Using the D.E. v i z . y" = g(x) y + s(x) , we get y'" =g*(x) y + g(x) y' + s' (x) . ...(3.24) We now substitute f o r y^' , y^' , y^' and . y ^ . i n the A.E.'s- (3.22) and thus these r e s u l t i n g equations contain "terms inv o l v i n g ,.y'^,, y^ and y ^ + ^ • In order to obtain expression for the f i r s t d e r i v a t i v e of y(x) at x x^, x^ and w e P r o c e e d as follows: We f i r s t consider A.E.'s of the form x 0' (i) I a y = h + h 2 I b y " + h 3 c y '" , 1=0 1 1 u i=0 3 3 ( i i ) I a y = h y ' + h 2 f d y" + h 3 c y'" i=0 i=0 • (3 The reason f o r including terms l i k e y V ' i n the above A.E.'s i s that the r e s u l t i n g T.E. associated with them can be made at l e a s t 0(h'') . The c o e f f i c i e n t s b's, c's, and d's are given as-follows: 'aQ = X + 2y - 1 a = -2X - 3y • + 1 where X and u are a r b i t r a r y r e a l constants, ients a's s a t i s f y the system a Q + a± + a 2 + a 3 = 0 Please note that the c o e f f i b„ = b, = a l + 2 a 2 + 3 a 3 = 1 (82X + 186 y + 271)/720 (94A + 192 y + 17)/120 (26.X + 258 y - 5)/240 (-2X + 24 y + D/360 43. d Q = ( 8A + 14 y - 7)/120 d1 = (194X + 402 y - 97)/240 d 2 = ( 16X + 138 y - 5)/120 d_ = ( -2A + 14 y + l)/240 c Q = ( 2A + 6 y + 7)/120 ^ = ( -6A - 18 y + 17)/120 and the T.E. associated with (3.25) is given by T.E. = 42(A-2y) - 53 ,7 (7) , - , Q O / . N \ n nn h Y > for(3.24i) 50,400 42(2A+y) - 53 , 7 (7) . 0 / . 50,400 h y , for(3.24n) Again substituting y^1 and y!^ 1 (using (3.24)) in A.E.'s (3.25) and solving for y^ and y^ , we get respectively 3 3 (i) y^ = h" 1 (l+h^go)- 1 [ I a y - h 2 I b yj - h 3c Q (g £y +s£) ] , i=0 i=0 1 9 •] 3 _ 3 _ ...(3.26) (ii) y» = h - 1 ( l + h ^ c g ) - 1 [ I a y - £ d y£ - h Jc ( g ^ + s ^ ] . i=0 i=0 In order to obtain expressions for y^ and y +^-^ > we consider A . E.'s (i) I * 3 - i yN-2+i = " h yN + h ' . 1 d 3 - i yN-2+i " h ' C l yN > i=0 i=0 . . .(3.27) 2 (ii) 7 a 0 . y.T 0 , . = - h y ' + h V b „ . y " . - h . rt 3 -1 •'N-2+i N+1 . L~ 3 -1 N-2+i i=0 1=0 c v'" 0 yN+l These A.E.'s (3.27) involving y^' and y^ .-^ are deduced from the A . E.'s (3.25) on replacing x by -x and then shifting the origin from x_^ to x T^_2 Again substituting y^1 aid y ^ (using (3.24)) in A - E.'s (3.27) and then solving for y^ and y^ +^ » w e 8 e t respectively 44. ( ± ) yN = h _ 1 I' I * 3 - i y N - 2 + i + h ' . t d 3 - i y N - 2 + i 1=0 i=0 3 h c f o ' v +sMl 1 N N N » (3.28) i l } y N + l = h _ 1 <1+h C 0 % + 1 ) _ 1 f-.I a 3 - i yN-2 + i + h ' X b 3 - i yN-2 + i i = 0 — - — i =o — - 2 + i ~ h 3 C0 ( gN+l yN+l + SN+1 } 3 • Thus the A.E.'s (3.221) and (3.22ii) are used at the two boundaries res-p e c t i v e l y , a f t e r y^ 1 , y^' , y^' and y^_ have been replaced using (3.24) and then the corresponding f i r s t d e r i v a t i v e terms v i z . y^ , y | , y^ and y ^ + 1 have been substituted using (3.26) and (3.28) r e s p e c t i v e l y . The optimum values of X and y can be obtained by solving the system "42 ( X - 2y) = 53 . ...(3.29) _42 (2X + y ) = 53 '', because the T.E. associated with the A.E.'s (3.25) and (3.27) w i l l be l e a s t i n that case. Solving (3.29), we get X = 53/70 , y = -53/210 , and with these choices of X and y the T.E. associated with the A.E.'s (3.25) and (3.27) becomes 0(h 8) . I f we now make use of these new modified A.E.'s at the boundaries which have been obtained by o v e r d i f f e r e n t i a t i o n and A.E. (3.8) at the i n t e r i o r points, then the l i n e a r boundary value problem of cla s s M i s replaced by A(c) y = b , ...(3.30) where A(c) d i f f e r s from A(c) given by (3.11) i n the f i r s t and the l a s t row only. The elements of A(c) that are d i f f e r e n t from the corresponding 3 elements of A(c) are given below, neglecting terms 0 ( h ) , a l l = 2 c - l i h2 3 1 § 1 + a l ( T0,0 + T ) a12 = 2-c + h 2 6 2 g 2 + a 2 ( T0,0 + T ) i , r a l 3 -1 + h 2 6 3 g 3 + a3 ( T0,0 + T ) 1,3/ *N,N- -2 -1 + h 2 5 3 % - -2 + a3 ( T1,N + T0,N+1 *N,N.--1 = 2-c + h 2 52§N--1.. + a2 ( T1,N + T O.N+l aN,N = 2 c - l + h 2 g l % + a l ( T1,N + T O.N+l T. . = ,2 (1 + h 2 c. X r 1 • ) ) ) , We would l i k e to e s t a b l i s h - c r i t e r i a f o r A(c) to be i r r e d u c i b l e and mono-tone such that A(c) > J(c) . We f i n d that the matrix A(c) defined by (3.30) i s i r r e d u c i b l e i f a. _ 4 0 , a.T X T . i> 0 and (2-c) + h 2 B.g„ + 0 1Z IN,JM—i j M. (j = 1,3) . These conditions w i l l be s a t i s f i e d for 0c>2 whenever 0 < h 2 3 2 g 2 + a 2 ( T 0 > 0 + T l 5 l ) < c-2 0 ^ ^ V N - 1 + a2 ( T1,N + T0 , N+1 ) < C " 2 * • ' ( 3 ' 3 ] h 2 3 J § M < c-2 (j=l,3) Since & > 0 for a l l values of c s a t i s f y i n g (3.17), the l a s t condition i s s a t i s f i e d for s u f f i c i e n t l y small values of h. We now proceed to examine the f i r s t two conditions of (3.31). We have seen above that X = 53/70 and U = -53/210. With t h i s choice of X and u we f i n d a x =17/70, a 2 = 53/70, a 3 = -53/210; c Q = 7/120, c± = 17/120 . The c o e f f i c i e n t s 3's and Yg are p o s i t i v e for any value of c s a t i s f y i n g (3.17). However y = (c-4)/20 < 0 , for c s a t i s f y i n g 2 < c < 4 . We now consider 46. 2 2 h 2 6 2 g 2 + a 2 ( T 0 j 0 + T 1 } 1 ) = h 2 3 2 g 2 + a 2 h 2 ^ ( 1 - ^ ^ + ^ ( 1 + % - h 23 2( g l+hg|+...)+a 2h 2[Y 0( g l-hg;+...)+y l g l] ^ *2 rlOc+115 , 53 ,2c-3 c-4,, - h g x [ — + — < — + — ) ] ^ 123c+646 .2 — h 840 °1 ' S i m i l a r l y • ' 2 2 2 2 2 171*1 J. 7h h * 2 % - l + a 2 ( T l , N + T 0 , N + l ) = h 3 2 g N - l + a 2 h [ Yl§N ( 1' f12b- % ) _ ^ O g N + l ( 1 4 1 2 o W ^ 123c + 646 ,2 — n o-840 °N ' neglecting a l l terms 0 ( h ) . Now i t i s evident that the f i r s t two conditions (3.31) w i l l be s a t i s f i e d i f 123c + 646 ,2 0 < h gM < c " 2 > and t h i s can be forced to.be s a t i s f i e d for. s u f f i c i e n t l y small values of h. Now we want the matrix A(c) > J(c) such that A ^(c) <_ J "'"(c) . This w i l l be the case i f the following s i x conditions are s a t i s f i e d , namely h h g i + a i \ o + V - ° ' h ' 3 4 - i gN-3 +i + V i ( T1,N + T0,N+1> ^ 0 > ..(3.32) ( i = 1,2,3) . Again through the usual analysis employed above i n e s t a b l i s h i n g the f i r s t two conditions (3.31) we prove that the conditions (3.32) w i l l be s a t i s f i e d i f " 654c + 61 2 ( 1 ) " W h g m > ° ( i i ) 123c + 646 ,2 •• . — ^ 4 0 h gm > ° ..(3.33) -60c + 327 . 2 47. at l e a s t for the values of c s a t i s f y i n g 2 < c <_ 14/3. But the above three conditions are obviously s a t i s f i e d for those values of c. Thus the conditions (3.32) are s a t i s f i e d at le a s t for c s a t i s f y i n g 2 < c <_ 14/3 , and for s u f f i c i e n t l y small values of h. Having proved that A(c) i s i r r e d u c i b l e and monotone such that 1 _1 A(c) > J(c) and A (c) < J ( c ) , we obtain the error equation i n the f o l -lowing form. : A(c) e = r (y) , W h e r e r 71C -349 (8) 3 0 2 , 4 0 0 y u r 31C -190 ( 8 ) ( f , 6 0 , 4 8 0 y r ( y ) = h e 31c-190 (8) , • . y -,) ..(3.34) 60,480 71c-349 (8) . . y v ^\T / -302,400 N-1"2 = for any value of c s a t i s f y i n g 2 < c _< 14/3 , and | |r(y) | | <_ h K Here K = max. h M o e a s i l y show I71c-349 71c-349 302,400 |31c-190. 60,480 , 349 190 One can K = 3 0 2 , 4 0 0 3 1 c - 1 9 0 l 6 0 , 4 8 0 , 6 M _ 1299 h M , 5.75 = 266 -601 _ . c < c < — — = 7.15 84 h M , otherwise. Thus for c s a t i s f y i n g 2 < c _< 14/3 , K turns out to be equal to 190-31c ,6, 6 0 h Mg . Now s u b s t i t u t i n g K instead of K i n (3.21), we get | |e|.j. <• h 6 ( 1 9 Q - 3 1 ^ ± ^ l 8 „ % ^ L . • ., 2- < c-.< 14/3 , " "• 483,840 (c-2) ~ » (3.35) 48. or. equivalently 0(h 6) . The above error bound i s analogous to the one obtained e a r l i e r (Theorem 3.3). The method described above w i l l henceforth be referred to as M-3(c). The r e s u l t i n g error i s based on (3.35). J 4 Recently an 0(h ) f i n i t e d ifference analogue was formulated by Bramble and Hubbard [3]. The c o e f f i c i e n t matrix was based on a f i v e point A.E. ( ± ) h \ = 12 (" yn-2 + 1 5 y n - l " 3 0 y n + 1 6 y n + l " yn+2> > at the i n t e r i o r points and by a three point A.E. 2 ( i i ) h y" = y , - 2 y + y , , ,„ „,N •'n n n+1 ...(3.36) at the boundary points of the range of i n t e g r a t i o n . In f a c t the problem considered by Bramble and Hubbard was the one given by (2.19). The main r e s u l t of t h e i r paper was that the d i s c r e t i z a t i o n error s a t i s f i e s the i n e q u a l i t y A - M M , • I h i | 1 h [ +'T§O ] • ...(3.37) In our case [ |e| | based on f i v e point A.E. (3.8)- i s -given by (3.21), which turns out to be a better bound than the one given by (3.37). We postpone comparison of the method proposed by Bramble and Hubbard (henceforth r e f e r r e d to as M-BH) with other methods t i l l the l a s t chapter. CHAPTER IV A CLASS OF HIGH ACCURACY DIFFERENCE FORMULAS OF LOW ORDER 4.1 In the present chapter we consider A.E.'s of order k = 2 which are of the form E a. y ,. = h 2E 8. y" . + h 3E y. y"* + . . .+ hfc E 9. y ( ^ , ...(4.1) • x n+x x Jn+x x Jn+x x Jn+x ' where t is a positive integer >_ 3 . In this chapter a l l summations extend from i = 0 to i = 2 . Since y" = g(x) y + s(x), in general y ^ ( x ) , n >_ 3, cannot be expressed in terms of y(x) and higher derivatives of g(x) and s(x) alone. The derivatives of y(x) higher than 2 w i l l contain y'(x). Therefore we need expressions for y n + ^ > i = 0,1,2. These expressions can be obtained by considering auxilliary A.E.'s of the form E a. y , . = h y' + h 2E b. y" +. . .+ h tE <j>. y ^ , ...(4.2) x Jn+x •'n+m x n^+x Yx yn+x (m = 0,1,2) . In general the coefficients b's, c's,... etc. depend on m. We obtain by successive differentiation of .the D.E. viz. .y" = g(x) y + s(x) (4) the expressions for higher derivatives y"' , y , ... etc. Assume y ( t ) ( x ) = 0 t(x) y + / ( x ) ,y' + / ( x ) , t >. 3 ...(4.3) t t t where the functions 9 (x) , <j> (x) and ip (x) depend on g(x) , s(x) and on their higher derivations. On substituting the expressions for higher derivatives given by (4.3) in A.E. (4.2), we get hy^ (1+h2 c Q ( ^ +...).+ h 3 y ^ + I •(c1< +^1+...-) "+ h 3y^ + 2(c 2 <^ + 2+...0 + n(m)=0 ...(4.4a) h 3yn ( c0*n +"- ) + h yn+l ( 1 + ^ • •) + h 3 7 n + 2 (c 2 <j>3+2+. ..) + n(m)=0 ...(4.4b) 50. h \ ( C 0 < +--' ) + h 3 y ; + l ( c l * n + l + ' " ) + hy; +2 ( 1 + h 2 c2*n+2 +'--> + ^ ) = ° • " ( 4 -2 3 where n(m) = - £ a. y ,. + h £ [(b. g ,. + h c . & ,. +...) y ,. l n+i x °n+i x n+x "'n+x + (b . s , . + h c. lb3-, . + . . .) ] . x n+x x rn+x We solve the above system (4.4) for the unknowns y^ , y^ +^ and y^ +2 • However we can simplify the above system (4.4) to a very great extent i f we consider the auxilliary A.E. (4.2) in the form Z a. y ,. = hy', + h 2£ b. y",.+h3c y"' + x n+x n+m x n+x m n+m h 4d y (t? +...+ht $ y (*2 , (m = 0,1,2) ...(4.5) m n+m m n+m Again on substituting in (4.5) the expressions for higher derivatives of y(x) as given by (4.3) and solving i t for y ^ ^ > we get h y* = -(l+h 2d <j>3 + h 3d <f>4 +. . . ) _ 1 n , , ...(4.6) -'n+m myn+m mYn+m n+m ' v J 2 where n , = -£ a. y ,. + h E (b. g ,. + s ,.) n+m x n+i i ton+i n+i + h 3 (c 0 3 + hd 64, +...) y , m n+m m n+m J n+m + h 3 (c i>3 + hd 4;4 +...) , (m = 0,1,2) . m n+m m n+m ' (4) We now obtain expressions for higher derivatives y"!., y ,., ... etc. in • , n+i Jn+i the form, (using (4.3) and (4.6)), = y ~ h _ 1(l+h 2c 4>3 +h3d ci 4 +...)"1 n _,+ ^ . ,...(4.7) n+m n+m n+m n+m m n+m m n+m n+m n+m (m = 0,1,2) , t _> 3 . We • substitute the expressions for higher derivatives y ^ ^ ( x ) , t >_ 3 as given by (4.7) into the A.E. (4.1) and get.the required A.E. to be used in the numerical.solution of the linear boundary value problem of class M. In what follows we w i l l consider the A.E. (4.1) for t = 3 and 4 in detail and the auxilliary A.E.'s (4.5) for t = 2 and 3 respectively. 51. It i s proved that the resulting error based on these methods i s 0(h q), q = 5,6,7 and 8. 4.2 THE CASE FOR t = 3 It can be shown (as in (2.3)) that the A.E. (4.1) for any pos-i t i v e integer m is necessarily of the form y — 2y . + y - = h 2 E 8. y " . +...+ h C Z 9. y ( ^ , ...(4.8) •'n •'n+l 'n+2 x •'n+i x yn+i ' A special case of (4.8) for t = 3 is y - 2y - + y _ = h 2 E 3. y",. + h 3 E 6. y"' , ...(4.9) Jn Jn+1 Jn+2 x J n+x x Jn+x ' where B Q = 6 2 = 2 / 1 5 '> 6 1 = 11/15 Y Q =~Y2 = 1/40 T l - o and T.E. of (4.9) i s given by T ' E - _ 302,400 h y Since y^ = 0 , therefore we only require auxilliary A.E.'s of the form (4.5) for m = 0 and 2. We f i r s t consider (4.5) for t = 2 , m = 0 namely E a. y , . = h y' + h 2 E b. y" . ...(4.10) l •'n+i Jn x Jn+x Using Taylor's formula for expanding y n + ^ and y^ +^ around x , one can show that a's used i n the A.E. (4.10) satisfy the following system: a 0 + a± + a 2 = 0 ax +2a2 = 1 (4.11) The general solution of the system (4.11) i s given by a Q = x-1 , a1 = -2\ + 1, a 2 = x, ...(4.12) where x i s an arbitrary real number. Using the above values of a^, a^ and a 2 we find through the usual analysis b Q = (2A+7)/24, b ="(10X+3)/12 and b 2 = (2A-l)/24, and T.E. associated with (4.10) is given by T.E.'-Tf h 5 y< 5 ) . 45 We might remark that the A.E.'s (4.9) and (4.10) are both such that the operators associated with them are definite and the generalized mean value theorem holds. The A.E. (4.5) for t = 2 and m = 2 is symmetric to the A.E. (4.10) and can be obtained by replacing x by -x or by flipping the inter-val [x , x + 2h] in the form n n . , Z a 2 - i yn+i " ~ h K+2 + h ' E b 2 - i Ci ' • • • ^13> The T.E. of (4.13) i s given by T.E. - - ^ h 5 y « > . , On substituting the value of y^ +^ using the D.E. i n (4.10) and (4.13) and solving for y^ , y ^ + 2 respectively, we get '(i) h y' = Za. y , . - h 2 Z b. (g . . y , . + s , .) ^n l •'n+x x '"'n+x yn+x n+x n / . 2 • • •(4.14) ( i i ) h y' = - Z a , . y , . + h Z b 0 . (g .. y + s , .) Jn+2 2-x n+x 2-x ton+x n+x n+x Similarly on substituting y ^ + i . y n + i (us:*-n§ D-E.) in the A.E. (4.9) and then substituting for h y^ , h y ^ + 2 using (4.14), we get the required . A.E. which i s used to approximate the boundary value problem in the form (-1 + A n) y n + (2 + B n + 1 ) y n + 1 + (-1 + C ^ ) y n + 2 + = 0 , ... (4.15) (n = 0,1,2,..., N-1) , 2 where A ± = h (&Qg± + Y 0 a 0 § i ~ Y 2 a 2 § i + 2 ) ( ± = O'1.---.1 -^!) B. = h 2 ( 3 l g i + Y 0 a i s i - 1 - Y 2 a l g i + 1 ) ( i = 1'2>"->N> 2 C ± = h ( 3 2 g ± + Y 0 a 2gi_2 " Y 2 a 0 g i ) ( ± = 2» 3>---> N + 1) 3 plus terms 0 ( h ) . 53. Here 6 Q = 2/15, 3 1 = 11/15, 3 £ = 2/15, 3 Q = - Y 2 = 1/40, and a. Thus X-1., a^ = -2X +1, a 2 A. = h (6 0g. + y0a0g± - Y 2a 2g. + 2) = h (6 0g. + Y n a n g , " Y,a 9(g, + 2h g! +...)) '0 0°i 2 2 V 6 i - h 2g.(S 0 + Y 0a Q - Y 2a 2) - h g i ( 1 5 + 40 ^ 6X+13 ,2 0 , n „ - 1 2 Q . h g . > 0 for X> -13/6 Similarly and ^ 47-^ 6X ,2 • _ , , -, ir B. = — — — h g. > 0 for X< 47/6 x 60 x c i = h 2 § i > 0 f o r A > " 1 3 / 6 ' • neglecting the. terms 0(h ) . Thus for sufficiently small values of h, A.'sv B.'s and C.'s are positive x x x v and also A^'s, C^'s are each less than unity provided (4.16) Now the associated matrix with the system of linear equations (4.15) has a l l the nice properties given by (2.8) provided the conditions (4.16) are satisfied. We can now obtain an error bound of the form (2.9). In fact we get | e i ' -h • M5 % ( x i _ a ^ ( b _ x i ) 1,800 e| | i h 5 M 5 g M (b-a)2/7,200 , (4.17) or equivalently ||ej | = 0(h ) . In this case we w i l l not be able to obtain the more accurate error bound of the form (2.15) for the simple reason that we w i l l not be able to associate with the matrix A, the matrix A (defined in Section m 2.4) with the desired properties v i z . A>A 'and A is symmetric and can m m be factored as a product of ? N ( v ) and Q . 4.3 MORE ACCURATE EXPRESSIONS FOR y', , (m = 0,1,2) —• : n+m — — —*~ We have seen that although the A.E. (4.5) is such that i t s 8 * 5 T.E. = 0(h ) , yet the error resulting from i t is 0(h ) . The reason is that the T.E. of the auxilliary A.E.'s is O(h^), and when we substitute y^ and y ^ + 2 using (4.14) in (4.9), the resulting A.E. turns out to have a T.E. 0(h 7) . Thus the accuracy, of the numerical procedure discussed above mainly depends on the accuracy of the auxilliary A.E.'s. Therefore we would like to obtain more accurate auxilliary A.E.'s. Hence we now consider (4.5) with t = 3 and m = 0 namely Za. y , . = h y' + h 2 z b. y" . + h 3 c n y"' , ...(4.18) x ^n+x -'n x yn+x O n ' where a's are given by (4.12) and b's are given as follows: b Q = (10X + 47)/120 b1 = (50A + 7)/60 b 2 = (10A - D/120 c Q - 1/15 ^ T ' E - ^ h 6 y ( 6 ) * 1 / 2 ' T.E. = 1,575 h7 y ^ , A = 1/2 We similarly develop a A.E. containing h y n +-^ in the form Za. y , . = h y' . + h 2 Z d. y", . + h 3 c, y"' , ... (4.19) x Jn+x Jn+1 x Jn+x 1 Jn+1 where d Q = (5X - 4)/60 d = 5(2A - 1)/12 d 2 = (5A - l)/60 c = 7/60 and T.E. = Wo * 6 y ( 6 ) .x*m 11 h7 (7) _ . 50,400 h y ' A ~ 1 / 2 However A.E. (4.5) for t = 3 and m = 2 can be obtained using A.E. (4.18). We replace x by -x, y 1 by -y', y'" by -y'" in (4.18), rearrange the suf-fixes and get Ea. y ,. = -h y' + h 2 £ b_ . y" - h 3 cn y"' , ...(4.20) i n+i yn+2 2-i yn+i 0 n+2 where the T.E. of (4.20) i s the same as that of (4.18) when X + 1/2 , and is of oppostie sign for X = 1/2. Now following the analysis used i n the previous section, we can prove that the resulting error in the boundary value problem is 0(h ), provided the conditions (4.16) are satisfied. This result w i l l also be a special case of the Theorem 2.4 for p = 6 and M = 29 Mg/302,400. The inequality (4.21) gives the error bound based on this method. , 29 h Mg (b-a) l U M i . 2,419,200 ...(4.21) 6 or e = 0(h ). The more accurate error bound of the form (2.15) does not apply here for the same reasons mentioned in the previous section. The numerical method in which we use A.E. (4.9) and auxilliary A.E.'s (4.18) and (4.20) w i l l henceforth be referred to as M-4. The error estimate i s given by the inequality (4.21). 4.4 THE CASE FOR t = 4 We now consider (4.8) for t = 4 namely y - 2y i + y , 0 = h 2 £ 8. y" . + h 3 E y • y " l - + h 4 E <S. y ( 4 ^ , ... (4.22) Ju yn+l 7n+2 I 'n+i 'i ' n + i l 'n+i ' where 3 o = 3 2 = 912Q " V= 8,256Q Y o =-Y2 = 177Q Y l = 0 = 6 2 = HQ 61 = 260Q Q = 1/10 ,080 . and The T.E. of (4.22) is given by T p -17 h12 (12) 5,588,352,000 y If we use D.E. to substitute expressions for y", y'" and y ^ i n (4.22) and then substitute expressions for h y', (m = 0,1,2) obtained n+m from auxilliary A.E.'s (4.18), (4.19) and (4.20), we get the required 9 A.E. The T.E. of this A.E. is 0(h ), again limited by the accuracy of the expressions approximating Y^^C 1 1 1 =0,1,2). The matrix associated with the system of linear equations w i l l be monotone and >J provided X and h are chosen so as to satisfy the following conditions. "(i) -245/118 < X< 1435/118 n , 118X+245 .2 ^ • ...(4.22) j x l ) ° < 3,360 h % " 1 We omit the details and state the error bound obtained for a linear bound-ary value problem of class M below in (4.23). 177 h 7 M? g M (b-a) 2 H e l l ^ : 63,504,000 ...(4.23) or e = 0(h 7). Please note that the error estimate given by (4.23) is a special case of the theorem 2.4 for p = 7 and M = 177 M., gM/7,938,000. The more accurate error bound of the form (2.15) does not apply here as well. 57. where 4.5 AN 0(h 8) FINITE DIFFERENCE ANALOGUE Another special case of A.E. (4.8) w i l l be V 2 y n + l + yn +2 - h ' E 6 i K+± + y n + i • (4-24) 8 0 = 6 2 = 660Q &1 = 13,800Q Y Q - Y 2 - -13Q Y l = 626Q 15.120Q = 1 and the T.E. associated with the above A.E. i s given by 10 (10) T.E. -59 h y' 76,204,800 The operator associated with the A.E. (4.24) can be easily shown to be definite. Again using the D.E. and the auxilliary A.E.'s (4.18), (4.19) and (4.20), we can obtain expressions for y ^ ^ (m = 0,1,2). Replacing y" and y ^ i n (4.24) by using D.E. and then substituting explic i t expressions for > t n e required A.E. is obtained with a T.E. 0(h"*"<">) . Here the D.E. i s approximated by a system of linear equations of the form (4.15) where now A. = h 2 60g.•+ 0(h 3) , B. = h 2 e i g ± + 0(h 3) , . and C. = h 2 B„g. + 0(h 3) . x 2 x Since 6. > 0 , therefore A., B., C. > 0 and i f h be chosen so that i ' x' x' x 0 < h 2 B i g M < 1 ( i = 0,2), then the associated matrix A w i l l be monotone and > J . As before, for sufficiently small values of h we state the error bound obtained i n the form h 8 (590 M 1 Q + 15,436 M? g^) (b-a) 2 ''e.''- 6,096,384,000 ' ...(4.25) The method given above w i l l henceforth be referred to as M-5, where the error estimate is given by (4.25). We have seen that the accuracy of the numerical method mainly depends on the accuracy of the auxilliary A.E.'s of the form (4.5). We have considered (4.5) for t = 2, and 3 in the previous discussion. By choosing A.E.'s (4.5) for t > 3, and combining them with a suitable A.E. of the form (4.8), we can produce numerical methods in which e = 0(h ), r >_ 9. The cases in which r = 5, 6, 7 and 8 have been discussed in this chapter in detail. 4.6 LINEAR BOUNDARY VALUE PROBLEMS OF CLASS M, WHERE g(x) = Constant We now consider a special class of two point boundary value problem of class M of the form (1.1), where f(x,y) = a y + s(x) . Here a is a real positive arbitrary constant. For this class of problems extremely accurate f i n i t e difference analogues can be developed by follow-ing the technique, discussed in this chapter, without using any auxilliary A.E.'s of the form (4.5). The reason, for this is the following: From v / o i \ ----(2v) v , r v-m (2m-2) . . y v = a y + l . a s (x) , m=l (2v+l) v , r v-m (2m-l), . and y = a y'.. + ... £ a- s (x) , m=l we notice that the higher derivatives of y of even order do not contain y', hence A.E. (4.8) can be modified in the form Y ~ 2y x l + y"V0" - h 2 z 3. y" . + h 4 Z y. y ( 4 ^ +...+ h 2 t Z 9. y ( 2 t ) ,... (4.26) •'n -^ n+l n^+2 pi> n+i ' i yn+i i •'n+i where t is a positive interger >1.' We w i l l consider two cases of (4.26) 59. in detail namely for t = 2 and t = 3. The A.E. (4.24) already considered is a special case of (4.26) for t = 2. . We w i l l not repeat i t here again. However, i f we approximate the two point boundary value problem of class M with constant coefficient by A.E. ( 4 . 2 4 ) , then the resulting error is 8 0 ( h ) . In fact the error bound obtained is as follows: 59 h 8 M Q (b-a) 2 " E " - 609,638,400 " . . . ( 4 . 2 7 ) The error bound for linear problem with constant coefficient obtained here is a special case of (2 .9) for p = 8 and M = 59 M Q/76,204,800. In this case however, we can obtain a more precise error bound using the . theory developed in Section 2 .4 . Without giving the computational details, we state a more precice bound of the form (2.15), below 8 5 9 h " M10 V y ) - 76,204,800 a (4.28) where y = (2+h2 a&'+ h 4 a 2 y 1 ) / ( l - h 2 a 8 Q - h 4 a 2 y ) > 2 for a l l sufficiently small values of h > 0. We now consider (4.26) for t = 3 viz. y - 2y + y 0 = h 2 E B. y" + h 4 E y. y ( 4 J + h 6 E 6. y ^ . , ...(4.29) 'n 'n+1 'n+2 x 'n+x 'x 'n+x x n+x where the coefficients B's, y's and <5's are obtained through the usual analysis and are given by eo - B 2 = 1,154,160Q B l = '36,943,200Q Y0 - Y 2 = -16,632Q Y l " 2,150,064Q 6o " 62 = 127Q • ^ 1 " 29,230Q and 39,251,520Q = 1 60. Also the T.E. associated with (4.28) is given by T E = ~45>469 h 1 4 J1^ 13.! x 272,580 y Then since the B's are positive, i t follows from (2.9) that 45,469 h 1 2 M (b-a) 2 l | e " ~ 131 x 2,180,640 ' ...(4.30) for a l l sufficiently small values of h. In this case as well we derive a more accurate error bound. It follows from (2.15) that 45,469 h 1 2 M 1 4 l y y ) H 8" - . 13! x 272,580 a ' ...(4.31) y = (2+h2 a $ 1 + h 4 a 2y^ + h 6 a 3 6 ; L)/(l-h 2 a B Q -^h 4 a 2y 0 - h 6 a 36 Q) . The numerical method based on A.E. (4.29) w i l l henceforth be referred to as M-6, where the error estimate i s based on (4.31) ... CHAPTER V LINEAR BOUNDARY PROBLEMS WHERE 8 £ ( ' X , Y ^ <0 . . 3 y 5.1 We know that the solution of the boundary value problem of class M namely y" = f(x,y), y(a) = y^, y(b) = y b is unique, provided the function f(x,y) , in addition to satisfying the usual conditions of the existence theorems i s such that ( x?y) >Q ^ n g 3y -8 f fx v) However i f — 3 y <®' t n e n t^ i e solution of the above boundary value problem i s not necessarily unique. For instance, the boundary value y" + T7 2y = 0 , y(0) = y(l) =0 has y(x) = C sin x as a solution for arbitrary values of C. Thus this problem has no unique solution. On the other hand, the solution of y" + T r 2y = x + T T 2 x 3/6, y(0) = 0, y(l/2) = 1/48 is unique. In fact y" + A 2y = g(x), y(a) = y(b) =0 has a unique solution whenever X =f mr/(b-a), n = ±1,±2,.... The purpose of this chapter i s to obtain error bounds for such problem where the solution is unique. , 5.2 Let us consider the problem in the form y" = -g(x) y + s(x) , g(x) > 0 in S, ...(5. y(a) = y a, y(b) = y b .' We assume that the above boundary value problem (5.1) has a unique solu-tion. We ,then replace the D.E. (5.1) by the A.E. (2.3) namely y n " 2 y n + l + yn +2 = ^ / 1 2 ) %K + S l y n + 1 + * 2 ^ 2 > 62. which has a truncation error determined through T.E. = C P+2 h P + 2 y ( p + 2 ) ( 0 j ? e [ X n > w ? where 8 Q + B 1 + 8 2 = 1, B Q = 8 2 , 8.. >_ 0 (j = 0,1,2) . We have already shown that the operator associated with the above A.E. is definite, see Section 2.1. The error equation is obtained in a usual manner- and can be put in the form Ae = T(y) ...(5.2) where now A = B, N-2 V l BN BN-1 \ so that A = 2 - h 6 g n I n B = -1 - h 2B ng , (n = 1,2,....N) n u n - d ||r(y)|| - |C p + 2| h P + 2 M p + 2 . We shall at f i r s t consider a special case of D.E. (5.1) in which g(x)= constant = g, say. Using the results of sthe section 2.4, we get . - l i N (1 + h BQg) sinG| sin(N+l)6| (5.3) where 2 cosO = y = (2 - h 2Qlg)/a + h 2B Qg) < 2 ~1~T Setting. sine = h Jg + hV (3Q - 3^/4) / (1 + h"BQg) , we f i n a l l y get |C | hP VL (b-a) l l ' l l , , P+2' P+2 /g+h2g2(B2-B2/4)|sin(N+l)6| ...(5.4) We state the above result in the form of a theorem as follows: Theorem 5.1 Let y(x) e C P + 2 be the true solution of the boundary value problem y" = -g y + s (x) , y(a) = y , y(b) = y, , where g is a real posi-cL D tive constant; also let y^ be i t s discrete approximation based on the A.E. (2.3), then I C o.o I h P M j.o (b"a> 11*11 < P+2 2±? : /g + h2 g2(02- B2/4)|sin( b 1^)G| ' . ~" '" or ||e|| = 0(h P) .where 2cos9 = (2 - h 2 3 1 g ) / ( l + h 2B Qg) < 2 Particular Cases of Theorem 5.1: Case (a) - A special form of A.E. (2.3) is of degree p = 2 with 8 Q = 3 2 = 0, B x = 1 and T.E. = ^ y ( 4 ) ( 5 ) , hence . . . h 2 M4(b-a) 12 /g-h^g^/4|sin(^#e) (5 2 where 2cos9 = 2-h g. Case (b) - Another and more accurate form of A.E. (2.3) is of degree p = with 6 0 = 3 2 = 1/12, 3 1 = 10/12 and T.E. = -h 6 y ^ ( £ ) 7 240, hence h 4 M,(b-a) o 240 /g-hzgz/6 | sin ( ^ 9 ) | where 2cos9 = (24 - 10 h 2g)/(12 + h 2g) . We know that the D.E. under consideration v i z . y" = -g y + s(x), y(a) = y a, y(b) = yfe does not have a unique solution whenever 0 < •/% (b-a) = n i r , n = l , 2 , 3 , . . . We also have 2cos9 = (2 - h 2 3 1 § ) / ( l + h 23 Qg) < 2 . ... (5 Expanding cos9 in powers of 0 and the expression on the right, hand side in powers of h, we get 2 4 2(1 - f j - + fr---> = ( 2 - h 2 hs)a - h 2 e 0g + h 4 B^g2-...) , = 2 - (23 Q + 3^ h 2g + 0(h 4) , = 2 - h 2g + 0(h 4) , or 0 ^ /g h . Consequently sin(N+l)G ^ sin [(N+l)/gh] ^ sin ((b-a)/g) ^ 0 i f 0<(b-a) /g = n 'ir , (n = 1,2,.. .) . Therefore the expression on the right hand side of•(5.4) which has sin(N+1)0 in i t s denominator could be very large ( i f not infinite) pro-vided the range of integration (b-a) is such that 0<(b-a)/g"2 n I T . (n = 1,2,3,...) , ...(5. and in that case mathematically the system of linear equations arising from the D.E. may have a unique solution but numerically the problem is ill-conditioned. We consider the boundary value problem . y" + E 2 y = sin x, (E 2 $ 0, 1 and E ^ where n is an integer). y(0) =1 , y(b) = sin(Eb) + cos (Eb) + sinb / (E 2 - 1), see (8.2). We choose those values of E and b whose proudct approximately equals nir, for some positive integer n. For example, for E = 2 and b = 11, Eb = 22 ^ 7tr. Note that 7TT is approximately equal to 21.99. The "max. absolute error" in the numerical solution with h = 1/8 is found to be 0.20 E - l and the theoretical estimate based on the inequality (5.6) is equal to 0.56 E - l . The ratio of "max. absolute error" to i t s estimate is 0'.35. On choosing E = 11, and b = 2, the "max. absolute error" and, the error estimate turn out to be 0.10 E+l and 0.28 E+l respectively, and their ratio being 0.36. On taking E = b = 5 (Eb = 25 2L 8TT) , the "max. absolute error" is equal to 0.68 E - l , while the error estimate based on (5.6) is 0.18. The ratio of "max. absolute error" to i t s e s t i -mate is 0.38. We tried several other combinations of E and b. The result are similar to those discussed above. The above experimental results confirm that the expression for error bound (5.4) is predominantly range dependent and whenever range (b-a) satisfies (5.7), the corresponding "max. absolute error" and i t s estimates both become large. 5.3 AN 0(h 8) ERROR BOUND We can obtain more accurate results of the form (5.4) for the boundary value problem (5.1) with g(x) = constant, by using A.E.'s of the form (4.26). We f i r s t consider (4.26) for t = 2 and show that under cer-g tain conditions the resulting error is 0(h ). Here the resulting system of linear equations w i l l be of the form (5.2) where 2 4 2 A n = 2 - h Z 3 1 § + IT YI8' , B n - -1 - h 2 3 Qg + h 4 Y ( ) g 2 , 59 M 1 Q h i U a n d l l r ( y ) " - 76,204,800 ' for the coefficients 3^'s and Y^'s > s e e A.E. (4.24). The constant u introduced in Section 2.4 w i l l now be given by u = (2-h233_g + h 4 Y ; Lg 2)/(l+h 2e 0g - h 4 Y ( )g 2) . The denominator of u is clearly greater than 1 because 3^ > 0 and y^ < 0. The numerator w i l l be'less than 2 provided 31 > h 2 y l g or on substituting the values of and y , we have h 2g < 6900/313 . Thus n < 2 i f h satisfies the above inequality. Finally using the theory developed in Section 2.4, we get MA"1!! •< N/(l+h 2 8 0g - h 4Y Qg 2) sine | sin (N+1) 9 | , where 2 cos6 = y. < 2 , and l l.l l., ^ N Mr(y)ll (l+h 2 8 0g- h 4Y Qg 2) sin9|sin(N+1)6 On substituting (l+h 2 8 0g- h 4y 0g 2) sine = h yg+h 2g 2(B 2-e 2/4)-h 2g 2(2Y 0+Y 1)-2h 4g 3(B 0Y 0 -8 1Y 1/4) • u6 4, 2 2 IIS + h g (YQ-YJ^/4) , = h /g + 0(h 2) , we get s 59 M h° ||e|| < , ..'.(5.8) 76,204,800|sin(N+1)6|/g + 0(h z) or | |e-| | = 0(h 8) . In a similar manner we can use A.E. (4.26) for t = 3 and show 12 that the resulting error in the corresponding case is 0(h ). We omit the details here for brevity. 67. 5.4 THE GENERAL CASE WHEN g(x) 4 constant In determining the error bounds (5.4) and (5.8) respectively, we have so far assumed that g(x) is constant. The general case when g(x) i s not a constant cannot be dealt with by the above method. We suggest that error bound should be obtained numerically. We consider the error equation (5.2) namely Ae = T(y) . Thus I M l 1 |C p + 2| h P + 2 M p + 2 MA" 1 ! ! . ...(5.9) In order to obtain the error bound given above, we w i l l compute A ^ and hence the norm | |A | | numerically. In actual practice, the error estimate obtained by the above technique is not very good, especially when the range of integration i s - ' ' " " " f a i r l y large. By the error estimate not being very good, we mean that the ratio of maximum absolute error to i t s estimate is very small. We tried the error estimate given by (5.9) in various situations. The results for the problem (5.10) are displayed i n Table IV on the next page. y" + x 2y = x + x 3 + 2/x3 , ...(5.10) y(l) = 2, y(b) = b + £ with y(x) = x + ^ , and <_ n ! . The Table IV shows that for the range of integration [1,2], [1,4] and [1,8], the error estimate is too crude so that the ratio of maximum absolute error to i t s estimate.becomes <<1. Also this ratio goes on decreasing as the range of integration increases and therefore the error estimate based on (5.9) does not reflect well the variation in the actual error. Besides the D.E. (5.10) we also tried the following ' TABLE IV EXPERIMENTS WITH PROBLEM (5.10) b h Max. Abs. Error Error Estimate based on (5.9) Ratio 2 1/4 0.18 E-3 0.36 E-2 0.050 1/8 0.13 E-4 0.29 E-3 0.045 • 1/16 0.82 E-6 0.19 E-4 0.043 4 1/4 0.21 E-3 0.56 E-2 0.037 1/8 0.14 E-4 0.46 E-2 0.030 8 1/4 0.38 E-3 0.15 E - l 0.025 two problems: y" + (1 - 2/x2) y = -2 cosx/x 2 ...(5.11a) with boundary conditions such that the true solution of the D.E. turns out to be y(x) = sinx/x. y" + ((3 + 4x)/16x2) y = 0, y(l) = 1, y(b) ...(5.11b) is so chosen that the true solution of the problem turns out to be y(x) = x"*"^ 4 [cos (/x - 1) - .5 sin (/x - 1) ] . The conclusions are identical to those obtained for the problem'(5.10). The experimental results based on theorem 5.1 w i l l be given in the last chapter.' • CHAPTER VI - — NON-LINEAR DIFFERENTIAL EQUATIONS 6.1 If the function f(x,y) introduced i n D.E. (1.1) is not linear in y, one cannot hope to solve the system of equations (2.4) namely J y + h 2 D f(y) = d , by algebraic methods. Some iterative procedure must be adopted. The following discussion is devoted to proving the existence and uniqueness of the solution of the system (2.4). At the same time, this solution w i l l be shown to be an approximation to the solution of the corresponding boundary value problem. The basis for the theorem which we shall state is a method which is a generalization of Newton's method. Kantorovich [16] proved a very general result which guarantees the convergence of Newton's method without even assuming the existence of a solution. We shall state Kantorovich's.result in a form applicable to our discussion. Theorem 6.1 > Let the n equations <j>± (y1,y2»---» y n) = °> ( i = i,2,...,n) for n unknowns y,,y„,...,' y be written in the vector form as <j>(y) = 0 . Let A(y) = (a..) denote the matrix with elements 1 J 8(f). •• , .then the Newton's method ,is written i n the form y(vf« = y(v) _ [ A - i ^ j , ( y ( v > ) f ( v = 0 > 1 > > > o m \ ; ; i 6 a ) Now i f the following conditions are satisfied: (i) For y = y ^ , the i n i t i a l approximation, the matrix A ( y ^ ) has an inverse and an upper bound for i t s norm is known i.e. I|r0l|. < B Q . ...(6.2) ( i i ) The vector y ^ ^ approximately satisfies the system of equations in the sense that the f i r s t correction vector <j> ( y ^ ) is such that 1 |rQ ( y ( 0 ) ) | |. < n0 • ...(6.3) ( i i i ) In the region defined by inequality (6.6) below, the compo-nents of the vector <j>(y) are twice continuously differentiable with respect to the components of y and satisfy n I j,k=l 9 2 <f>. l 3 V y k - K Q ( i = 1,2,..., n) ...(6.4) (iv) The constants B Q , T ) Q and K Q introduced above satisfy the inequality 0 < h Q = B Q n Q K Q < 1/2 , ...(6.5) then the system of equations <Ky) has a solution y* which is located in the hypercube . . 1 - / l - 2h My - y ( 0 ) l l i r—- % . ...(6.6) 0 ( v ) Moreover, the successive approximation y defined above exist and con-verge to y*, and a bound of ||y^-y*|| is given by the inequality My ( v ) - y * l l i e f e r ( 2 h 0 ) 2 V - 1 n 0 ...(6.7). For the proof of the above theorem, see [12], Theorem 7.6. We would like to remind ourselves that at present we are attempting to solve the system of nonlinear algebraic equations (2.4) by Newton's method. Now for a given y, the residual vector r(y) is defined as 71. r(y) = J y + h D f(y) - d . ...(6.8) For our problem, the matrix A(y) introduced in the Theorem 6.1 is given by A(y) = J + h 2 D F (y) , ...(6.9) (the matrices J and D have been introduced in the Section 2.2) . where F(y) = f y(x 2,y 2) and we assume that f (x,y) exists and is bounded in S. Now y ^ ^ i s a yy ' • vector believed to be close to the actual .solution of (6.1) so that the residual vector r ( y ( 0 ) ) = J y ( 0 ) + h 2 D f ( y ( 0 ) ) is small. We now solve the linearized system d' r ( y ( 0 ) ) + ( J + h 2 D F ( y ( 0 ) } ) A y ( 0 ) = Q $ (6.10) and i t s solution i s given by Ay(0> - - [ J + h 2 D F ( y ^ ) ] " 1 r(y<°>> - - [ A ^ ) ] " 1 r<y < 0 )) , provided that the inverse of the matrix A(y) defined by (6.9) exists for y = y ^ . Then the vector y ^ = y ^ + A y ^ w i l l be a better approx-imation to the exact solution, the residual vector r ( y ^ ^ ) w i l l be smaller, and the process can be repeated with y ^ ^ taking the place of y ^ \ etc., u n t i l the convergence is achieved. Note that the solution of the linearized system (6.10) is • greatly f a c i l i t a t e d by the fact that the matrix A(y^^) i s again t r i -diagonal. In fact, i f A ( y ^ ) - (a ), we have w mn " an sn-l = - 1 + ^ y K - l ^ . <n= 2>3,-..,N) a = 2 + h V f (x ,y ( 0 )) (n = 1,2,...,N) ...(6.11) n,n 1 y n y n ' ' ' *n,n+l - - 1 + ^ V y ^ n + l ' O (n = 1,2,. . . ,N-1) and a l l the other elements are zero. Now the method described earlier for the linear boundary value problem in Chapter II is immediately appli-cable. The only extra work that i s to be required is to evaluate the residual vector r ( y ^ ) and the partial derivatives f (x ,y ) . y n y n 6.2 CONVERGENCE OF NEWTON'S METHOD FOR THE SYSTEM (2.4) We shall verify that the conditions of the Theorem 6.1 can be satisfied for the system of the nonlinear equations (2.4) based on the A.E. (2.10b). Since $ ± >_ 0, f >_ 0, the matrix A(y) given by the equa-tion (6.9) has a l l the right kind of properties stated in (2.8). Thus we have 0 < [A(y)]" 1 < J " 1 . But || <_ (N+1) /8, i t follows that • (6.2) i s satisfied for >0 Since B n = (N+l)2/8 = (b-a) 2/8h 2 . ...(6.12) * i ( y ) = " y i - l + 2 y i _ y i + l + ^ 2 ' 1 2 ) [ f ' ( x i - l ' y i - l ) + 1 0 f< xi»yi> + f ( x . + 1 , y . + 1 ) ] , hence I Z 1 = (h /12) [f (x. . ,y. ..) + 10 f (x.,y.) j,k=l 3 y j ^ k yy x - 1 " i - l ' yy v x'/x + f y y ( x i + l ' y i + l ) ] > ( ± = 1 ' 2 - - " N ) 73. or N I j , k - l 324>±(y) 3 y j8y k < h max. f (x,y) yy Thus the condition (6.4) of the theorem 6.1 is satisfied for K Q = h 2 L 2 . (6.13) where L„ = maximum f (x,y) (x,y)eS 1 7 7 .(0) Now let the i n i t i a l approximations y be defined by where Q(x)eC 6 satisfying Q(a) = y , Q(b) = y, , and let Q = max. o/n^ (x) 3. D • n yi R = max. |Q"(x) - f(x,Q(x))| , xc [a,b] . X It then follows (see (2.10b)) that |-Q(x) + 2Q(x+h) - Q(x+2h) + (h2/12) [Q"(x) + 10Q"(x+h) + Q"(x+2h)] (6.14) <_ h° Q6/240 Hence l*^/' 0*)!. - |-Q(x._^ 2Q(x.) - Q(x. + 1) +(h 2/l2) [f(x._ 1,Q(x._ 1)) + 10 f(x.,Q(x.)) + f ( x . + 1 , Q ( x i + 1 ) ) ] | = l-Qtx.^) +2Q(x.) -Q(x. + 1) +(h 2/12)[Q"(x._ 1) +10Q"(x.) +Q" (x. + 1) ] - (h2./12)[(Q,,(x._;L) -f(x._ 1,Q(x._ 1))) +10(Q , ,(x.)-f(x.,Q(x i))) + (Q"(x i + 1) - f ( x ± + 1 , Q ( x i + 1 ) ) ) ] | < h6Q6/240 + (h2/12) (R + 10R + R) < h6Q6/240 + h2R, ( i = 1,2,...,N) . Thus (6.3) i s satisfied for 2 n Q = -^f- (R + h 4 Q6/240) ..(6.15) 74. The main condition (6.5) of the Theorem 6.1, which guarantees convergence of Newton's process, thus turns out to be satisfied i f (b-a) 4 L 2 (R + h 4 Q6/240) < j . ...(6.16a) If as an i n i t i a l approximation Q(x) a polynomial of degree <6 i s chosen, then = 0, and the method w i l l converge i f the quantity R defined by the equation (6.14) satisfies the inequality R< 3-| . ' ...(6.16b) (b-a)* L 2 Having thus demonstrated the existence of a solution of the system (2.4) uniqueness can now be proved. If y* and y** are any two solutions of the system (2.4), viz. J y + h 2 D f(y) = d , then we may put f(x ,y*) - f(x ,y**) = (y* - y**) n , ...(6.17) n y n n ' n . v n . n n 9 f where n is a value of — .• The vector (y* - y**) is easily seen to n 3y J J J satisfy the system (J + h 2 D B) (y* - y**) = 0 , ...(6.18) where J and D are defined by equation (2.4), and B is the diagonal matrix 2 with poiitive elements n . The matrix J + h D B is known to be non-n singular and monotone. It follows that the vector y* - y** = 0 . 6.3 The conditions of convergence of Newton's method namely (6.16a) and (6.16b) can be improved by using a better estimate B^ for the norm of the matrix TQ using techniques of the section 2.4. In fact, one proves that (6.2) is satisfied for where u = (24 + 10 h 2 L)/(12 - h 2 L) > 2 , L = min. 8^^ * y ^ »• (x,y)eS and ^ ( u ) has been introduced by (2.14) Thus.(6.16a) can be rewritten i n the form - f H^(U) (R + _ ^ ) < I , and (6.16b) i n the form ,2 R < 2 L 2 H 2(u) 6.4 Now we w i l l consider the solution of the nonlinear D . E . (1.1) by A.E.'s of higher order. In particular we consider the use of A . E . (3.8) at the interior points and A.E.'s (3.9) near the boundary points. When the two point nonlinear boundary value problem of class M i s replaced by A.E.'s (3.8) and (3.9), we get the system J(c)y + h 2 D(c) f(y) = d ...(6.19) where J(c) and D(c) have been introduced by (3.12). Now we w i l l develop c r i t e r i a for the convergence of Newton's method, as applied to the system (6.19). We w i l l estimate B Q, n Q and K Q as introduced in the Theorem 6.1. For our purpose A(y) = J(c) + h 2 D(c) F(y) , ...(6.20) where F(y) is defined i n (6.9). We also have [A(y)]" 1 ^ J _ 1 ( c ) , provided c, h satisfy the inequalities (3.16) and (3.17). Thus the cond-i t i o n (6.2) i s now satisfied for B 0 = ( B I A ) 2 H f c ) » ...(6.21). 8ti (c-2) W see inequality (3.20). Again using notations of the Theorem 6.1, we have 76. .H I 92 0. 9 y j 9 y k ti L 2 (c+1) , ( i = 1,N) h 2 L 2 (c+2) , ( i = 2,3 N-1) . Thus we have KQ = h L 2 (c+2) . Also following a usual analysis, as i n the previous section (6.22) r ± ( y ( 0 ) ) h 2 R (c+1) + h 6 G 1 Q6 , ( i = 1,N) h 2 R (c+2) + h 8 G 2 Qg , ( i = 2,3,...,N-1). c+1 for some real constants G. and G„. In fact G, = -r-rr and G„ = 1 z 1 240 z 31c-190| 60,480 c ¥ 1 9 0 3 1 , the other symbols having been defined in the previous section. Thus ( 6 . 3 ) would be satisfied for N 0 = 8^?2T V c ) [ R ( G + 2 ) + ° ] . . . (6.23) 4 6 where a = max. [h • G- Q,, h G„ Q„] . Hence the condition for the con-1 o z o vergence of Newton's method in the present case becomes H 0 " (b-a) * (c+2) L 2 [R(c+2) + a] ^ ( c ) 1 ± 2 ' (6.24) 6 4 (c-2) Also as before, i f t h e . i n i t i a l approximation Q(x) is a polynomial of degree <6, then a = 0 and the method w i l l converge i f the quantity R satisfies the inequality R < 32 (c-2)' (b-a) 4 ( c + 2 ) 2 L, ..(6.25) 6.5 We also used a A.E. of order k = 2 obtained by overdifferentia-tion for solving nonlinear boundary value problem of class M. In pa r t i -cular we used the A.E. (4.9). We used auxilliary A.E.'s (4.18) and (4.20) to obtain y^ and y n + 2 respectively. We f i r s t substituted y 3x 3y y 77. using the D.E., and then the expressions for (m = 0,2) in the A.E. (4.9) . We w i l l present numerical results using the above technique as well as the techniques discussed in Sections 6.2 and 6.4, respectively. We might remark that we w i l l not try to use A.E.'s of the form (4.1) for t _> 4 to solve a nonlinear boundary value problem of class M. The main reason is that y ( 4 ) = f + f f + 2 f y ' + (y') 2 f , xx y xy yy therefore when we substitute y", y"1 ...... .etc. using the D.E. and then the expressions for y n + m (m = 0,1,2) obtained from auxilliary A.E.'s (4.18), (4.19) and (4.20) in A.E. (4.1), the resulting A.E. becomes quite g complicated and at the same time i t retains a T.E. = 0(h ), which i s not at a l l an improvement on the previous situation. CHAPTER VII BOUNDARY VALUE TECHNIQUES AS APPLIED TO.INITIAL VALUE PROBLEMS ' 7.1 GENERAL In the present chapter i t is shown that the boundary value •techniques developed to solve the two point boundary problem can be used to solve ia certain class of i n i t i a l value problems which are usually solved by step-by-step method. Allen and Severn [1] suggested a method whereby a f i r s t order system, normally solved by i n i t i a l value techniques, might be replaced by a second order boundary value problem. No indica-tion was given of the merit of this procedure, except when extended to . a certain partial differential equations. Later Fox [10] proposed to treat a derived second order equation; but remarked that either i n i t i a l value or boundary value techniques might be used, the choice depending on the nature of the solution. No great details were given regarding this choice. Fox and Mitchell [11] in 1957 investigated these two methods by which a f i r s t order D.E. can be transformed into a second order D.E., solvable either by i n i t i a l value or boundary value techniques. Their i l l u s t r a t i v e examples were based on linear D.E.'s with constant coeffic-ients . In what follows, we w i l l consider a certain family of f i r s t order linear D.E.'s and solve them by the boundary value techniques, developed i n Chapters II and IV. . It w i l l be proved that the resulting t error is 0(h ), t = 4,6 or 8 . 7.2 LINEAR D.E. OF.THE FIRST ORDER Consider a D.E. of the f i r s t order in the form y' = <Kx,y) , y(a) = y o . ...(7.1) Si Fox [10] produces a second order D.E. by direct differentiation of the original problem (7.1) and subsequent elimination of the f i r s t derivative. For instance (7.1) is replaced by y" - If + £ • For boundary conditions, he takes the given condition y(a) = y and the 3. f i r s t order equation (7.1) applied at x = b viz. y(a) = y a , y'(b) = <Kb,y b) . ...(7.3) On approximating (7.2) - (7.3) by a A.E. of the form (1.2) or (1.3), we get a system of algebraic equations, linear i f <j>(x,y) is linear. If <)>(x,y) is nonlinear, then the resulting system of algebraic equations must be solved by some iterative method. However in this chapter, we want to show how the techniques developed in Chapters II and IV for solving two point linear boundary value problem of class M (M-1, M-4 and M-5) can be modified to solve a linear D.E. of the f i r s t order (7.1) and error bounds can be conveniently obtained provided certain conditions (to be stated later) are satisfied. Thus we assume <j>(x,y) = f(x) y + g(x) , hence the derived second order boundary value problem takes the form ~ y" = [f 2(x) + f*(x)] y + [f(x) g(x) + g'(x)] = F(x) y + G(x) ' o ...(7.4) y(a) = y a , y'(b) = f(b) y(b) + g(b) , a. which we want to. solve by modifying the direct methods M-1, M-4 and M-5, respectively. This time the points x n are given by = a + nh (n = 0,1,2,...,N) with ' h = (b-a)/N . We f i r s t use the A.E. (2.10b) to replace the boundary value problem (7.4). This A.E. is such that the operator associated with i t i s definite and the generalized mean value theorem (1.11) holds, see Section 2.1. We then obtain the following system of (N-1) linear algebraic equations in N unknowns y^, y^,..., y^,' (- 1 + I I V y n + <2 + ^ W VKL + ^ + Tl W yn+2 h 2 + 12 ( G n + 1 0 Gn+1 + Gn+2} = °» ( n = 0,1,...,N-2). ...(7.5) Thus we have (N-1) equations in N unknowns. We develop another equation by considering a A.E. of the form . k k . k T a. y = h Y b. y' . +...+ h J 7 . (j) . ^ . ,_ . L n l 'n+x . L n x 'n+x . L n 9. y ^ . , j. > 1 , ...(7.6) x=0 x=0 x=0 i 'n+i J ' which involves the boundary condition at x = b. The reason for this choice of the above A.E. is that the matrix associated with the resulting system of linear equations becomes tridiagonal for k = 1, j > 1 . The following Table V gives the values of the coefficients associated with the A.E. (7.6) for k = 1 and j =.3,4 and 5. TABLE V j COEFFICIENTS 3 OF THE A.E. (7.6) 4 5 -1 -1 -1 a l 1 1 .1 1/2 1/2 1/2 b l 1/2 1/2 1/2 co 1/10 3/28 1/9 c l -1/10 -3/28 -1/9 do 1/120 1/84 1/72 d l 1/120 1/84 1/72 eo 1/1,680 1/1,008 e l -1/1,680 -1/1,008 1/30,240 g l 1/30,240 T.E. -h 7 y(7Vl00,800 h 9 y(9)/25,401,600 - h 1 1 y ^ m x i ; The coefficients li s t e d above for j = 3 and 4 have also been given by Lambert and Mitchell»[18]. The A.E. (7.6) for k = 1, j = 3 and n = N-1 becomes " y N - l + yN = ( h / 2 ) ( y ' N - l + yN> + ( h 2 / 1 0 ) ( y N - l " yN> ••• ( 7' 7 ) + (h3/120) ( y ^ + y»' ) . The numerical solution of the system (7.5) along with the equation (7.7) is assumed to approximate the solution of the boundary value problem (7.4). 83. 7.3 ERROR ANALYSIS The error equation i s obtained in the usual manner in the form Ae = F(y) , where we now have (7.8) A l B2 B l A2 B3 B2 A3 BN-2 \ - l BN BN-1 \ , r(y) = 'N-1 with n = 2 + (lOh /12) F - l - k + ^ F , . n (n = 1,2,...,N-1) , 3 (n = 1,2,...,N) , 3 **N 1 ~ 2 2 F N + 10 F N " 120 * N B — 1 + ^ - F n 12 n 2 R* = -, _ h f h_ _ n N - 1 6 2 N-1"10 N - 1 120 ^ -ho/6' ^n + V > Y N = IboTsoo y ( ? ) u ) » " V i * N (n = l , 2,...,N-l ) , | 0 n|<l Here y' R(x) y + S(x) . Please note that the A . E . (7.7) also satis-fies the generalized mean value theorem. This can be proved by using techniques discussed in Section 2.1. Now for a l l sufficiently small values of h, we have l | r < y ) " h 6 M, - 240 However, i f f(x) <0 and F(x) >_ 0 over [a,b] (note that these two conditions w i l l automatically be satisfied in case i f f(x) i s a real negative constant), then the matrix A has a l l the right kind of proper-ties mentioned in (2.8) and i t can be shown that 0 < A - 1 < J _ 1 where J differs from J in the last element on the main diagonal which is 1 rather than 2. In fact '2 -1 - l r 2 -1 -1 2 -1 J = -1 2 -1 -1 1 and J ^ can be determined expl 1 1 1 1 2 2 1 2 3 J - 1 3 3 c i t l y (see Rutherford [23]) in the form 1 2 3 N-1 N-1 N-1 N Now we prove-the following lemma. Lemma 7.1 Let Jt = co, where t i s a vector with components t^, t^y co being a column vector with each component unity, then (i) ( i i ) t. = x i(2N+l-i) |t| = N(N+1) 2 85. The proof of this lemma is t r i v i a l , since t. is just the sum x of the elements of the i t h row in J ^ . Now using (7.8), we get e = A" 1 r(y) h 6 M I | 6 i(2N+l-i) • l e J 1-240- * 2 • h 6 M,.N(N+1) •and Hell <- 6 480 h 4 M f i 2 l l e l l 1^37/ [(b-a) Z + h(b-a)] . ...(7.9) Thus we have proved the following theorem: Theorem' 7.2 Given the boundary value problem "y" = [f 2(x) + f ( x ) ] . y + [f(x) g(x) + g*(x)] = F(x) y + G(x) , 7(a) = y . y'(b) = f(b) y(b) + g(b) 3. (derived from the i n i t i a l value problem y' = f(x) y + g(x), y(a) = y„) , we replace i t by the A.E.'s y n " 2 y n + l + yn +2 = K + 10K+1 + K+2> (n = 0,1,...,N-2)and -v + v = — (v* +v') + - — (v" - v") + - — (v 1" + v'" ) a n Q yN-l ' yN 2 ^ yN-l yN ; 10 ^ yN-l yW 120 ^ yN-l yN ' 2 3 h , , . , s . h , .. ... . h Then i f (i) f(x) <_ 0, F(x) _> 0 over [a,b] ( i i ) y(x) e C 6, f(x), g(x) e C 3 , h 4 M I hi I l ~ 4 8 o ^ t(b-a) Z + h(b-a) ] . or equivalently ||e|I = 0(h 4) . 86. The boundary value technique described above to solve an i n i t i a l value problem of the f i r s t order where the error estimate i s given by Theorem 7.2 w i l l henceforth be referred to as M-1 . 7.4 AN 0(h 6) AND 0(h 8) FINITE DIFFERENCE ANALOGUE FOR (7.4) We approximate the boundary value problem (7.4) by the A.E. (4.9), viz. y - 2y ,. + y ,. = h 2 £ 3. y" + h 3 E y. y"' (n = 0,1,...,N-2) •'n •'n+l n^+2 l •'n+x x •'n+x ' ' ' and 2 3 -%-i + % =! + y$ + % ( y N-i - yS> + h + y? + _ J l 4 _ . . ( y ( 4 ) _ ( 4 ) , _ 1,680 k y N - l yN ' We prove in a manner analogous to the one used in previous section that 29 h M l l e M < 6 Q 4 > 8 0 0 8 [(b-a) 2 + h(b-a)] .. " ...(7.10) This method will.be designated as M-4 for reference purposes. 8 In order to develop an 0(h ) f i n i t e difference analogue to solve (7.4) we use the A.E. (4.24) viz. 2 4 (W) y - 2y ., + y 0 = h £ 3. y". + h E y. y\J. (n = 0,1 N-2) J n Jn+1 Jn+2 l y n+x x Jn+x ' and 2 3 " yN-l + yN = ! ( y N - l + * + f ( y N - l " ^ + 71 ( y N - l + yN > + ( y ( 4 ) - y ( 4 ) ) + ^ ( / 5 ) + y ( 5 ) ) . 1,008 > yN-l yN ; 30,240 v yN-l. yN ' Under the same conditions as mentioned in Theorem 7.2, we prove that 59 h 8 M H£M ^ 152,409,600 t(b-a)2 +h(b-a)] . ....(7.11) 87. 0/ This method w i l l henceforth be referred to as M-5. Experimental results to be given in the last chapter, w i l l reveal that the method M-5 is best so far to solve (7.4). 7.5 THE CASE WHEN f(x) g f (constant) The error, bounds (for the boundary value problem (7.4)) given by (7.9), (7.10) and (7.11) are obtained on the assumption that 'f(x) <_ 0, F(x) _> 0 over [a,b] . We have also remarked that the above conditions are automatically satisfied i f f(x) is a real negative constant. We w i l l now show that error bounds of the form (7.9) and (7.11) could be obtained i f the function f(x)'is constant = f (say) over [a,b] (regardless of the sign of f ) . The matrix A introduced in (7.8) could now be expressed as a product of P^ (u) and Q. The matrix Q has already been introduced in Section 2.4 and V - l . -1 y -1 P N(y) = - l -a -1 v b ..(7.12) Here X y = a = and b .= = 1 - h f 2-< 1 (2 + h 2 f 2)/X > 2 , h , . h 2 2 - h 3 ,3\ ( 1 " 2 f + 10 f - 120 f ) A ...(7.13) 88 Let denote the determinant of the nth principal minor of the matrix P*(y) , (n = 1,2,...N-1), then D n satisfy the A . E . (2.13a). Hence _ sinh(n+1)9 ' D = •' -un ' » y = 2 cosh9 > 2 , n smh9 ' (n = 0,1,...,N-1), see (2.13b) . * Also, i f <J>N denotes the determinant of the matrix P N(y), then *N " b DN-1 " a DN-2 Let P'' "'"(y) = (p..), then we have N r i j D . , <j> . / ( j ) , j < i j-1 yn-j' yn J -P i , N , j > i ( i , j = 1.2,. a D , (j = 1,2,. , ( i = 1,2, •,N-1) .,N-1) .,N) .,. • (7.15) The above results can be obtained by following the usual analysis used in Section 2.4, for proving the corresponding results for the matrix PAT(y) . Now let R. denote the sum of the absolute values of the elements N l of the i t h row in the matrix P^ (u) i.e. N R i = X | p i j ' ' j=l (note that the inverse of the matrix PJJ(P) 1 S n o t necessarily a monotone matrix unless f <_ 0) , then we have 89. N-1 R l = . j * N - j=0 R. = N-i 'N -L. I D- 1 + i - l N 'N' j=i+l ( i = 2,3,...,N-1), <J,0 = 0 , RN= ( a T D j + D N - 1 ) 7 I * N I 3=0 ...(7.16) We can also prove using (2.13e) that * (D N_1 " DN-2 " « + ^"2> °N-1 (P-2) |*NI * - l , We f i n a l l y obtain the norm of the matrix P^ (u) v i z . ||P* _ 1(u)|| = max. (R A ) = R (say) . . ... l<i<N The quantity R can only be computed numerically, using the expressions for R^ ( i = 1,2,...,N) given by (7.16). Now using (7.8), we obtain (7.17) < h 6 M£ R/240 — 6 ...(7.18) The boundary value technique described in Theorem 7.2 to solve an i n i t i a l value problem of the f i r s t order where the error estimate is given by (7.18) w i l l henceforth be referred to as M-1.' Numerical results given in the next chapter reveal that the error estimate based on M -1 is better than the onebased on M-1 . CHAPTER VIII EXPERIMENTAL RESULTS AND CONCLUDING REMARKS 8.1 LINEAR BOUNDARY VALUE PROBLEMS OF CLASS M In this section, we consider linear boundary value problems of class M. The methods M-1, M-2(c), M-3(c), M-4, M-5 and M-6 introduced earlier w i l l be used for obtaining numerical solutions and theoretical error bounds. M-6 w i l l only be used to solve linear problems with con-stant coefficients. The following D.E.'s are chosen to present the num-erica l results in order to study the comparative merits of methods from M-1 to M-5. (A) (B) x 2y" with x 2y" - 2y = y(x) -= 3 - V = -x, y(2) = y(3) = 0 1 / - i n c 2 36,. 38 ( 1 9 x " 5 x ~Ty x 2, y(l) = -2.2, y(2) = with y(x) = 0.8 x 2 + (x 2-4)//x . (C) y" - (l+x 2) y = 0, y(0) = 1, y(2) = with y(x) = exp (x2/2) (D) y" - 2 y = x • -2, y(0) = 0, y(l) = 1 with y(x) = 2 (2 sinh x/sinh 1) - x (E) y" - y = -4 x e x, y(0) = y(l) = 0 with y(x) = x(l-x) e X (F) y" + y = -x , y(0) = y(b) = 0 with y(x) = (b sinx/sinb) - x . Please note that the problem (F) r-is not of class M, but has been included for the sake of comparison. We now discuss problem (A) in detail. The true solution of (A) is used to check the accuracy of the numerical solution. We also have (n) . 18 , T Vn+l , , n+1 . y (x) = (-1) n!/x , n >_ 3 . Thus Mn = 18xn!/(19><2n+1) , n >_ 3 . Also experimentally i t has been observed that the optimum choice of 'c' in M-2(c) and M-3(c) depends on 'h' to a very great extent. Sometimes the behaviour of this dependence varies from problem to problem. In some cases i t has been observed that optimum 'c' remains constant while in others i t either decreases or i n -creases with decreasing 'h', t i l l we - h i t the round-off region. Table VI for instance, displays how the optimum choice of ' c" varies with 'h' i n method M-3(c)» l o g 2 (l/2h) TABLE VI OPTIMUM CHOICE OF c IN M-3(c) = c m c in problem m (A) (B) (C) (D) (E) (F) 2 3 4 5 6 7 4.1 4 4 4 5.9 5.3 4.6 3.9 5.8 6.3 4 4 4 • 4 4 4 4 4 4 3 Round off Region 4 4 4 4 4 4 4 5 The computational results are summarized in Table VII for problem (A). TABLE VII EXPERIMENTS WITH PROBLEM (A) WITH h=l/8 Method Max. Abs. Error Estimated Error Ratio M-1 0.17 E-6 0.66 E-6 0.27 M-2(5) 0.49 E-7 0.14 E-5 0.04 M-3(4) 0.97 E-9 0.19 E-7 0.05 M-4 0.65 E-9 0.34 E-8 0.19 M-5 0.15 E - l l 0.11 E-10 0.14 These computations were performed on IBM 7040 at the University of British Columbia, Vancouver. . In order to reduce round-off error to a minimum, a l l calculations were made using double precision arithmetic. The tabulated results in this chapter were a l l based on double precision arithmetic. We plotted curves of log^Q|Error|versus log2(l/2h) and also log^Q|Error|versus log^ Cost under two different assumptions given below by (8.1). The curves which we have plotted to study the variation of log^|Error|versus log^ (l/2h) (or log^^|Error|versus l o g 2 Cost) are based on methods M-1, M-3(c) and M-5 respectively. Also in order to estimate the cost of computation, the. following two different assumptions are made. 8.1 (i) Assume that the time taken by the machine in calculating the elements of the associated matrix and solving the system of linear algebraic equations is negligible as compared with the time taken i n evaluating g(x) , s(x) and'their higher derivatives at the grid points x (n = 0,1,...,N+1) . or 8.1 ( i i ) Assume that the time taken by the machine i n evaluating g(x), s(x) and their derivatives at the grid points is negigible as com-pared with the time taken in evaluating the elements of the associated matrix and solving the system of linear algebraic equations. The accompanying graphs show the variation of l o g ^ l Error | versus log^ (l/2h) for the problems (A) and (D), (see graphs 1 and 7, respectively). The graphs of log^Q]Error|versus l o g 2 Cost are also given for the same problems. In this case, the only difference is that the graphs of methods M-3(c) and M-5 are shifted to the l e f t . The graphs for the other problems are similar to those of (A) and (D). These curves g indicate that in a l l cases that we tried, the method M-5 with e = 0(h ) is superior to the rest of the techniques. Thus in comparison to M-1, M-2(c), M-3(c) or M-4 we w i l l recommend M-5 to solve a linear boundary value problem provided analytic expressions are given for functions g(x) and s(x), so that the required higher derivatives of these functions at grid points may be computed easily. But i f we do not have analytic expressions for g(x) and s(x), then the only choice of method w i l l be M-1 or M-2(c). For the sake of comparison, we also solved the above problems using single precision arithmetic. We have also plotted the curves for the problem (A), (see Figure 4, 5 and 6 respectively). The conclusions are slightly different from those obtained above using double precision arithmetic. Here the curves log1„|Error|versus log_ (l/2h) show (see -10 - I I 1x1 o -12 -13 •14 -15 •H6I J L FIGURE 2 Log 2 Cost J L Problem (A) Assumption 8.1(i) J L f £ + \ fl+2 fJL+4 {JL+5 LL+6 LL+7 fJL+8 fJL+9 IJOJJ3 Bo] 97. - 5 FIGURE 4 Problem (A) (Single Precision) .4 rf A A-6 o o -8 u \ \ \ _M- I \ \ ' ' X T * ' M-3(4) A - 5 3 4 log«(l/2h) * FIGURE 5 Assumption 8.1(i) 8 §-6f UJ o -8h \ \ \ \ q \ A \ ,\\ 1-1 • 3 ( 4 ) ^ . X T , ~ o — -.XT' / JST' M - 5 ^ o -• \ - V -9 /ir!-4 / x + 5 Log2 Cost -/£+6 fA+7 FIGURE 6 Problem (A) . Assumption ( i i ) ^ ~ (Single Precision) i — • 1 1 1 1 1 : ! i L _ /i+0 fj.+\ /x+5 /i+6 /i+7 Log 2 Cost > 00 102. Figure 4 for problem (A)) that as h decreases, a l l methods become equally accurate. M-5 has no advantage either over M-3(c) or M-1. We hit the round-off region even before reaching h = 1/8, see Figure 4. The minimum in a l l the three curves i s attained for values of log2 (l/2h) <_ 3 . However, the curves of log^Q|Error|versus l o g 2 Cost (see Figures 5 and 6) indicate that M-5 is better than M-1 or M-3(c) under either of the two assumptions. It should be noted again that i f func-tions g(x) and s(x) are given i n tabular form, then the only choice of method w i l l be either M-1 or M-2(c). In both these methods, the result-4 ing error i s 0(h ). Since M-1 is easy to use, so i t i s preferred in comparison to M-2(c). the number of function evaluations i s 2N+4, 2N+12 and 6N+12 in methods M-1, M-3(c) and M-5 respectively. In M-1 we evaluate g(x), s(x) at x^ (n = 0,1,.. . ,N+1) ; in M-3(c) we evaluate g(x), s(x) at x^ ' .. (n = 0,1,...,N+1) and g'(x), s'(x) at the points x.. (j = 0,1,N,N+1); fi n a l l y in M-5 we evaluate g(x), g'(x), g"(x), s(x), s'(x) and s"(x) at a l l the N+2 grid points. 8.2 LINEAR PROBLEMS OF CLASS M WITH CONSTANT COEFFICIENT We choose problem (D) for i l l u s t r a t i o n , since the true solu-tion of (D) is Before we close this section, we would like to mention that y(x) (2 sinh x/sinh 1) -x 2 hence 2 sinh x/sinh 1 , m even and > 2 (x) 2 cosh x/sinh 1 , m odd and > 2 . Thus 2, n even and > 2 M n 2 coth 1 = 2.626, n odd and > 2 over the interval [0,1] . Table VIII shows the ratio between the maximum absolute error in num-eri c a l solution to i t s theoretical estimate. TABLE VIII EXPERIMENTS WITH PROBLEM (D) WITH h = 1/8 Method Maximum Abs. Error Estimated Error Ratio M-1 M-2(3) M-3(4) M-4 M-5 0.11 E-6 0.23 E-7 0.66 E-10 0.38 E-10 0.50 E-14 0.23 E-6 0.97 E-6 0.51 E-9 0.91 E-10 0.12 E-13 0.47 0.02 0.13 0.42 0.44 The conclusions are obvious enough that in case of two point linear boundary value problem of class M with constant coefficients, we w i l l recommend M-5 (see Figures 7, 8 and 9). In fact M-6 is superior to M-5 and yields very good results even for large values of step-size 12 8 h. The resulting error i n M-6 is 0(h ) while in M-5 i t i s 0(h ). We solved the problem (D) by taking the other boundary at the point b = 1, 2,4,8 and 16 with h = 1/2 i n each case using methods M-5 and M-6 respec-tively. The maximum absolute errors are tabulated in Table IX. The figures listed do" justify the claim that method M-6 yields better results than M-5 for large step-size and over greater range of integration. 104. TABLE IX EXPERIMENTS WITH PROBLEM (D) WITH h = 1/2 USING M-5 and M-6 b Maximum Absolute Error in Method M-5 M-6 1 0.30 E-9 0.39 E-15 2 0.23 E-8 0.49 E-14 4 0.25 E-7 .0.54 E-13 8 0.14 E-5 0.27 E - l l 16 0.41 E-2 0.82 E- 8 8.3 We w i l l now compare methods M-AH (developed by Aziz and Hubbard [2] and referred earlier in Section 2.4) and M-BH (developed by Bramble and Hubbard mentioned in Section 3.4) with our techniques M-1, M-2(c), M-3(c) and M-5. In M-AH the D.E. (2.21) i s replaced by a set of n sim-ultaneous linear equations 2 2 2 4 ( 1 " 12 8 n ) y n + ( 2 + I12"' gn+l ) yn+l + (- 1 + 12 gn+2)yn+2 = h * Si+1 + 12 Si+1> (n = 0,1,...,N-1) . 4 The resulting error i s 0(h ) as mentioned earlier, (2.22). Now we tabulate maximum absolute errors in the entire range of integration for the problem (B) for a series of values of h and x^ +2 =b. A glance at Table -X shows that method M-5' should be preferred over M-1, M-AH, M-BH, M-2(c) or M-3(c). Even M-3(c) is better than M-1, M-AH, M-BH and M-2(c). The methods M-1 and M-AH become identical i f s(x) i s a constant. That is why in Table X the entries in columns M-1 and M-AH are the same. 105. However M-AH gives better results than M-1 i f s(x) is not a constant. In Table XI we display the results obtained by solving the problem (E) for a series of values of b and h using methods M-1 and- M-AH. TABLE X EXPERIMENTS WITH PROBLEM (B) VIZ. x 2y" - j y = x 2, y(l) = 2.2, y(b) = .8 b 2 + (b2-4)//b~ WITH y(x) =' .8 x 2 + (x2-4)//x" . b h Maximum Abs. Error in Method - M--1 M--AH M--BH M-2(c) M-3(c) M--5 2 1/4 0.15 E-3 0.15 E-3 0.18 E-6 1/8 0.10 E-4 0.10 E-4 0.26 E-2 0.42 E-5 0.16 E-6 0.83 E-9 1/16 0.65 E-6 0.65 E-6 0.27 E-3 0.10 E-6 0.31 E-8 0.34 E-11 1/32 0.40 E-7 . 0.40 E-7 0.23 E-4 0.22 E-8 0.53 E-10 0.18 E-13 3 1/4 0.22 E-3 0.22 E-3 0.18 E- l 0.14 E-3 0.85 E-5 0.21 E-6 1/8 0.14 E-4 0.14 E-4 0.27 E-2 0.44 E-5 0.21 E-6 0.99 E-9 1/16 0.89 E-6 0.89 E-6 0.28 E-3 0.11 E-6 0.40 E-8 0.41 E-11 1/32 0.56 E-7 0.56 E-7 0.24 E-4 0.86 E-8 0.68 E-10 0.11 E-12 5 1/4 0.25 E-3 0.25 E-3 0.19 E - l 0.14 E-3 0.97 E-5 0.23, E-6 1/8 0.16 E-4 0.16 E-4 0.27 E-2 0.45 E-5 0.24 E-6 0.11 E-8 1/16 0.10 E-5 0.10 E-5 0.29 E-3 0.11 E-6 0.45 E-8 0.44 E-11 1/32 0.65 E-7 0.65 E-7 0.24 E-4 0.35 E-7 0.76 E-10 0.13 E-11 9 1/4 0.27 E-3 0.27 E-3 0.19 E- l 0.14 E-3 0.10 E-4 0.24 E-6 1/8 0.17 E-4 0.17 E-4 0.28 E-2 0.45 E-5 0.25 E-6 0.11 E-8 1/16 0.11 E-5 0.11 E-5 0.29 E-3 0.16 E-6 0.47 E-8 0.43 E-11 1/32 0.69. E-7 0.69 E-7 0.24 E-4 0.13 E-6 0.78 E-10 0.19 E-10 17 1/4 0.27 E-3 0.27 E-3 0.19 E- l 0.14 E-3 0.10 E-4 0.24 E-6 1/8 0.18 E-4 0.18 E-4 0.28 E-2 0.45 E-5 0.25 E-6 0.11 E-8 1/16 • 0.11 E-5 0.11 E-5 0.29 E-3 0.51 E-6 0.48 E-8 0.60 E-10 1/32 0.71 E-7 0.71 E-7 0.24 E-4 0.49 E-6 0.12 E-9 0.18 E-9 NOTE: The '-' i n Table X indicates that the method i s not applicable for that value of h. 106. TABLE XI EXPERIMENTS WITH PROBLEM (E-) VIZ. y" - y = -4 xe X, y(0) = 0, y(b) = b (1-b) e b WITH y(x) = x(l-x)e X . Maximum Abs. Error in Methods b h" ""' M--1 M--AH 1 1/4 0.93 E-4 '0.12 E-5 1/8 0.59 E-5 . .. .- 0.81 E-7 1/16 0.37 E-6 . 0.53 E-8 i 1/32 0.23 E-7 0.33 E-9 2 1/4 0.69 E-3 0.82 E-4 . 1/8 0.43 E-4 0.52 E-5 1/16 0.27 E-5 0.33 E-6 1/32 0.17 E-6 0.20 E-7. —• 4 1/4 . 0.11 E- l . - ~ 0.33 E-2 1/8 0.68 E-3 0.21 E-3 1/16 0.42 E-4 0.13 E-4 1/22 0.27 E-5 0.80 E-6 8.4 We w i l l now consider boundary value problems which are not of class M. In particular, we consider y" + E 2y = sinx, (E 2 <f> 0,1 and E + ^ 2L • ...(8.2) where m is an integer) y(0) =1 y(b) = sin(Eb) .+ cos(Eb) + S i ^ b ( E - l ) 107. with y(x) = sin (Ex) + cos (Ex) + sinx/(E -1) . It is easy to verify that Mn <_ /2 E n + 1/(E2-1) , over [a,b] . We solve (8.2) for E = 2,3,4,5 and 6; b = 1,2,4,8 and 16 with h = 1/4, 1/8, 1/16 and 1/32 respectively. We use methods M-1 and M-5. The error estimates are based on Theorem 5.1. The results of D.E. (8.2) for E = 2 and h = 1/8 are given in Table XII using method M-1, and for h = 1/4 using method M-5. Method TABLE XII EXPERIMENTS WITH PROBLEM (8.2) b Max. Abs. Error Estimated Error Ratio M-1 1/8 1 0.18 E-4 0.45 E-4 0.40 2 0.28 E-4 0.12 E-3. 0.24 4 0.64 E-4 0.18 E-3 0.35 8 0.44 E-3 0.13 E-2 0.34 16 0.46 E-3 0.13 E-2 0.35 M-5 1/4 1 0.37 E-8 0.73 E-8 0.50 2 0.70 E-8 0.20 E-7 0.35 4 0.11 E-7 0.33 E-7 0.33 8 0.88 E-7 0.24 E-6 0.37 16 0.86 E-7 0.25 E-6 0.34 108. Method M-5 gave better results as was expected even for large values of the step size h, therefore we have tabulated values of the max. abs. error in case of method M-5 for h = 1/4, E = 2. Our conclusion i s that the formula (5.6) gives a useful e s t i -mate of the discretization error. The ratio of "max. abs. error" to the estimated error i s less than 1. Also this ratio stays f a i r l y constant as b increases, so that (5.6) can be regarded to reflect well the varia-tion in the actual error (even in large range of integration). We have 2 2 not tabulated y= (24 - lOh g)/(12 + h g) because i t is'not sensitive to changes in 'h'. Also u does not depend on 'b' at a l l , hence we are sat-i s f i e d with tabulating b and h in Table XII. Thus i t i s clear that the inequality (5.6) for the estimation of ||e|| in method M-1 is quite useful, and assists one in the selection of 'h' and we can rely upon i t for a f a i r l y large range of integration. 8.5 NONLINEAR PROBLEMS Iterative procedures based on A.E.'s referred to in methods >Hy M-2(c) and M-4 w i l l be studied for solving nonlinear boundary value problems of class M. The following D.E. (8.3) i s chosen for experiment-ation. y" = yy 2, y(0) = 4, y(i) = l , y > o . ...(8.3) The closed solutions y(x,y) of D.E. (8.3) are not known except for y = 0 and u= 3/2. We notice that y(x,0) = 4-3x and y(x,3/2) = 4/(l+x) 2 . The solution y(x,0) i s taken to be an i n i t i a l approximation i.e. Q(x) = 4-3x according to the notations introduced i n Chapter VI. We solve 109. problem (8.3) for p = 3/2 . The criterion for stopping iterations i s that the remainder as defined by (6.8) be such that I|r(y ( j ))|| < 1 0 - 1 0 , for j = 1,2,....N . The results are tabulated i n Table XIII for h = 1/10. In method M-2(c) results are given for c = 3, although we tried the values of c in the range -5 c <_ 18 and varied c each time by an amount of .1. The best results were obtained for c = 3. The results based on M-2(3) are more or less identical to those based on M-1. Now we solve the same problem using method M-4, as described bri e f l y i n Section 6.5. The results given i n Table XIII also confirm that the order of the convergence of Newton's method is 2 i.e. the method is quadratically convergent or that the number of correct decimal places is doubled at each step. The figures given also indicate that the results are best using method M-4 in the sense that "max. abs. error" is least in that case. We notice that the maximum absolute error attains i t s minimum value and then does not change, however the quantity ||r(y)|| , the norm of the remainder vector goes on decreasing as the number of iterations increases, t i l l we hit the round-off region. In practical cases however, we w i l l have no way of obtaining "max. abs. error", we have to stop iterations when ||r(y)|| becomes less than a preassigned small number e>0. Finally we notice that M-2(c) has no advantage over M-1. In M-2(c), we have to decide the optimum value of the parameter c and then the amount of computation in M-2(c) i s more than M-1. Therefore the method M-2(c) i s dropped in comparison to M-1. Between M-1 and M-4, M-4 should be used to solve a nonlinear boundary value problem of class M, because "max. abs. error" is smaller in that case, although the numer of iterations required in both cases i s 4 before ||r(y)|| _< lO -"^ . The only advantage which M-1 has over M-4 is that i t i s very easy to apply. Besides the problem (8.3), we also tried the problem y" = -2 + ysinh y , y(0) = y(l) = 0 ,. y > 0 . The closed solution y(x,y) of the above problem is not known except for y = 0 and y(x,0) = x(l-x). This solution was taken to be an i n i t i a l approximation to the actual solution of the-problem. After only two steps the iterations were stopped in this case. The conclusions drawn were identical to those drawn for the problem (8.3). TABLE XIII EXPERIMENTS WITH NONLINEAR PROBLEM.(8.3), h = .1; y = 3/2 Method No. of Iterations || r(y)|| ...Max. Abs.. Error M-1 0 0.21 1 • 0.76 E-2 0.43 E - l 2 0.27 E-4 0.11 E-3 3 0.35 E-9 0.64 E-4 4 0.46 E-15 0.64 E-4 M-2(3) 0 0.87 1 0.37 E-l 0.43 E - l 2 0.13 E-3 0.15 E-3 3 0.17 E-8 0.32 E-4 4 0.12 E-15 0.32 E-4 M-4 0 0.20 1 0.75 E-2 0.42 E - l • 2 0.26 E-4 0.15 E-3 3 0.31 E-9 0.75 E-6 4 0.16 E-14 0.75 E-6 8.6 INITIAL VALUE PROBLEMS We f i n a l l y present numerical results of i n i t i a l value problems solved by boundary value techniques. The numerical techniques employed are M-1, M-4 and M-5 (introduced in Chapter VII). Most of the experi-ments we performed in this connection were based on D.E. y' = sy + sin(tx) , y ( 0 ) = — r - ^ — , ...(8.4) ( s Z + 0 2 2 with y(x) = -[s sin(tx) + t cos(tx)]/(s + t ). We also have Mn = | t | n / / s Z + x} . We chose h = 3/32, and the "max. abs. error" given in Table XCV is over the range 0 _< x <_ 7.5 for s > 0 and over the range 0 <_ x <_ 21 for s < 0 . The results are listed in Table "XIV below: TABLE XIV EXPERIMENTS WITH PROBLEM- (8.4) Method s t Max. Abs. Error Estimated Error . Ratio x 100 M-1 • -2 1 0.29 E-7 0.32 E-4 0.01 M-4 0.59 E-11 0.64 E-8 0.01 M-5 0.31 E-14 0.45 E-12 0.68 M-1- -1 2 0.27 E-5 0.20 E-2 0.13 ^ . M-4 0.22 E-8 0.16 E-5 0.13 M-5 0.63 E-12 . • 0.46 E-9 0.13 M-1 2 1 ' 0.13 E - l -^ . M-4 0.32 E-2 - -M-5 0.26 E-2 - -M-1 1 2 3.725 - -^ , M-4 0.34 E-2 • • -M-5 0.10 E-5 — . — Note when s > 0, the condition of the Theorem 7.2 viz. f(x) _< 0, is not satisfied and the expressions for error bounds no longer hold good. That is why in Table XIV the corresponding columns for estimated error and ratio have been l e f t blank. We notice that the max. abs. error in M-1 (being less accurate) for s = 1, t = 2 with h = 3/32 over the range [0,7.5] is quite large so that the corespond-ing numerical results are of no practical importance. However the max. abs. error equals 0.32 E-3 over the range [0,3]. This d i f f i c u l t y can also be overcome by using more accurate techniques v i z . M-4 or M-5, as is apparent from Table XIV. We mention that we have obtained results for a number of other situations besides those lis t e d in Table XIV. In addition to D.E. (8.4) we also experimented with the following D.E.'s: y» = S y + e S X , y(0) = 1 with y(x) = (1+x) e S X , ...(8.5) y' = 12y - 11 e X, y(0) = 1 with y(x) = e X, ...(8.6) xy' = -y + x 3, y(l) with y(x) = (x 3 + x _ 1)/4. ...(8.7) We now <use method M -1 and M -5 to get the error bounds for the same problem viz. (8.4). We tabulate the numerical results in Table XV. We have obtained error bounds for f(x) = constant (regardless of the • ... sign of the function) and also the error estimates are better than those given in Table XIV. Note that the analysis of the Section-7.5 i s not applicable to the method M-4 because the A.E. Ea. y = h 2 E3. y" .+ h 3 E y. y"' I Jn+i I Jn+1 l n+i i s such that YQ = 1/40, y^ = 0, y^ = -1/40 , hence again the associated matrix A cannot be factored as a product of PJJ( u) and Q» 113. TABLE XV EXPERIMENTS WITH PROBLEM (8.A), USING M*-l and M*-5 WITH h = 3/32. Method s t Max. Abs. Error Estimated Error Ratio * M -1 -2 1 0.29 E-7 0.32 E-5 .009 * M -5 0.31 E-14 0.46 E-13 .067 M*-l -1 2 0.27 E-5 0.80 E-3 .003 M*-5 0.63 E-12 0.18 E-9 .0035 M*-l 2' 1 0.13 E - l 0.24 .055 M*-5 0.26 E-2 0.25 E - l .104 M*-l 1 2 3.725 0.19 E+3 .002 M*-5 0.19 E-5 0.54 E-4 .002 Now we proceed to compare the methods M-1, M-4 and M-5 with some known techniques usually employed for solving (7.1). For the numerical solution of a l i n e a r . f i r s t order D.E. of the form y' = f(x) y + g(x) , y(a) = y , Allen and Severn introduce an auxilliary function z(x), connected with y by the equation y = ij)(z,z') . The form of jj;(z,z') i s to some extent arbitrary, but is generally chosen to be a simple linear combination of z and z' such that the substitution y = \j>(z,z') i n the above linear i n i t i a l value problem produces a second order D.E. for z i n which the coefficient of z' vanishes. This second order D.E. i s solved with boundary conditions y0 = ^' z0' z0^ a n d Zn = X where X is an arbitrary value of z at any other point x^ in the range. The value of X does not affect the resulting solution for y, which is f i n a l l y recovered from y = i|>(z,z'). As mentioned earlier in Section 7.1 that no indication was given by Allen and Severn regarding the choice of ip(z,z'), except when applied to certain partial D.E.'s. Fox's method is direct and more convenient than that of Allen and Severn. He approximates the derived second order system (7.4) by the A.E. (2.10b) (n = 0,l,...,N-2) along with the A.E. which involves the boundary conditions at x = b - y N - i + % + i • . ( h / 3 ) ( y ' N - i + 4 y N + W • • • ( 8 - 8 ) with a T.E. = (-1/90) h 5 y ( 5 ) . Obviously the A.E. (7.7) which we use to approximate the second boundary condition of (7.4) at x = b i n M-1 is more accurate than A.E. (8.8). With the result M-1 turns out to be better than Fox's method. D.E. (8.6) was considered by Fox [10] to demonstrate the fact that i n i t i a l value problems could be conveniently solved by boundary value techniques. The maximum absolute error for h = .2 was found to be 0.10 E-3. However in our case, the max. abs. error turns out to be 0.71 E-5, 0.45 E-7 and 0.39 E-9 using M-1, M-4 and M-5 respectively. Also in Fox s method, no a. error bounds are available, hence we drop i t in favour of M-1. Now we w i l l compare the results of the problem (8.5) for s = -1 with h = 0.1. using some of the standard i n i t i a l value techniques. 1 . . The methods li s t e d in Table XVI and referred to as w^Cc) f° r the solution of linear D.E. of the f i r s t order are based on / ( 1 - h e k W • • • • ( 8 - 9 ) Thus we need k starting values viz. y , y y ,, -, only. If we n yn+l •'n+k-l assume that the D.E. is such that the functions f(x) and g(x) are quite complicated, then the total cost of computation w i l l be proportional to the (2N+2) evaluations of functions f(x) and g(x) respectively in method w^(c). However the example (8.5) whose results we summarize i n Table XVI is such that we require only (N+1) evaluations of the exponential func-tLon e in a l l the methods listed therein, except in the fourth order. X Runga-Kutta method, where a total of 4N evaluations of e are needed. A close study of the Table XVI shows that the results obtained by using M-1 are better than step-by-step methods w^(.6) and W2(l), but the fourth order Runga-Kutta method gives the best results. However, M-4 yields better results than w^(l). The methods that we propose viz.. . M-1, M-4 and M-5 are direct methods while the remaining are step-by-step methods (or iterative in the sense that they are predictor-corrector methods provided we are solving a nonlinear differential system of the f i r s t order) and suffer from a disadvantage that they require formulas for starting the solution. Thus, in short, the methods that we propse compare very favourably with the standard i n i t i a l value techniques. They are better in most cases and are very easy to apply. 'n+k k-1 k -I (a.-hS.f ..) y . . +' h" 7 S.g . L n i i n + i n + i . L n i i=0 i=0 n+i. 116. TABLE XVI COMPARISON OF M-1, M-4 and M-5 WITH STANDARD INITIAL VALUE TECHNIQUES, PROBLEM (8.5) Order of the • .. Maximum No. of function No. Method: Resulting .Error Absolute '_. Evaluations in the Method Error 1 w^(.6), see Hull and Newberry [15] 0(h 4) ' .81 E-6 2N+2 2 w 2 ( l ) , Milne method 0(h 4) .82 E-6 2N+2 3 M-1 o(h 4) .22 E-6 2N+10 4 Fourth order Runga-Kutta method [13,pp.237] 0(h 4) .18 E-6 4N 5 w^(0), Adams type 0(h 5) .25 E-6 2N+2 6 w 4(l),.[13,(6.6.10)] 0(h 6) .53 E-8 2N+2 7 M-4 0(h 6) .73 E-10 4N+4 8 Method based on A.E. (7.6), j = 4 o(h 8) .93 E-10 3N+3 9 M-5 0(h 8) .81 E-14 4N+4 117. BIBLIOGRAPHY 1. D.N. ALLEN and R.T. SEVERN, The Application of Relaxation Methods to • the Solution of Differential Equations in Three Dimensions; Quart. J. Mech. App. Math. 4 (1951),. pp. 199-208. 2. A.K. AZIZ and B.E. HUBBARD, Bounds for the Solution of the Sturm-Liouville Problem with Application to Finite Difference Methods; J. SIAM, Vol. 12, 1964, 163-168. 3. J.H. BRAMBLE and B.E. HUBBARD,; On a Finite Difference Analogue of an E l l i p t i c Boundary Value Problem which is Neither Diagonally Dominant nor of Non-Negative Type; J. Math, and Phys., Vol. 43, 1964, 117-132. 4. R.R. BROWN, Numerical Solutions of Boundary Value Problems Using Non-Uniform Grids; J. SIAM, Vol. 10, 1962, 475-495. 5. E.A. CODDINGTON and N. LEVINSON, Theory of Ordinary Differential Equations; McGraw-Hill, New York, Toronto and London, 1955. 6. L. COLLATZ, The Numerical Treatment of Differential Equations; Springer, Berlin, 1960. 7. V.N. FADDEEVA, Computational Methods of Linear Algebra; Dover, New York, 1958. 8. G.E. FORSYTHE and W. WASOW, Finite Difference Methods for Partial Differential Equations; Wiley, New York, London, 1960. 9. L.. FOX, The Numerical Solution of Two-Point Boundary Value Problems in Ordinary Differential Equations; Oxford University Press, 1957. 10. , A note on the numerical integration of f i r s t order d i f f e r -ential equations; Quart. J. Mech. App. Math. 7 (1954), 367-378. 11. L. FOX and A.R. MITCHELL, Boundary Value Techniques for the Numerical Solution of I n i t i a l Value Problems in Ordinary Differential Equations; Quart. J. Mech. App. Math.,-Vol. 10 (1957), 232-243. 12. P. HENRI CI, Discrete Variable Methods in Ordinary Differential Equations; Wiley, New York, London, 1961. 13. F.B. HILDEBRAND, Introduction to Numerical Analysis; McGraw-Hill, New York, Toronto and London, 1956. 14. , Advanced calculus for applications; Prentice-Hall, 1964. 118. 15. T.E. HULL and A.C.R. NEWBERRY, Integration Procedures which Minimize ' Propogated Errors; J. SIAM, Vol. 9, No. 1, 1961. 16. L.V. KANTOROVICH, Functional Analysis and Applied Mathematics; Nat. Bur. Standards Report No. 1509, translation by CD. Benster, 1952. 17. Z. KOPAL, Numerical Analysis; Wiley, New York and. London, 1955. 18.. J.D. LAMBERT and A.R. MITCHELL, On the Solution of y' = f(x,y) by a Class of Higher Accuracy Difference Formulae of Low Orders; ZAMP, Vol. 13, Fasc.3, 1962, 223-232. 19. A.N. LOWAN, The Operator Approach to Problems of Stability and Con-vergence; Scripta Mathematica, New York, 1957. 20. W.E. MILNE, Numerical Solutions of Differential Equations; Wiley, New York and London, 1960. 21. F.D. PARKER, Inverses of Vandermonde Matrices; The American Math. Monthly, Vol. 71, 1964, 410-411. 22. C.E. RIDLEY, A Numerical Method for Solving 2nd Order Linear Differ-ential Equations with Two Point Boundary Conditions; Proc. Cambridge Phil. Soc, Vol. 53, 1957, 442-447. 23. D.E. RUTHERFORD, Some Continuant Determinants Arising in.Physics and Chemistry II; Proc. Roy. Soc. Edinburgh,"A63, 1952, 232-241. 24. J. TODD, The Condition of a Certain Matrix; Proc. Cambridge Phil. Soc, Vol. 46, 1950, 116-118. 25. R.S. VARGA, Matrix Iterative Analysis; Prentice-Hall, N.J., 1962. 26. , Iterative methods for solving matrix equations; The American Math. Monthly, Vo. 72, No. 2, Part II, Computers and Computing, Feb., 1965.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical solution of boundary value problems in ordinary...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical solution of boundary value problems in ordinary differential equations Usmani, Riaz Ahmad 1967
pdf
Page Metadata
Item Metadata
Title | Numerical solution of boundary value problems in ordinary differential equations |
Creator |
Usmani, Riaz Ahmad |
Publisher | University of British Columbia |
Date Issued | 1967 |
Description | In the numerical solution of the two point "boundary value problem, [ equation omitted ] (1) the usual method is to approximate the problem by a finite difference analogue of the form [ equation omitted ] (2) with k = 2, and the truncation error T.E. = O(h⁴) or O(h⁶), where h is the step-size. Varga (1962) has obtained error bounds for the former when the problem (1) is linear and of class M . In this thesis, more accurate finite difference methods are considered. These can be obtained in essentially two different ways, either by increasing the value k in difference equations (2), or by introducing higher order derivatives. Several methods of both types have been derived. Also, it is shown how the initial value problem y' = ϕ(x,y) can be formulated as a two point boundary value problem and solved using the latter approach. Error bounds have been derived for all of these methods for linear problems of class M . In particular, more accurate bounds have been derived than those obtained by Varga (1962) and Aziz and Hubbard (1964). Some error estimates are suggested for the case where [ equation omitted ], but these are not accurate bounds, especially when [ equation omitted ] not a constant. In the case of non-linear differential equations, sufficient conditions are derived for the convergence of the solution of the system of equations (2) by a generalized Newton's method. Some numerical results are included and the observed errors compared with theoretical error bounds. |
Subject |
Boundary value problems Differential equations |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-09-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080615 |
URI | http://hdl.handle.net/2429/37402 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1967_A1 U8.pdf [ 5.8MB ]
- Metadata
- JSON: 831-1.0080615.json
- JSON-LD: 831-1.0080615-ld.json
- RDF/XML (Pretty): 831-1.0080615-rdf.xml
- RDF/JSON: 831-1.0080615-rdf.json
- Turtle: 831-1.0080615-turtle.txt
- N-Triples: 831-1.0080615-rdf-ntriples.txt
- Original Record: 831-1.0080615-source.json
- Full Text
- 831-1.0080615-fulltext.txt
- Citation
- 831-1.0080615.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080615/manifest