LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAICN 
 
 510.84 
 
 U(pr 
 
 ho. 293-300 
 aob.2 
 
CENTRAL CIRCULATION AND BOOKSTACKS 
 
 The person borrowing this material is re- 
 sponsible for its renewal or return before 
 the Latest Date stamped below. You may 
 be charged a minimum fee of $75.00 for 
 each non-returned or lost item. 
 
 Theft, mutilation, or defacement of library material! can be 
 causes for student disciplinary action. All materials owned by 
 the University of Illinois Library are the property of the State 
 of Illinois and are protected by Article 16B of Illinois Criminal 
 Law and Procedure. 
 
 TO RENEW, CALL (217) 333-8400. 
 University of Illinois Library at Urbana-Champaign 
 
 'JUN2 81999 
 ff 6 1 AH 
 
 When renewing by phone, write new due date 
 below previous due date. L162 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/computergraphict295dill 
 
rA/ Report No. 295 /T^^^TJ COO-1U69-0106 
 
 Lit? As 
 
 y.Zf> 
 
 A COMPUTER GRAPHIC TECHNIQUE FOR 
 FINDING NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS 
 
 by 
 Carl Dill 
 
 January 1969 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN • URBANA, ILLINOIS 
 
Report No. 295 
 
 A COMPUTER GRAPHIC TECHNIQUE FOR 
 FINDING NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS* 
 
 by 
 
 Carl Dill 
 
 January 1969 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, Illinois 6l801 
 
 * This work was supported in part by Contract U. S. AEC AT(ll-l)lU69 and 
 was submitted in partial fulfillment of the requirements for the degree 
 of Master of Science in Computer Science, January 1969. 
 
Ill 
 
 ACKNOWLEDGMENT 
 
 The author wishes to thank Professor C. W. Gear for his 
 valued guidance and support during the preparation of this thesis. 
 The technical assistance provided by Professor Gear is especially- 
 appreciated. Thanks are also extended to Miss Barbara Hurdle for 
 typing the final manuscript. 
 
IV 
 
 TABLE OF CONTENTS 
 
 Page 
 
 1. INTRODUCTION 1 
 
 2. METHOD OF TESTING k 
 
 3. SYSTEM DESCRIPTION lU 
 
 k. TESTS AND RESULTS 19 
 
 5. CONCLUSIONS 22 
 
 BIBLIOGRAPHY 23 
 
 APPENDIX 2h 
 
1. INTRODUCTION 
 
 Recent interest in the numerical integration of stiff 
 ordinary differential equations has led to a search for methods of 
 higher order than those currently known. This paper describes a 
 computer graphic technique which was developed to make a systematic 
 search possible. New high order methods, found using this technique, 
 will be presented. 
 
 Stiff equations are characterized by the presence of widely 
 varying time constants. Terms corresponding to very small time 
 constants decay rapidly and are of no interest because of their small 
 magnitude. Conventional methods are unstable if the step size used 
 is much greater than the smallest time constant. Stiffly stable 
 methods allow a much larger step size while maintaining stability. 
 
 Known stiffly stable methods are predictor-corrector 
 methods. The corrector may be written 
 
 (1 " 1) \ J » t Vi J «-i , '" n oVi 
 
 where a ^ 0, |a | + -g J > 0. The discretization error of the method 
 
 K. U (J 
 
 is defined as the difference between the value of y calculated from 
 
 n 
 
 (l-l) and the value of the exact solution y(t ). Define 
 
 e - y - y(t ) = y - F 
 
 n J n J n ^n n 
 
Then the errors e obey the difference equation 
 
 (1-2) 
 
 I (a. + hXB.) e _. + T(F, t , h) = 
 . - 1 i n-k+i * n' 
 i=0 
 
 where T is the truncation error for the exact solution. Convergence 
 
 requirements make it necessary to hound the solution of (1-2). This 
 
 [2] 
 
 requires the stability of the homogeneous difference equation 
 
 obtained from (1-2) with T set to zero. This equation is stable if 
 
 and only if the modulus of no root of the polynomial equation 
 
 (1-3) k 
 
 £ (a. + hAB. ) r = 
 i-0 1 
 
 exceeds 1, and roots of modulus 1 are simple. 
 
 The order p of a multistep method may be defined as the 
 
 highest degree of polynomial functions y for which the truncation 
 
 [3] 
 error T(y, t , h) vanishes identically in t and h. The truncation 
 
 error is then directly proportional to h , so p should be as large 
 
 as possible. Dahlquist has shown that if (l-3) is stable for all 
 
 A with negative real parts, then the maximum order of the method is two, 
 
 Gear, C. W. , "The Numerical Integration of Stiff Differential 
 Equations," University of Illinois, Department of Computer 
 Science, Report No. 221, p. k. 
 
 [2] 
 
 Henrici, P., Discrete Variable Methods in Ordinary Differential 
 
 Equations , John Wiley and Sons, Inc., New York, 1962, pp. 218-220. 
 
 [3] 
 
 Gear, C. W. , ojo_. cit . , p. 5. 
 
 [h] 
 
 Dahlquist, G. , "A Special Stability Problem for Linear Multistep 
 Methods," BIT 3, 1963, pp. 27-^3. 
 
Gear has relaxed these requirements and defined those necessary 
 for stiff stability. The requirements are shown in Figure 1. 
 
 hA Plane 
 
 ICURAT 
 
 4v 
 
 RELATIVELY STABLE 
 
 -0 
 
 Figure 1, 
 
 Since accuracy is usually defined in terms of order, where order is 
 defined as above, we again see that a large p is desirable. Previous 
 results show the existence of stiffly stable methods of order as high 
 as six for proper D, 0, and a. The object of the current investigation 
 was to find new methods of order greater than six. 
 
 [5] 
 
 Gear, C. W. , op_. cit . , pp. 5-7. 
 
2. METHOD OF TESTING 
 
 k k 
 
 Define p(?) = E a.^ 1 and o(c) = £ 8.5 1 . Then (1-3) 
 i=0 X i=0 x 
 
 may be rewritten 
 
 (2-1) pU) + hXa(e) = 
 
 Letting hX ■*■ °°, the roots of (2-1) approach those of 
 
 o(£) = 0. This implies that the roots of o(£) must also lie inside 
 
 the unit circle or be on the unit circle and simple, and that a(c) is 
 
 of degree k. Given any o(£), the unique method of degree k may be 
 
 \ f> 1 
 found as follows. For a method of degree k, we have 
 
 (2_2) P(0 + log(0 0(e) ~ c(C-D k+1 as^l 
 
 k , 
 
 Set £ - 1 = x. From (2-2), p(l+x) = - [ Z 6. (l+x) J ] log(l+x), 
 
 j=0 J 
 truncated to terms in x and lower. 
 
 rn 
 
 Let ■$■ = [a Q , a ±9 , a fc ] , then 
 
 L ■* Henrici, P., op_. cit . , pp. 226-228. 
 
r 
 
 1 -l 
 
 -2 3 
 
 1 -3 
 
 1 
 
 (-if 
 (-l) k_1 k 
 
 < 
 
 k 1 
 
 1 k 
 
 ( k ) 
 
 (*) k 1 
 
 + 3 
 
 k-1 
 
 k-1 1 
 ff) k-1 
 
 1 k-1 ... 1 
 
 + . . .+ 
 
 o 
 
 Q..i 
 
 o 
 
 i 
 
 -1/2 
 1/3 
 
 (-D k /k 
 
 The last matrix consists of the coefficients of the logarithm series 
 
 for (l+x). Each matrix in the inner sum of matrices corresponds to a 
 
 k 
 multiplication by (l+x) in £ 8. (l+x) . The first matrix provides 
 
 j=0 J 
 
 the conversion back to £ = l+x. 
 
We want hX values such that (2-l) has roots inside the 
 unit circle or on the unit circle and simple. This region is bounded 
 
 by the locus of - ^y^f in the hX plane for £ = e 10 , 0e[0,2tt]. If 
 
 a(c) is stable, the method will be stable at hX = infinity, so by 
 
 continuity arguments any region connected to hX = infinity will be 
 
 [7] 
 
 stable. In order to obtain values for plotting this locus, we 
 
 let y(Q) = - (fy witn £ = el " and calculate y(0) for = 0(^)tt, 
 
 where M is chosen to provide a small enough increment to obtain a 
 suitable number of points. The calculation of y(Q) is as follows 
 
 y(e) = - »*& 
 
 o(C) 
 
 kiG 
 k 
 
 10 
 
 kiG 
 
 3 k e + ... 
 
 ... + B / 9 + B 
 
 Letting R = £ a cos(jQ) , S = E a sin(jO) 
 
 j=0 J j=l J 
 
 k k 
 
 A = T. B. cos(j0) , B = I B, sin(jO) 
 
 we get 
 
 j=0 J j=l j 
 
 fn\ R + iS 
 
 Y(0) = " aTTb 
 
 (RA + SB) + i(SA - RB) 
 
 2 2 
 A + B 
 
 w 1 + (RA + SB) , _ (SA - RB) 
 Now let x = - -^ — ^— and y = - — 5 p— 
 
 A + B A + B 
 
 [7] 
 
 L,J Gear, C. W. , op_. cit . , pp. 9-10. 
 
The desired locus is obtained by plotting the curve 
 passing through the x,y values calculated for © varying over the range 
 to f, giving the lower half, and then rotating this curve 180* to 
 
 obtain the upper half, since the locus is symmetric about the real 
 
 v 
 axis. Figure 2 shows a sample plot for o(C) = K with k = 6. 
 
 -«.o 
 
 Figure 2. 
 
Given a k and the coefficients of the polynomial o(£), 
 the stiff stability of the unique method of degree k may thus he 
 investigated by visually examining the plotted locus. The 
 computer graphic approach used in the current study takes advantage 
 of this fact. A plot of the locus is created on a CRT display and 
 an interactive system is provided for manipulating k and the 
 coefficients of a(£). The great advantage is that of the speed 
 of obtaining a new plot. Conventional plotting techniques are 
 too slow to make a study of a large number of methods possible in 
 a reasonable length of time. 
 
 While speed was a major goal in designing the graphic 
 system, it should be noted that implementation of the system is 
 on a DEC PDP-8 computer with an attached Model 338 CRT display. 
 The PDP-8 is a 12 bit per word, fixed point binary machine whose 
 basic arithmetic capability is two's complement fixed point addition. 
 The necessary floating point arithmetic is thus accomplished entirely 
 by means of software. It was therefore desirable to minimize the 
 number of calculations necessary to obtain the points used in 
 plotting the locus . 
 
 The calculation of p- given in (2-3) may be optimized in 
 terms of fewest floating point operations as follows. Given k 
 and 6 , B-, » , B, generate the Pascal triangle matrix A, where 
 
 A. . = (.), or as shown in (2— U ) and where the element in the upper 
 
 left corner of A is indexed A and so on. 
 
(2-U) 
 
 A = 
 
 1111 
 
 12 3 
 
 1 3 
 
 1 
 
 ( k ) 
 
 k 
 1 
 
 To generate the matrix B representing Z 8. (l+x) , for 
 
 j = 0, l*,.. t k multiply the j-th column of A by the corresponding 6. 
 
 J 
 
 1 
 
 r "\ 
 
 i 
 
 6 o { }h< 
 
 
 
 
 v J 
 
 r ^ 
 
 i 
 
 )'" ^ ( 
 
 ( k ) 
 \ 2 j 
 
 ( k ) 
 
 k 
 
 1 
 
 
10 
 
 Sum across the rows to get the column matrix M, 
 (2-5) 
 
 M = 
 
 l 6. 
 j=0 J 
 
 k 
 
 £ J 6. 
 3-1 J 
 
 M 
 
 
 
 M 
 
 "k 
 
 Since A is an upper triangular matrix, indexing on the implemented 
 versions of the above operations is constructed in such a manner 
 that additions and multiplications are performed only on upper 
 triangular elements. In terms of (2-5), B is now given below by 
 
 B = 
 
 M, 
 
 M l M o 
 
 M 2 M x 
 
 M. 
 
 \ \-l \-2 • • • M l M 
 
11 
 
 The large savings in number of multiplications and 
 additions over straightforward calculation as in (2-3) is obvious. 
 In addition, a significant savings in storage during calculation 
 and for B itself is realized. Because of the nature of B, only 
 M need be retained. A transformation on the indexing function 
 on B may be used to obtain corresponding elements from M. 
 
 The first matrix in (2-3) is obtained by negating each 
 element of A the sum of whose subscripts is odd. Call the new 
 matrix A' . The next step in the calculation of p- is the matrix 
 multiplication 
 
 A* x B = C 
 Taking advantage of the fact that A 1 is upper triangular and B 
 lower triangular, and using M to obtain the desired elements of 
 B, the most efficient method of obtaining C is summarized in the 
 following section of ALGOL program 
 
 for i := step 1 until k do 
 
 for ] ;= step 1 until k do 
 begin sum := 0.0; 
 
 for n := i step 1 until k do 
 if j < n then 
 
 sum := sum + A(i ,n)*M(n-j ) ; 
 C(i,j) := sum; 
 end ; 
 Multiplication of the matrix C by the last matrix in (2-3) then 
 gives us - jy = - [o Q , o , , a k J . 
 
12 
 
 In order to increase the efficiency of the calculation 
 
 of y(Q), tables of values of sin(O) and cos(O) for © = 0(— )n are 
 
 stored permanently. A table look-up technique is employed to 
 obtain the values of sin(jO) and cos(jO) as needed. The technique 
 works in the following manner. Each will be of the form 
 
 ItTT 
 
 — , k = 0, 1, 2, , M. The table entries for sin and cos 
 
 may be thought of as being numbered from to M. Thus the value 
 
 of cos(— ) is found in entry number k in the cos table. To obtain 
 
 values of sin(jO) and cos(jO) we note that j0 = *•— : — — . If we 
 
 form i = ( j •k)mod(M) , then the desired value will be * entry 
 number i in the sin or cos tables respectively. The correct 
 sign may be determined by counting the number of times M must 
 be subtracted from j*k to get i. An odd number indicates that 
 
 a sign change is necessary. This is because jQ = *•— : - — falls 
 
 in a quadrant diagonally opposite to that — falls in. As an 
 
 example let M = 100, k = 30, and j = h, so that j0 = = — — . 
 
 Then i = 20 and one subtraction of 100 is necessary to get i, so 
 the desired sin or cos is minus entry 20 in the respective tables. 
 It should be obvious from the above that incrementing 
 
 over = 0(— V simply corresponds to incrementing k over k = 0(l)M. 
 
13 
 
 Values of sin(jO) and cos(j0) for use in calculating any particular 
 point on the locus are thus obtained using an i based on j and the 
 current value of k to index the tables. 
 
Ik 
 
 3. SYSTEM DESCRIPTION 
 
 The graphic system was designed to make a rapid search 
 for new stiffly stable methods possible. Easy to use interactive 
 capabilities were thus included to control the display and the basic 
 input parameters k and o(£). Implementation is on a PDP-8 with four 
 U096 word memory banks. The attached Model 338 display runs off 
 display files constructed in the computer memory. An on-line 
 teletype and a bank of 12 push-buttons are used for interactive 
 communication by the current system. 
 
 Computer storage allocation is shown in Figure 3. 
 
 Bank : 
 
 1 
 
 2 
 
 3 
 
 SYSTEM 
 PROGRAMS 
 
 LOG MATRIX 
 
 
 S's 
 
 a's 
 
 A MATRIX 
 
 SIN 
 TABLE 
 
 DISPLAY 
 FILE A 
 
 M MATRIX 
 
 COS 
 TABLE 
 
 C MATRIX 
 
 DISPLAY 
 FILE B 
 
 FLOATING 
 
 POINT 
 
 SOFTWARE 
 
 
 LOCUS 
 COORDINATES 
 
 TELETYPE 
 MESSAGES 
 
 ; .... _ 
 
 Figure 3. 
 
 The majority of the system programs in bank are concerned with 
 the calculation of -p and y{&) using the methods described in the 
 previous section. The others calculate an appropriate display file 
 for each locus. 
 
15 
 
 The display is based on a user defined display window 
 specified in terms of left (X T ) and right (X^) boundaries relative 
 
 to the real (X) axis, and upper (Y.J and lower (Y ) boundaries 
 
 relative to the imaginary (Y) axis. Only the section of the locus 
 the coordinates of which fall within these boundaries is displayed. 
 Axes are displayed as appropriate to the specified window. Each 
 display is scaled to cover the entire usable screen area. The 338 
 screen size is 102U units by 102U units, where a unit is an 
 intensifiable point, and points are numbered through 1023. A 
 coordinate x will be displayed if X < x < X , and the calculation 
 
 x - X 
 
 *102U 
 
 performs the required scaling. A similar calculation produces the 
 desired y' for Y T '- y < Y . The display consists of vectors 
 
 connecting consecutive points (x' ,y' ) on the locus. If the 
 increment used in calculating y(Q)» and thus the (x' ,y' ) coordinates, 
 is small enough, the locus will appear as a smooth curve. 
 
 In order to increase the efficiency of the display file 
 building, an appropriate switch is set whenever a point on the locus 
 first falls above, below, left, or right of the display window. When 
 each new point thereafter is tested, the switch on causes the 
 corresponding out of window condition to be tested immediately. Thus 
 only a single test, rather than a four-way test, is necessary on each 
 new point as long as the condition occurs. When the locus re-enters 
 
16 
 
 the window, a non-intensified vector is drawn from the last point 
 in the window to the new point, and the switch is turned off. 
 
 The most obvious requirement of the system is to allow a 
 user to input a k, the coefficients 8 , B, , , 6, of the 
 
 polynomial a(£), and the display window boundaries, from which a 
 
 display of the - / ~\ locus will be created. Since it is much too 
 a(fcj 
 
 cumbersome to respecify all of this data each time a specific 
 
 parameter is to be changed, interactive facilities for selectively 
 
 changing the parameters have been included in this system. The 
 
 amount of recalculation necessary to create a new display based on 
 
 the change varies significantly depending on the nature of the change. 
 
 Two display files are maintained so that the old display is left on 
 
 until the new display file, resulting from the latest parameter 
 
 change, has been built. Display control is then transferred to the 
 
 new file. 
 
 User intentions with respect to specific changes are 
 indicated by the use of the push-buttons, each of which has an implied 
 meaning. Messages asking 'for new parameters are printed on the teletype 
 as required, after which the user inputs the appropriate data. Push- 
 buttons are associated with the following actions. 
 
 1. Respecify k, S Q , 6^ , 6 k and the 3^, X^ Y y , Y L 
 
 display window boundaries. A new locus is 
 automatically displayed. 
 
17 
 
 2. Change g. . User specifies i, then new 3.. No new 
 
 locus is calculated. This allows the user to change 
 more than one 3. before getting a new display. 
 
 3. Calculate a new locus and create a new display based 
 on one or more 3. changes. 
 
 h. Specify an increment for 8. . User specifies i and 
 
 the increment amount by which 8. is to be incremented. 
 
 1 
 
 No incrementing is done here. 
 
 5. Increment 8. based on the data specified in the last 
 
 use of the button associated with action h. A new 
 locus is calculated and displayed using the new B.. 
 
 This is separated from action k so that a 3. may be 
 
 incremented without requiring new teletype input each 
 time it is incremented. 
 
 6. Respecify the display window boundaries. A new display 
 is created based on the new window. 
 
 7. Double display window size and create new display. Since 
 no teletype input is required, this enables a very fast 
 change in the display when a larger viewing area 
 
 is desired. 
 
 8. Halve display window size and create new display. 
 
 9. List the (x,y) coordinates of the points on the 
 currently displayed locus. 
 
18 
 
 10. List the a , a , a coefficients of the 
 
 U _L K 
 
 polynomial p(£) calculated in connection with the 
 currently displayed locus . 
 
 11. List the current 6^, 3 n , , 6, coefficients. This 
 
 1 k 
 
 is particularly useful after incrementing or changing 
 
 several 3. . 
 
 l 
 
19 
 k. TESTS AND RESULTS 
 
 Using the system described in the previous section, a 
 systematic search for stiffly stable methods of order greater than 
 six has been conducted. Suitable coefficients for a(£) have been 
 found which indicate the existence of new methods of order as high 
 as eight. 
 
 In addition to the stability requirement, the relationship 
 between the 3's and the number of items which must be saved to 
 provide information about the history of the solution influenced 
 the choice of polynomials o(£) tested. The corrector formula (l-l) 
 can be written more generally as 
 
 a y (m+l) = - E a y - h E 1 6 y' (m+l) 
 k y n+l k-j y n-j+l n . =Q \-j y n-j+l 
 
 By defining v as 
 
 ^ = [ V W ' y n-k + l' h y n> h K-V ' h Ck^l 1 
 
 r 8 1 
 
 and defining suitable A, &_, and F, Gear has shown that the numerical 
 predictor-corrector process may be written as 
 
 (m+l) (m) . / (m) x 
 
 y n + l = y n+l + ^ F( ^n + l > X n+1' h) 
 
 [8] 
 
 Gear, C. W. , "The Numerical Solution of Ordinary Differential 
 Equations of Various Orders," Argonne National Laboratory Report 
 ANL 7126, Argonne, Illinois, January 1966, pp. 7-9. 
 
20 
 
 For general o(£) as defined in this study, we would have k = k in 
 
 [k-l). However, if 3. = 0, i = 0, 1, , j , then k = k-j-1. 
 
 Since y contains information to be saved at step n, it is obvious 
 
 from (k-2) that small k require the saving of less information at 
 
 each step of the corrector iteration. 
 
 In order for (k-2) to make sense, we must have k >_ 1. 
 
 The least amount of solution history data must be saved when k = 1. 
 
 Looking at (4-l), this corresponds to allowing only B and B to 
 
 have non-zero values. Initial search efforts were thus directed 
 to methods involving a(£) satisfying this condition. 
 
 Since methods of order up to six were known, the most 
 immediate goal of the current investigation was to find a stiffly 
 stable method of order seven. For convenience, B was set equal 
 
 to one throughout the study, since the coefficients of o(£) may 
 always be scaled in such a manner as to make this hold. Fixing 
 B to one, and imposing the above restriction, made it possible to 
 
 study the most optimal possibilities of order seven by varying only B/-. 
 
 Stability requires the roots of o(C) to lie within the unit circle, or 
 be onthe unit circle and simple. Another restriction on a(£) is that 
 it have no common factor with p(£). Since p(^) will always have a 
 root at +1, a(C) must not have a root at +1. In terms of the above 
 restriction, with B =1, this implies that B must not equal -1. For 
 
 k K—-L 
 
 optimal methods of order seven, possible B/- are thus limited to the 
 
21 
 
 range -1 < 3,- f_ 1. By setting 3,- equal to -.99, a locus was 
 
 obtained which indicated the existence of suitable parameters 
 D and for a stiffly stable method. Thus the unique method 
 determined by the polynomial a(£) = £ - • 99£ was found to be a 
 stiffly stable method of order seven. As 3/- was incremented toward 
 
 +1, the parameter was seen to decrease toward zero until 3,- passed 
 
 -.8, at which time a suitable ceased to exist. Details of this 
 and other new methods are summarized in the Appendix. 
 
 Using this same technique, a similar study of the most 
 optimal possibilities for orders up to 15 was made. No locus was 
 found which indicated the existence of a stiffly stable method. 
 Results indicate that no parameter exists for methods of these 
 orders determined by stable a(£) satisfying the 3, 
 
 k-2 M k-3 
 
 & = condition. 
 
 By removing the restriction on non-zero 3's, less optimal 
 possibilities for methods of order greater than six were investigated. 
 Efforts were concentrated on methods of order seven and eight. For 
 these orders, each coefficient 3. i = 0, 1, , k-1 was incremented 
 
 over the range -1 < 3. £ 1 with 3, = 1 and 3. = 0, j^i,j^k. 
 
 1 k J 
 
 For k = 7, loci indicating stiffly stable methods were found for 3> 
 
 in the range .6 <_ 3> <. .8. Thus, for example, the unique 
 
 method 
 
 7 h 
 determined by the polynomial o(£) = £ + .75 is stiffly stable. No 
 
 stable methods were found when this approach was applied to methods 
 of order eight. 
 
22 
 
 5 . CONCLUSIONS 
 
 Although no simple method for theoretically determining 
 stiffly stable numerical integration methods is known, the computer 
 graphic system described here has enabled the efficient examining of 
 large numbers of possible methods. Systematic search procedures 
 have led to the discovery of stiffly stable methods of orders seven 
 and eight. 
 
 Defining optimality in terms of least amount of solution 
 history data to be saved at each step of the corrector iteration, 
 the optimal possibilities for orders seven through fifteen were 
 tested, with negative results for all orders except seven. Results 
 of the tests indicate that no optimal methods of order greater than 
 seven exist. 
 
 Removing the optimality restriction allows the unique 
 method determined by any stable a(c) to be a valid candidate for 
 satisfying stiff stability requirements. No general results can be 
 given, since the number of potential methods excludes the possibility 
 of any type of exhaustive testing. Testing of large classes of methods 
 tends to indicate the scarcity of stiffly stable methods of orders 
 greater than six. Methods of orders seven and eight were found, however, 
 while examining some of the more likely classes of possible methods. The 
 greatest advantages of the graphic system seem to be the ability to 
 rapidly accept or reject whole classes of a(£), and the ability to narrow 
 in on a method in a promising region by interactively manipulating the 
 
 coefficients of a(c) while watching the changes in the locus of - / \ . 
 
23 
 
 BIBLIOGRAPHY 
 
 1. Dahlquist, G. , "A Special Stability Problem for Linear Multistep 
 
 Methods," BIT 3, 1963. 
 
 2. Gear, C. W. , "The Numerical Integration of Stiff Differential 
 
 Equations," University of Illinois, Department of Computer 
 Science, Report No. 221, 1967. 
 
 3. Gear, C. W. , "The Numerical Solution of Ordinary Differential 
 
 Equations of Various Orders," Argonne National Laboratory 
 Report ANL 7126, Argonne, Illinois, 1966. 
 
 k. Henrici, P., Discrete Variable Methods in Ordinary Differential 
 Equations , John Wiley and Sons, Inc., New York, 1962. 
 
2k 
 
 APPENDIX 
 
25 
 
 This section -will present details of the new methods 
 of order greater than six. Two methods of order seven and one of 
 order eight will be given. 
 
 The coefficients for an optimal method of order seven are 
 
 iO 1 2 3 k 5 6 
 
 3. 
 
 i 
 
 
 
 
 
 
 
 
 
 
 
 -.99 
 
 1 
 
 a. 
 
 l 
 
 0.166 
 
 -1.361+ 
 
 U.9^2 
 
 -10.1+00 
 
 lU.lUl -13.U70 
 
 8.1+35 
 
 -2.1+51 
 
 The locus of - / A is given in Figure A-l. From Figure A-l we see 
 
 that -6.1 is an approximate upper hound on the parameter D and that 
 .5 is an approximate upper hound on the parameter 0. 
 Another method of order seven is given by 
 
 
 
 .7 
 
 a. 0.138 -1.120 3.990 -8.050 11.1+92 -10.920 7.070 -2.560 
 
 The locus of - t ' t for this method is given in Figure A-2. An 
 approximate upper bound for D is -12.1+ and for 0, .25. 
 
26 
 A method of order eight is specified by 
 
 i 1 
 
 2 
 
 3 h 
 
 5 
 
 6 
 
 7 
 
 8 
 
 S. 
 
 i 
 
 
 
 
 
 
 
 .81 
 
 -1.8 
 
 1 
 
 a. -.±62 1.U892 
 1 
 
 -6.129 
 
 1U.890 -23.763 
 
 26.587 
 
 -21.070 
 
 IO.636 
 
 -2.U78 
 
 The locus of - ~JTT ^ s gi ven i- n Figure A-3. The approximate upper 
 bound for D is -6.9 and for 0, .25. 
 
27 
 
 | I I I I | I I I I | I I 
 
 100 
 
 Figure A-l 
 
28 
 
 - 10.0/ = ; 
 
 Figure A-2. 
 
29 
 
 Rigure A- 3. 
 

 ; ^ 
 
 ? 
 
UNIVERSITY OF ILLINOIS-URBAN* 
 
 3 0112 045402051