i NUMERICAL METHODS FOR THE SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS by ARTHUR CHRISTOPHER ROLLS NEWBERY A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF ARTS i n the Department of MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September, 195&© ABSTRACT Families of three- and four-point corrector formulae are derived, which d i f fer from standard formulae i n that they express y n i n terms of more than one previously computed ordinate. It i s shown that the standard formulae are spec ia l cases of the more general formulae derived here. By theoretical argument and by numerical experiments i t i s shown that the standard formulae are often in fer ior to others which are developed i n this thes is . The three-point family, with i t s associated truncation error , i s given i n (7) and (9) of Chapter 2 on page 12. The four-point family i s given i n (41) on page 25. With the help of Rutishauser's method each family i s examined for s t a b i l i t y . In the four-point case a procedure i s described, whereby the magnitude of the coefficient i n the error term can be minimized subject to the restr ict ion that the formula sha l l remain stable e Also a theorem i s proved, which states that no stable four-point formula can have a truncation error of degree higher than f i f t h i n the step-size h„ I hereby cert i fy that th i s ahstract_is satisfactory. In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree a t the U n i v e r s i t y of B r i t i s h Columbia, I agree t h a t the 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 study. I f u r t h e r agree t h a t p e r m i s s i o n f o r ext e n s i v e copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s r e p r e s e n t a t i v e . I t i s under-stood that copying or p u b l i c a t i o n of 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 not be allowed without my w r i t t e n p e r m i s s i o n . Department of ^^A<^\^jJt^i The U n i v e r s i t y of B r i t i s h Columbia, Vancouver'^, Canada. i i i TABLE OF CONTENTS CHAPTER ONE. HISTORY OF NUMERICAL METHODS FOR SOLVING DIFFERENTIAL EQUATIONS 1 CHAPTER WO. THREE- AND FOUR-POINT CORRECTOR FORMULAE . . . . . . . . 8 1. A one-parameter family of three-point formulae 8 2s Table of numerical results obtained by various three-point formulae 14 3« A two-parameter family of four-point formulae.. . 15 4» Theorem: No stable four-point corrector formula can have a truncation error of degree higher than f i f t h i n h 19 5o Rules for minimizing the truncation error of four-point formulae while preserving the s t a b i l i t y . . . . 24 6. Graph indicating the results to be expected from various choices of parameter pairs 27 7. Table of numerical results obtained by various four-point formulae. 29 BIBLIOGRAPHY 32 ACKNOWLEDGEMENTS The author wishes to acknowledge his indebtedness to Dr. T, E. Hull for encouragement and advice received throughout the investigation which forms the topic of this thesis. Thanks are also due to the Computing Centre of the University of British Columbia for the many hours of machine time, without which the practical verification of results would necessarily have been far more limited i n scope. - 1 -CHAPTER ONE. History of Numerical Methods for Solving Dif ferent ia l Equations The earl iest numerical methods for solving ordinary d i f f erent ia l equations date back almost as far as the calculus i t s e l f . The f ie lds of astronomy and physics were among the f i r s t i n which the calculus found pract ica l application, and both these f ie lds gave r ise to ordinary d i f ferent ia l equations for which no analytic solutions could be found. Examples of such d i f ferent ia l equations are those aris ing from the computation of eccentric anomalies of planetary orbits , and those encountered in connection with simple and compound pendulums. The la t ter problems gained particular importance i n Newton's time, since maritime explorers and traders required accurate chronometers for astronavigation. The simplest and crudest method for solving f irs t -order i n i t i a l -value problems i s Euler ' s . Using the standard notation as used, for example, i n [1], the problem i s th i s : Given a d i f ferent ia l equation, y» = f(x,y) and an i n i t i a l point ( X Q , V Q ) , calculate y^r i . e . f ind the value of y when x = nh, where h i s a small quantity. Euler 's solution i s y^ = Jq +hf(xQ,yQ) where the bars denote approximations. Now suppose the d i f f erent ia l etc - 2 -equation possessed an analytic solution, y = F ( x ) , then i t would be possible to express y^ i n the form of a Taylor expansion thus: 7 L = F C X Q • h) - F C X Q ) + h F ' ^ ) + i i - F"(z) 2 ' where x^ ^ z 4 x^, and since F ' ( X Q ) i s the derivative of the function at (xQ,y 0), we may write y ± - y 0 + hf(x 0,y 0) . + | ! F " ( z ) . The value of y^ obtained by Euler's method i s what one would obtain by truncating the Taylor series after the second term. It would be exact i f F(x) were linear, because the second and higher derivatives would then be zero. The error, therefore, of Euler's predictor formula i s 0(h )• Since h i s taken to be a small quantity, attempts were made to produce formulae with truncation errors of higher degree i n h, and i t was found that various existing numerical quadrature formulae could easily be adapted for the purpose. For example the trapezoidal quadr-ture formula, ^ f(x)dx - F C ^ ) - F ^ ) - f ^ ) + tixJJ , X 0 where f(x) i s the derivative of F(x), became the trapezoidal corrector formula, - 3 -Thus, nothing more than a simple change of notation was necessary i n order to convert this quadrature formula into a formula for solving d i f f e r e n t i a l equations numerically. The above formula could be used it e r a t i v e l y to correct the values predicted by Euler's formula. This i s discussed i n [1, p«26J. The truncation error for this formula i s O(h^). In order to obtain corrector formulae with truncation errors of higher degree than third, some of the more elaborate Newton-Cotes quadrature formulae were adapted. The best known case i s that of Simpson's rule, v i z . The truncation error i s now OCtr), but this improvement was bought at a high price, because, before attempting to satisfy the above equation by iteration, i t i s necessary to know y^ and y» . Without this inform-ation Simpson's rule cannot be used. Moreover, as w i l l be explained later, this formula i s conditionally unstable. After this stage was reached i n the middle of the eighteenth century, no further progress was made u n t i l 1883 when the Adams method was published. Adams* formulae are of interest h i s t o r i c a l l y for two reasons: f i r s t l y they are the earliest of the more elaborate formulae devised specifically for the solution of dif f e r e n t i a l equationsj earl i e r After changing the notation i t reads: y2 "* 70 = 5 ^y0 + 4 y l + y P " mathematicians had been content to adapt quadrature formulae for the purpose. The second point of interest about the Adams formulae i s that they are stable. In Adams' day i t i s doubtful whether even the concept of s t a b i l i t y of formulae existed; yet i n 1952 when Rutishauser published the f i r s t comprehensive theory on the subject, he showed that the Adams formulae were not only stable but they possessed optimum s t a b i l i t y properties. Although Adams clearly recognized the fact that a good quadrature formula does not necessarily make a good corrector formula, he did not completely sever the long-established bond between the two types of formula© His corrector formulae can easily be transformed into quadrature formulae by applying i n reverse the notational changes mentioned above. However, the quadrature formulae so derived might have l i t t l e application, since they would involve evaluation of the integrand at points outside the range of integration. It was not u n t i l 1895 that a complete break with quadrature formulae was made. In that year Runge published his celebrated method [4], and since then numerous variants on his method have also been published. The relative merits of the so-called Runge-Kutta procedures and the predictor-corrector methods have been widely discussed e.g. [2, pp„ 247,248]o With the advent of d i g i t a l computers the former methods have gained i n popularity at the expense of the l a t t e r . Never-theless this thesis i s devoted to a study of predictor-corrector methods and to an investigation of means by which they may be made more competitive. The view taken i s t h i s : A highly accurate procedure must inevitably be complicated; complexity increases more rapidly with accuracy i n the case of Runge-Kutta methods than i n the case of predictor-corrector methods; the former class of methods w i l l reach 'saturation point' before the latte r ; therefore the la t t e r should not be neglected i n favour of the former. The concept of s t a b i l i t y i s of fundamental importance to the study of numerical methods for solving d i f f e r e n t i a l equations. It may be defined verbally thus: A formula i s stable i f i t i s insensitive to small errors i n the data to which i t i s applied. Since frequent reference w i l l be made to Rutishauser's paper on the subject [3], the relevant parts of his argument are reproduced here for convenience, Rutishauser's s t a b i l i t y analysis. Let the dif f e r e n t i a l equation be y' = f(x,y), and l e t y,y' be approximations to the true values of y,y', so that y = y + s and y 1 = y' + s 1 . Now y' = f(x,y) = f(x,y + s) - f(x,y) + sf y(x,y) = y' + sf y ( x , y ) . (1) Hence s ' ~ sf (x,y) . If Simpson's rule i s used, then y k + l - y k - l + f C y k + 1 + 4 y k + yk-l ]> ° r yk + i + s k + i = y k - i + s k - i + 1 C y k + i + S U + 4 ( y i + s i ) + y i - i + s i - i ] • This yields the difference equation for s: s, _ [1 - § f 1 - 4£f , 8. - [1 + § f . _ ] s. . - 0. k+1 3 y,k+l 3 y,k k 3 y,k-l k-1 If we make the simplifying assumption that f i s constant, vre obtain - 6 -the result s^ = p , where p i s the root of a quadratic equation whose roots are ,22 kf P l - 1 + h f y + | - f y + ... - e 1 It may be seen that p^ approximates to a solution of (1), for by writing kh = x and s^ = p^ we have £ s e °3 = e °i , therefore sk " P l s,' ~ f e X ^ l = f s, . k y y k The other root, p i s parasitic, and when f i s negative the magnitude of i s greater than one; hence i n this case i t causes an exponentially increasing osci l l a t i o n . Moreover this situation cannot be remedied by reducing the size of h« Next the Adams four-point corrector i s examined, y k + l - y k * feC9yk>l + 1 9 y k " H-l + yk-2 ] + 0 ( h 5 ) ' Proceeding as before, the difference equation for s i s obtained: ^ ' f f y , k + l ] s k + l + I ? * , k ] sk + t f y , k - l s k - l " 24 f7,k-2 sk-2 = °* On writing s^ = pk a s before, and taking f as constant, we obtain a cubic equation i n p with two parasitic roots. This case differs from the three-point case just discussed, i n that the parasitic roots tend to zero with h. When h = 0 the cubic equation i n p reduces to p^ - p^ = o, and therefore the parasitic roots vanish with h. Hence Rutishauser concludes that the Adams method i s stable for sufficiently small h. In the following chapter a study of three- and four-point corrector formulae i s made. In each case a new family of formulae has been developed, and experimental results are given to aid assessment of the relative merits of the new formulae and the established ones. The restriction to corrector formulae i s justi f i e d by the fact that i n practice i t i s these formulae rather than the predictors which determine the value of a method. In the theoretical part, attention i s focussed upon single first-order i n i t i a l value problems. In practice the extension to simultaneous and higher-order differential equations i s quite simple, though the s t a b i l i t y analysis i s harder. No assumptions concerning linearity of the equations are made, but i t i s assumed that the solution function has continuous derivatives of fourth and f i f t h order i n the three- and four-point cases respectively. - 8 -CHAPTER TWO. Three- and Four-Point Corrector Formulae A study of predictor-corrector procedures for solving ordinary dif f e r e n t i a l equations reveals the fact that a l l the formulae i n common use give an expression for the value of a new ordinate i n terms of one previously computed ordinate and a linear combination of previously computed derivatives. A collection of such formulae may be found i n [1, pp. 48,49]« This fact suggests the following questions: ( i ) Is i t possible to derive a generalized n-point corrector formula of which a l l known n-point corrector formulae are special cases, and i f so what would be the procedure for derivation? ( i i ) Given such a generalized formula could we carry out a s t a b i l i t y -and error analysis on i t i n such a way that this analysis would also apply to the special cases and enable us to assess their various merits? ( i i i ) Why i s i t that the established formulae only make use of one known ordinate? Would It not be possible, by using a l l the available ordinates, to produce formulae which are superior to those i n common use? The above questions w i l l be investigated for the cases n = 3 and n = 4 . Attention i s concentrated on corrector formulae only, because i t i s upon the corrector formula that the st a b i l i t y and accuracy of a given procedure depend. We w i l l start with the derivation of a generalized three-point corrector formula. Let the equation to be solved be y' = f(x,y), and - 9 -let. i t s solution- be y = :F(x). The f i r s t two rows, *6 yo y6 * i ?! y{, are given. The entries for the third row have been estimated by means of a predictor formula i n conjunction with the given d i f f e r e n t i a l equation. Denote these values by y 2 , f'2 . The generalized three-point corrector formula i s required to give a corrected value for y^ (denoted by y 2 c ) ± n terms of y Q, y^, y£ and y£ . Since the f i n a l formula may involve three ordinates, i t i s no longer adequate to adapt a quadrature formula. Instead we shall use a method of undetermined coefficients, and the formulae so derived w i l l include as special cases the results obtained by standard procedures. Let the required formula be of the form ( 1 ) y2c = a0 y0 + a l y l + h C V 6 + V i + V P ' where the a's and b's are coefficients to be determined i n such a way that when y^, y^ are substituted for y^ a n d 7^, "*"n ^ e resulting formula shall have the highest order of accuracy that i s consistent with s t a b i l i t y . The 'order' of accuracy i s defined i n the usual way thus: If (1) i s exact when y = F(x) i s any polynomial of degree ^ n, but not when F(x) i s some polynomial of degree n+1, then the accuracy of (1) i s of n t h order. In order to determine the a's and b's we express (l) i n terms of the shift operator E and the dif f e r e n t i a l operator D, so that - 10 -E F ( X Q ) - F(x Q + h) - FCx^ = y1 etc. (See [2], Chap, V). We obtain (2) E2y 0 = Ca Q + a ]_E + hD(bQ + b^E + b ^ 2 ) ] y Q . Now E => e ^ [2, p.134]. Write hD = u, E = e u, and substitute into (2). On cancelling the and proceeding formally, we obtain (3) e 2 1 1 = a Q + a 1 e u + u(b Q + b^ + b^) . Now expand each side of (3) as a power series i n u, and equate corres-ponding powers of u. We obtain an i n f i n i t e set of simultaneous linear equations i n the a's and b's. + a 1 a 1 *1 + b 0 + b l + b2 = 2 2 2 (4) 21 - + b l + 2 b2 = 21 b l 2 2 b 2 2 3 31 21 21 31 + — + a. _1 + V £ ^ 2 - ^ 41 31 31 = 41 etc. Since we only have five degrees of freedom, the five coefficients a^, a^, b0* b l * b2* W e c a n n o t bope to satisfy more than five equations of (4). The n t n equation of (4) i s derived by equating coefficients of u n ~ l i n (3)» If this equation i s not satisfied, therefore, we shall have an error term involving un**^ i . e . hn~-*,Dn""^ . In general we want the error term to have - 11 -a high degree i n h, so i f we decide to satisfy p of the equations (4) we shall always choose the f i r s t p of the equations. It w i l l be shown that i f the coefficients are chosen so as to satisfy the f i r s t five equations of (4), then the corrector formula so derived turns out to be Simpson's rule. However Rutishauser [3] points out that Simpson's rule i s unstable under certain conditions, and therefore i t may be asked whether i t i s possible, by sacrificing the f i f t h equation, to derive a stable corrector formula from the f i r s t four equations of (4). Since we now have only four constraints, but s t i l l have five degrees of freedom, an i n f i n i t e number of solutions i s possible, and these can be expressed parametrically. Using the second, third and fourth of the equations (4), we can obtain a l l the b's i n terms of a^ . The f i r s t equation of (4) then determines a.Q i n terms of a^ . Hence a^ may be regarded as a parameter which, once a value i s assigned to i t , determines the values of a , b , b and b . 0* 0* 1 2 The solutions of these equations are % = 1 - ^ , b0 = 1 2 ^ ~ 5 a l ) » b l = 12 ( 2 * a l ) » b 2 = i ( 4 + a i } • On substituting these values into ( l ) we obtain (6) y 2 = (1 - a 1 ) y Q + + ^ 7 ^ ( 4 - 5^) + 8y^(2 - s^) + y^U+a^], (5) - 12 -or (7) 72 - ( l - a 1 ) y 0 + + i L r ^ + l 6 y.. + ^ . a1(5y(J+8y|-yp] . Equation (7) represents, i n terms of a single parameter a ^ a l l possible three-point corrector formulae with truncation error of fourth or higher degree i n h. In particular, when a^ = 0, (7) reduces to Simpson's rule, and when a^ = 1 we have Adams' three-point formula. For other values of a^ we have a family of new formulae whose properties are to be investigated. Before preceding to an empirical investigation of these properties, two questions must be raised. F i r s t : What i s the truncation error associated with (7) i n terms of a^? Second: For what range of values of a^ i s formula (7) stable? Truncation Error, In deriving formula (7) we ensured that the f i r s t four equations of (4) were satisfied identically. We must therefore examine the f i f t h equation i n order to obtain the truncation error. Let the truncation error be R, and add R to the R.H.S, of (7). The f i f t h equation of (4), when written out i n f u l l becomes (S) ajtA b ^ 2%^ 2V+ 4 •'• R a | j [ 2 4 - ^ - 4bx - 4«2 3b 2] . On substituting for b^, b^ from (5) we obtain (9) R = - M 4 -a jh^DMz) 4J = 41 I f we choose a^ = 0, then the fourth order error i s zero. This 4J 41 v*6™ x0< z ^ x 2 * - 13 -i s to be expected, since (7) has now become Simpson's rule, and the truncation error for this formula i s known to be of f i f t h degree i n h. St a b i l i t y . Rutishauser's criterion for s t a b i l i t y , when applied to (7) requires that the parasitic root of the equation p 2 - a^p - (1 - a^) =0 should be less than one i n absolute value. The parasitic root i s 1 - a^. We therefore require (10) 0 < a 1< 2. When a^ = 1, as i s the case with Adams' method, the parasitic root i s zero, thus we have optimum s t a b i l i t y conditions. From (9) and (10) i t w i l l be seen that by taking a^ such that 0 < a^<£ 1, i t i s possible to produce stable formulae with less truncation error than that associated with Adams' method. Several tests were carried out i n order to establish the properties of various members of the three-point family. As might be expected, the results confirm that Simpson's rule, when i t i s stable, i s the best formula of the family. However, when Simpson's rule i s not stable, i t may be far inferior to the other formulae. This i s well exemplified by the equation y' « -2xy2, y ( l ) = | , whose solution i s y = - ~ - . * x +2 This equation was solved on the Alwac III-E computer at the University of British Columbia by each of the six three-point formulae corresponding to six equally spaced values of between 0 and 1 inclusive. At regularly spaced points on the abscissa the errors were computed and output. A step-size h = ~ was taken] 31 binary digits after the point - 14 -were carried, and X Q was taken at 1^ . From the results shown i n Table 1 the i n s t a b i l i t y of Simpson's rule i s at once apparent. The errors arising from computation by this method were i n fact alternating i n sign, but this does not show i n the table, because the spacing of the outputs shown i s an even multiple of the step-size. x a. 10 9E 3.8125 15.4375 27.0625 3^.6875 50.3125 1. 7 .8 -12 .6 -20 . 4 -24 .2 - 2 4 0 . -66 1. -13 +S - 7 .6 -7 . 4 -5 .2 -3 0 . -238 1. - 4 . 8 - 4 .6 - 2 . 4 - 4 .2 - 3 0 . -500 1. -1 .8 0 .6 0 . 4 0 .2 -3 0 . -815 1.. 6 . 8 6 .6 5 . 4 5 .2 5 0 . -1131 Table 1 showing the errors corresponding to each of six values of a^, at five equally spaced points on the abscissa. The equation i s y' = -2X3T2. The results shown i n Table 1 strongly suggest that most values of the parameter a-^ y i e l d better results than those obtained by Simpson's rule (a^ = 0). The question now arises: Which value of a^ should be chosen i n a particular case? Both theoretical investigations and some preliminary experimental results suggest that values of a^ close to zero are best, i . e . one should choose the smallest value of a^ which yields a stable formula. A joint paper by T. E. Hull and the author i s planned, i n which this subject i s treated i n more d e t a i l . The four-point case. The procedure adopted for finding the general expression (7) for a l l possible three-point corrector formulae, i n terms of a single parameter a^, may be extended to the four-point case. However, for reasons to be explained, i t i s now desirable to use two parameters. Let the four-point corrector be of the form (11) y 3 c - • a l 7 l + a ^ • h C b ^ t a ^ * V ' ^ y j ] . On rewriting (11) i n terms of the operators E and D we obtain (12) E2y 0 a (aD+a1E+a2E 2)y 0 + hDCbg+b^+b^+b^lr . Now write hD = u, E=e h D = e u i n (12), and drop the y^ . We. obtain (13) e3u = a Q + a-Leu + a ^ 2 * + uCbg+bjeU+bge^+b^e^] . After expressing each side of (13) as an i n f i n i t e power series i n u, and equating the coefficients of successive powers of u on the L.H.S. and R.H.S., we again have an i n f i n i t e set of simultaneous equations i n the - 16 -a's and b's: (14) a 0 + a l + a 2 a + 2 a + b + b + b + b 1 2 0 1 2 3 a l 2 a 2 2i + TT + \ + 2b 2 + 3b 3 1 3 f 21 h 2\ 3 3b 3 31 31 + 21 + 21 + 2 i etc. 3i The Nth row, for N> 3 i s \ 2^\ b x 2 ^ % 3 N- 2b 3 3 N - 1 (N-l ) i + (N-l ) i + (M-2)i + (N-2)J + (N-2)i " (N-l)J The f i r s t five equations of (14) enable us to express a-^ and the b's i n terms of two parameters a^ and a.^ thus: (15) a l = ! l - a 0 - a 2 ' b 0 = S ^ 9 a 0 + a 2^ > - rj^B + 19*0 - 13a2] , b 2 " " 5 a 0 " 1 3 a 2 ] • b 3 = 2 4 C 8 + a 0 + a 2 ] * When a^ and the b's are so chosen, the f i r s t five of equations (14) are identically satisfied, regardless of how a^ and a^ are chosen. The •question may now be raised: Is i t possible to choose a^, a^ i n such a way that (a) The f i r s t six equations of (14) are satisfied, and (b) The resulting corrector formula i s stable? - 17 -The sixth equation of (14) i s 5 (16) a. 2 5 a b 2^b 3 4b 3 — + £ + — + + 2. - — 51 51 4i 4i 4i 5 1 On simplifying and writing a^ and the b's i n terms of a.Q, a 2 with the help of (15)> this reduces to (17) 19a Q + l l a 2 + 8 = 0 To test for s t a b i l i t y we have to determine the upper and lower bounds of the parasitic roots of the equation (18) p3 - a 2 p 2 - a xp - - 0 . If we divide (18) by p - 1 , the remainder i s l-a^-a^a^, which by (15) i s zero. On removing the factor (p-l), from.- (I8)rwe~are l e f t with a quadratic equation whose roots are the parasitic roots of (18). This quadratic i s (19) p 2 + ( 1 - a 2)p + a Q = 0 . Eliminating a 2 between (17) and ( 1 9 ) we have ( 2 0 ) p 2 + ^ ( 1 + ao)p + aQ = 0 , ( 2 1 ) .\ p = ^ [ ^ 1 9 ( 1 + 3 Q ) + / L 9 2 ( l + a Q ) 2 - 4 X l l 2 a Q The expression under the root sign i s positive, therefore we have real roots. Denote these by p^, p 2 . From ( 2 0 ) p + p 2 = ^ ( l + a^ ,) and p ^ = a Q , .'. p 2 + p 2 = 32?(l + a^) 2 - 2 a Q . 1 2 21 u u 18 2 2 The equation Whose roots are p^, p^ i s + aj = 0. The equation whose roots are pj - 1, p 2 - 1 i s (22) (x + 1) - (x + 1) i . e . x + x - 2 ^ 2 C1 + a f / " 2 a C 2 + 2 a0 - ^ ^1 + a o ) ] + 1 + a o + 2a r t -19^ ? ^ 5 ( l + a 0 ) 2 = 0. Now Rutishauser showed that a sufficient condition for s t a b i l i t y was p^, p 2 < 1, and a necessary condition was p^, p^ ^ 1. He also showed that i n some cases the condition p^, p^ < 1 was necessary as well as sufficient. Since we do not know the nature of the equations on which our formulae may be used, we shall regard the condition p^, < 1 as being necessary and sufficient for s t a b i l i t y . It i s therefore necessary that p 2 - 1 and p 2 - 1 shall both be negative, hence the constant term of (22) must be positive. Now this constant term i s ( l - s u ) 2 ( l - ^^)» 11 Clearly this expression cannot be positive, therefore our question i s answered: It i s not possible to derive a stable four-point corrector formula satisfying the f i r s t six of equations ( 1 4 ) . Now the sixth equation of (14) was obtained by equating the coefficients of u 5 i n (13) . Since we have shown that this equation cannot be satisfied by a stable formula, and since u = hD, i t follows that a stable four-point formula must have a truncation error of degree five or less i n h. This result - 19 -may be expressed as a THEOREM. A stable four-point corrector formula cannot have a truncation error of degree higher than f i f t h i n h. It can now be seen why we chose to satisfy only five equations of (14) and to use two parameters, for i f we had satisfied six equations and used only one parameter, then this single parameter could only have generated unstable formulae. Next we wish to determine whether i t i s possible to derive a stable four-point formula with less truncation error than the Adams four-point formula, (in view of the above theorem we cannot hope to obtain a truncation error of higher degree i n h, but we may be able to reduce the numerical coefficient of h ). For this purpose we require an expression for the error term i n terms of aQ, a^. Let the error term R be added to the L.H.S. of (16), and rewrite the whole equation i n terms of a^, a 2, using (15)• Then (23) R - .19a0 + 13*2 + 8 u5 In order to investigate s t a b i l i t y we must return to (19) and consider separately the cases of real and complex roots. It w i l l be helpful to plot a graph on which each point represents a parameter-pair. We shall then be able to determine which areas on the graph contain points which generate stable formulae. The graph i s given on page 27. Case I. p 1 and p 2 are real. From (19) the condition for real roots i s (24) (1 - a 2 ) 2 > 4 a Q .. - 20 -2 2 The equation whose roots are p£ - 1, p 2 - 1 i s x 2 + x [2 + 2a Q - (1 - a 2 ) 2 ] + ( l + a Q ) 2 - ( l - a p 2 = 0. For s t a b i l i t y we require both roots of this equation to be negative, and for this i t i s necessary and sufficient that (25) (1 + a Q ) 2 > (1 - a 2 ) 2 and (26) 2(1 + a Q) > (1 - a 2 ) 2 . Any parameter-pair satisfying the inequalities (24) to (26) w i l l generate a stable formula. For the purpose of graph-plotting i t i s convenient to make the transformation: (27) x = l + aQ, y = l - a 2 . The above three inequalities then become (28) y 2 > 4 ( x - l ) , (29) x 2 > y 2 , (30) 2x > y 2 . By substituting equality for inequality signs we can find the boundary of the area which contains points satisfying the three inequalities. From the graph i t i s seen that the required area i s enclosed by the l i n e -pair (29) and the parabola (28). - 21 -Case II . and p 2 are complex conjugate. For s t a b i l i t y we require that the modulus shall be less than one. From (19) this condition reduces to 0 < aQ < 1, or after the transformation (27), (31) 1 < x < 2. From the graph i t i s seen that the required area i s enclosed by the line x = 2 and the parabola (28). By applying the transformation (27) to equation (23) we can find the locus of points which generate formulae with zero f i f t h degree error, thus 19aQ + + ® = ®> or> ^Ster transformation, (32) 19x - l l y = 0 . Since i t i s seen that this line does not pass through any area which contains, stable points, we have a graphical verification of the theorem. We can also plot the locus of points which generate formulae with the same error term as the Adams four-point formula. Adams chose a 2 = 1, = 0, therefore from (23) the error term for his formula i s ~g2. . The equation of this locus w i l l be 19aQ + l l a 2 + 8 = 19, or, after transformation, (33) 19x - l l y - 19 - 0 . Adams' choice of parameters i s represented on the graph by the point A(l,0). It can be seen from (19) that this choice makes both parasitic roots zero, therefore his formula has optimum st a b i l i t y properties. - 22 -One more question relating to four-point corrector formulae may be raised: Given that the modulus of the larger parasitic root equals c, where 0 < c < 1 , what i s the smallest possible error term subject to this restriction? Applying transformation (27) to (19) and ( 2 3 ) , we have (34) p 2 + y p + x - l = 0 , p ^ j / y 2 - ^ - ! ) _ (35) R = - 1 9 x 6 T ^ u5 . Consider the case of complex roots f i r s t . The modulus of the roots i s •fx. - 1, .*. x - 1 = (cp . This gives a straight line parallel to the y axis, and the condition for complex roots requires that we consider only that segment of the straight line which l i e s within the concave portion of the parabola (28). The question now reduces to the following: What point on the line segment l i e s closest to the line (32)? It can be seen from the graph that the required point would have to be indefinitely close to the point where the line segment intersects the parabola (28) and where y i s non-negative. Having thus located the point,., at least for practical purposes, the minimum error can be determined from (35). Now consider the case of real roots. From (34) we have max J ijTy + i / y 2 - 4(x - l ) (3.6) ly| + / y 2 - 4(x - 1) = 2c . = c, or equivalently - 2 3 -On rationalizing (36)'we obtain y 2 - 4(x - 1) » 4 c 2 - 4c|yJ + y 2 , ( 3 7 ) /. y - ± Equation ( 3 6 ) imposes the restriction y 4 2 c , so ( 3 7 ) , subject to the restriction ( 3 6 ) , represents a pair of line segments intersecting at 2 1 ( l - c , 0 ) and with gradients + — . In order to minimize the error c corresponding to the given c, we have to determine the point on the line segments which i s closest to the line ( 3 2 ) . In order to determine this point we need to know what are the extremities of the line segments. Firs t take the negative sign i n ( 3 7 ) , .'. -(x - 1) = c 2 + yc . Substituting for -(x - l ) i n ( 3 6 ) we obtain + 4 c 2 + 4 y c = 2 c . On rationalizing, this reduces to y 2 + 4 c 2 + 4 y c = 4 c 2 - 4 c|yl + y 2 , i . e . y = -|y| , so y must be non-positive, and one extremity of this line segment must be the point ( l - c , 0 ) . Since the vectorial angle of this line i s negative acute, i t can be seen graphically that the line i s directed away from ( 3 2 ) . It therefore contains no points closer to (32) than the point ( l - c , 0 ) . We have shown that the negative sign i n ( 3 7 ) i s associated with non-positive y-values. In the same way i t can be shown that the positive - 24 -sign i n (37) i s associated with non-negative y-values. In this case the vectorial angle i s positive acute and the gradient i s — -.' The c gradient of (32) i s i£ # ^ e m u s t now distinguish between the two cases In case (a), i f we start from the point ( l - c , 0) and proceed along the line segment, then we are moving away from the line (32)j therefore i n this case the point on the line segment closest to (32) i s the point (1 - c 2, 0). In case (b), as we move along the line segment we are approaching the line (32), therefore i n this case we should proceed as far as restriction (36) permits. The largest permissible y-value from (36) i s y = 2c, and for this i t i s necessary that y 2 - 4(x - 1) = 0 . Hence the coordinates of the point on the line segment nearest to (32) are given by (38) y - 2c, x = c 2 + 1 . It i s now possible with the help of (35) to determine the minimum error coefficient i n terms of c where 1 > c ^ 0 thus: (39) R = - (1 - c 2 ) u 5 , c ^ g (Case a ) R = - £JT C 19(c 2 + 1) - 22c ]u 5 , c « | i (Case b ). The values of ag and a 2 which give minimum error for a given c are: (40) f Case (a), c ^ i i , a Q = - c 2 , a 2 = 1 , 11 2 Case (b), c j§ , a.Q = c- , a 2 = 1 - 2c . - 25 -If c = ^ , then aQ and. a 2 may be determined by either of the equations (40)> since for this value of c the two expressions for R i n (39) are equal* The f i n a l form of the four-point corrector formula i n terms of two parameters i s : (41) y 3 c = a 0 ( y 0 - y i ) + y± + a 2 ( y 2 - y ^ + 2^ C8y{ + 32y 2 + 8y£ + agC^+lfrj-Sy^+yj) + a2(y0-13y.'-13y2^y3,):i-£J- [19a Q + l l a 2 + 8] h 5 F ( V ) ( Z ) . The quadrilateral OBCD contains a l l the points from which one may derive stable corrector formulae with less error than the Adams formula. (28) Parabola. Points on the convex side generate real values of parasit ic roots. (29) l ine pa i r . Points must be on the arrowed side i n order that the formula shal l be stable, (30) Parabola. For s tab i l i ty i t i s necessary that points be chosen on the concave side. However this condition i s dominated by ( 2 9 ) . (31) Straight l ine . For s tab i l i ty i t i s necessary that points be chosen to the left of this l i n e . This condition together with (29) forms a necessary and sufficient condition for s t a b i l i t y , provided h i s suff ic iently small. (32) Locus of points giving zero f i f t h degree error. (33) Locus of points giving the same error as the Adams formula. i i 1 1 - -r -1 - — _ i ... - 4 I ---- — - ... - -- .... ... ... - -- -- ... - — 1 •i -- - - • - D X V -----j t \ o Y 1 1 - - 7 i 1/ / / / \ —rf- i / J ! / > j ^ - V -1 >> 1 i -h - - -L - ) w i t i 1 / -1/ I •Jt— I -... -% — I Ay --- — ... - - - -( n y • V, . i -i 5 — X i I ! ] \ \ -• ; 1 •-• I 1 - i — - -- I ! - -- - 1 ..... ... ... - - \ - - 1 i - — -— -— -• -- — -- - - j-i { - -- 4. - -f-i i j i ^ i i i I - 28 -To i l l u s t r a t e the properties of the various members of the family of four-point corrector formulae, three second-order equations were chosen with respective solutions (a) y = sin x (b) y = sin | (c) y = sin 2x . In a l l three cases the modulus of the maximum truncation error i s bounded. In (a) the upper bound of the f i f t h derivative i s one, i n case (b) i t i s 2"^ and i n case (c) 2^. It may therefore be expected that the truncation error i s greatest i n case (c) and least i n case (b). The examples were computed on the Alwac III-E d i g i t a l computer at the University of British Columbia. A step-size h = ^ was used i n each computation. In each case the values of a^ and were automatically computed from the given c by use of formulae ( 4 0 ) . An extract of the results i s given i n Table 2, followed by an explanation of the table and some conclusions suggested by these results. - 29 -( r a d i a n s ^ c lO^E y" = -y y" - - J y» = -4y y = sin x y = sin | y = sin 2x 0, -44 419 5636 •25 -24 463 2725 5 .5 -13 486 1607 .75 -8 502 1204 s 1. -1947 0. 332 -557 -6335 .25 164 -657 -2842 10 .5 86 -70S -1597 .75 60 -738 -1174 s 1. -8481 f 0. 421 361 -1074 .25 214 535 -1056 15 .5 113 617 -809 .75 80 659 -660 1. -11628 s 0. -349 86 15114 .2$ -169 -148 8065 20 .5 -85 -252 5011 .75 -58 -298 3830 1. -13507 y 0 . -951 -557 -29928 «25 -478 -327 -15152 25 .5- -250 -230 -9166 .75 -172 -193 -6932 1. -9981 S 0. -123 792 37599 .25 -72 667 18419 30 .5 -44 618 10937 .75 -32 606 • 8217 1. -3688 Table 2. - 30 Explanation of Table 2. The three right-hand columns each indicate true-minus-computed errors which arise when solving the diff e r e n t i a l equation at the head of the column. These errors are multiplied by 10 . The column headed c indicates the largest parasitic root permitted i n the formula used. The column headed x gives the values of x i n intervals of five radians. Since each problem was worked by four or five different formulae corresponding to different values of c there are four or five entries i n each error column corresponding to each value of x. For example i n the f i r s t error column the f i r s t figure corresponding to x = 10 i s 332. The column i s headed y" = -y, and the corresponding entry i n the c column i s 0 . This means that when solving the equation y" = -y by a formula which permitted a maximum parasitic root of zero ( i . e . by Adams' formula), the true-minus-computed error was 10" X 332 at the point where x = 10 radians. Conclusions drawn from Table 2. Consider f i r s t the second error column of Table 2, giving the errors arising from computation of the equation y»» = - ^ y. This was the only equation on which an unstable formula was tr i e d . The errors arising from use of the unstable formula corresponding to c = 1 were such that a repetition of the experiment did not seem to be jus t i f i e d . .It may further be noted that this equation i s the only one of the three i n which the Adams method, corresponding to c = 0 , compares favourably with the other methods. Even here the superiority of the Adams method i s not consistent, though i t lasts for 20 radians. An explanation of the success of the Adams method applied to this equation i s not hard to find. The f i f t h derivative of the solution i s less than 2~5 , therefore the truncation error i s bound to be small. In these circumstances - 31 -i t i s clearly better to use a formula with optimum s t a b i l i t y properties, rather than one .which i s designed, at some sacrifice of st a b i l i t y , to reduce the already small truncation error. The results of the other two series of experiments display the Adams method as being consistently inferior to a l l other methods of the family, therefore i t i s hard to escape the conclusion that the Adams method can generally be bettered. The foregoing investigations may only be regarded as a f i r s t step towards the improvement of corrector formulae. Of the questions which remain unanswered the most important would appear to be: What choice of para-meters w i l l give the 'best' four-point formula for application to a given differential equation? This question and analogous questions relating to. f i v e - and six-point formulae would form a good f i e l d for further investigation. BIBLIOGRAPHY Milne, W.E. Numerical Solution of Dif ferent ia l Equations, New York - London 1953. Hildebrand, F . B . Introduction to Numerical Analysis, New York 1956. Rutishauser, H. Uber die Instabi l i tat von Methoden zur Integration gewdhnlicher Differentialgleichungen, ZAMP v o l . 3 , 1952, pp. 65-74. Runge, C. Uber die numerische Auflosung von Di f ferent ia l -gleichungen, Math. Ann. v o l . 46, 1895, pp. 167-178.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical methods for the solution of ordinary differential...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical methods for the solution of ordinary differential equations Newbery, Arthur Christopher Rolls 1958
pdf
Page Metadata
Item Metadata
Title | Numerical methods for the solution of ordinary differential equations |
Creator |
Newbery, Arthur Christopher Rolls |
Publisher | University of British Columbia |
Date Issued | 1958 |
Description | Families of three- and four-point corrector formulae are derived, which differ from standard formulae in that they express yո in terms of more than one previously computed ordinate. It is shown that the standard formulae are special cases of the more general formulae derived here. By theoretical argument and by numerical experiments it is shown that the standard formulae are often inferior to others which are developed in this thesis. The three-point family, with its associated truncation error, is given in (7) and (9) of Chapter 2 on page 12. The four-point family is given in (41) on page 25. With the help of Rutishauser's method each family is examined for stability. In the four-point case a procedure is described, whereby the magnitude of the coefficient in the error term can be minimized subject to the restriction that the formula shall remain stable. Also a theorem is proved, which states that no stable four-point formula can have a truncation error of degree higher than fifth in the step-size h. |
Subject |
Differential equations |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-01-19 |
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.0080641 |
URI | http://hdl.handle.net/2429/40168 |
Degree |
Master of Arts - MA |
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_1958_A8 N3 N8.pdf [ 2.63MB ]
- Metadata
- JSON: 831-1.0080641.json
- JSON-LD: 831-1.0080641-ld.json
- RDF/XML (Pretty): 831-1.0080641-rdf.xml
- RDF/JSON: 831-1.0080641-rdf.json
- Turtle: 831-1.0080641-turtle.txt
- N-Triples: 831-1.0080641-rdf-ntriples.txt
- Original Record: 831-1.0080641-source.json
- Full Text
- 831-1.0080641-fulltext.txt
- Citation
- 831-1.0080641.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080641/manifest