"LI B R.AR.Y OF THE UN IVERSITY Of ILLINOIS 510. 84 M 271-278 cop. 2 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN HOV ?2 ^ JAM 2 8 W2 L161 — O-10Q6 Digitized by the Internet Archive in 2013 http://archive.org/details/comparisonoftech274ratl Xi^ Report No * 2TU /)iaM t C00-1U69-0082 JUL 2 4 1S6< A COMPARISON OF TECHNIQUES FOR THE NUMERICAL INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS by K. Ratliff July 8, 1968 Report No. 27^ A COMPARISON OF TECHNIQUES FOR THE NUMERICAL INTEGRATION* OF ORDINARY DIFFERENTIAL EQUATIONS by K. Ratliff July 8, 1968 Department of Computer Science University of Illinois Urbana, Illinois 6l801 This work supported in part by Contract No. US AEC AT(ll-l)lU69 I. ;ntroduction and Description of Tests Performed The purpose of this report is to present the results of various tests made concerning the numerical integration of differential equations. The routine used in the tests is a Fortran program which integrates a set of N first order ordinary differential equations using predictor-corrector multistep methods. The theoretical basis for the methods used is given in Report No. 221, Numerical Integration of Stiff Ordinary Differential Equations , by C. W. Gear. The coefficients for these methods are given in Appendix I. The main program uses two subroutines which are related specifically to the equations being integrated. The first evaluates the function which is given in the following form: i sy r hf i (x ' y r F, = F n " y „ " h f „ l**V »y. Where h is the current step size. The second of these subroutines calculates the following matrix: I + ah 3f -1 Where a is the leading coefficient for the method being used. This matrix is used in obtaining the corrections to the predicted values of the function. The initial values of the function must be given, along with the upper bound on the error (absolute error is used), the number of corrector iterations to be allowed in one attempt to integrate, and the original step size to be used. The step size is thereafter automatically controlled. The criteria for convergence and doubling or halving of the step size are derived from the error bound provided initially. The tests involved the use of two versions of the integration program, one in which the order of the method to be used is specified, and another in which the order is automatically chosen and changed depending on the course of the integration. The order of the method is subject to change if the step size is halved or doubled. The decision to change order or not is based on the comparison of estimates of the error in various order methods . The purpose of the tests which were made was as follows: 1. The investigation of the efficiency of various numbers of corrector iterations allowed with various order methods. 2. The comparison of stiff and non-stiff methods on both (stiff and non-stiff) types of equations in terms of accuracy and time cost. 3. The comparison of various fixed order methods and automatic order selection with respect to cost and accuracy. The description of the four systems of equations used for the tests may be found in Appendix II. II. The Question of the Number of Correction Iterations It can be seen from the graph in Appendix III that it is in general not advisable to use four corrector iterations, as this is as costly or more costly than using two or three iterations. The decision whether to use two or three iterations is a rather difficult one to make. The results of this comparison made on the other equation sets are essentially the same as equation set I (shown in the graph), though the use of three iterations with the other equations seems to be more advantageous in some cases. The only marked differences between two and three iterations occurs in the high order method where three iterations i s preferable for large and medium errors and two iterations for small errors. The advantage in using two iterations in the high order-small error case however was not indicated in tests on the other equation sets. In general two iterations is more efficient for low order-small error and three iterations is better for high order-large error situations. The value of using two or three iterations seems to depend on the equations being integrated to some extent, but on the basis of these results, the remaining tests were performed using three corrector iterations exclusively. A remark is in order concerning the errors indicated on this graph. These are not the actual errors, but the error bound provided for the integration program. The actual error and desired error correspond closely for the cases where; l) low order is used with a large error bound, 2) medium order is used with a medium error bound, 3) and high order is used with a small error bound. When low order methods were used with a small error bound, the actual error exceeded the bound given. However when high order methods were used with a large error bound, more accuracy was obtained than was requested. It can be seen that the two iteration process uses generally the fewest number of function evaluations and the greatest number of matrix evaluations, while four iterations uses the greatest number of function evaluations and the fewest matrix evaluations. The use of three iterations usually falls between the other two but sometimes requires the fewest function evaluations or matrix evaluations . It was also seen that the number of corrector iterations used is not consistently related to the actual error produced. That is, three or four corrector iterations would not in general yield more accurate results than two iterations for example. III. Stiff Versus Non-Stiff Integration Methods The four graphs in Appendix IV show, for each equation set, the relationship between the actual error and the cost of the integration in terms of function evaluations required. The results are shown for stiff and non-stiff methods of order two, four, and six, and automatic order selection. Three corrector iterations were allowed in every case. In general it can be remarked that the graph of the method of order two has the greatest slope, while the graph of the sixth order method has the smallest. That is, it is preferable to use low order methods when a large error can be tolerated, and high order methods when more accuracy is desired. The graphs of the stiffly stable equations show that non-stiff methods used on stiffly stable equations are in all cases inferior to the stiff methods. The integration was not actually performed over the entire intervals of integration indicated in the description of the equation sets, since it was found to be very impractical in terms of computer time expended. The results shown are an extrapolation based on the assumption that the actual error would not have been substantially larger if the integration interval had been completely traversed. The graphs of the non-stiff equations however, show that the non-stiff methods are superior to the stiff methods of the same order, but that some stiff methods are better than the non-stiff methods of a different order. It is also evident that in most cases the differences in the two types of methods for the same order are not appreciable. It is therefore reasonable to choose stiffly stable methods over non-stiff methods when integrating an arbitarary system of equations which may be stiff or non-stiff. IV. Fixed and Varying Order Integration Techniques In making the comparison between fixed order and automatic order selection, it will be helpful to discuss the stiff and non-stiff equation sets separately, since the behavior of each group is dependent on the coefficients (stiff or non-stiff), used with the automatic order selection program. Equation sets I and II (stiffly stable) were integrated more successfully in general using the stiff automatic order selection version than when the fixed order methods were used, considering the range of error as a whole. The exception to this is that stiff fixed order methods were sometimes more effective when a relatively large error was allowed. With equation set I, the non-stiff automatic order selection version was only slightly better than the fixed order six method, and was more costly than all other methods. The behavior of equation set II is more complicated. The last remarks concerning equation set I apply to equation set II only in the small error range. For medium-sized errors, the non-stiff automatic order selection method is more effective than non-stiff orders four and six, and slightly less effective than order two. The result is a plateau- like graph. There is another sharp drop (not shown on the graph), corresponding to the transition from medium to large error. The explanation for this phenomenon seems to be that the high orders (five and six) are chosen by the order-selecting algorithm when small error is requested, while lower orders predominate when medium or large error is permitted. The integration of the non-stiff equations indicates that automatic order selection (using stiff methods) is more effective than the fixed order process using stiff methods, and approximately equivalent to or better than the non-stiff methods. The non-stiff automatic order selection program is superior to the stiff automatic version for both equation sets III and IV, especially in the large error range. However the non-stiff automatic method is in general not quite as good as the fixed order six method. 5 These results would suggest that for the integration of a general system of ordinary differential equations where it is not known whether the system is stiffly stable or not, it is more efficient to use stiffly stable integration methods. Also it is indicated that it is preferable to vary the order of the method being employed during the course of the integration, especially when using stiff methods. V. Some Additional Results Since the tests described thus far have been completed, the automatic order selection integration program has been revised. The significant changes and a comparison between the results using the old and new version are found in Appendix V. Results are shown for the new version for two different cases. The first case uses a constant ,10 , as the criterion for convergence •of the corrector. The second case uses the variable (-P-3) E = 2 * requested error where P is the order of the method being used. In the first case it can be seen that there is little correspondence between the requested and actual errors, and much more accuracy is achieved than is requested in many instances. This is true of course because the bound does not depend on the requested error. When the new version was used with the variable criterion for convergence (which is the same as that used with the old version) the correspondence between requested and actual error is essentially restored. This would seem to be a desirable feature, even though in a few cases the constant convergence criterion seems to be more economical. Overall, it can be seen that the present program is less costly then the previous version for a comparable amount of accuracy. The partial derivative matrix has, in all of these tests, been evaluated using the values of the solution at the last successful step in the integration. It has been found that a decrease in the number of function evaluations necessary can be achieved by evaluating the matrix using the predicted values of the solution. Appendix VI contains a comparison of these two alternatives. Stiff Methods APPENDIX I COEFFICIENTS FOR METHODS USED R : E ?. a a l 2 a 3 % a 5 a 6 2 -2/3 -1 -1/3 3 -6/ll -1 -6/11 -1/11 h -12/25 -1 -7/10 -1/5 -1/50 5 -60/137 -1 -225/27*+ -85/27U 1 -15/27U -I/27U 6 -180/UUl -1 -k06/khl -735/176UJ -175/176U -21/176U 1 -1/176U Non-Stiff Methods a o a i a 2 a 3 a l+ a 5 a 6 2 -5/12 -1 -1/2 1 1 1 3 ft -3/8 -1 -3/U -1/6 D h -251/720 -1 -11/12 -1/3 -1/2U R 5 -92/288 -1 -25/2U -35/72 -5/1+8 -1/120 6 -19087/60U80 -1 -137/120 -5/8 -17/96 -1/1+0 -1/720 APPENDIX II DESCRIPTION OF SYSTEMS OF EQUATIONS TESTED y^ = ((R - Ay 3 -E^ + Ey 2 ) / V y^Q) = 1 J 2 = W( yi - y 2 ) y 2 (0) = 1 y^ = Py x / C y 3 (0) = A = .0001 E = .0065 V = .0001 C = 200. P = 20. W = .0785 This system is stiffly stable and was integrated from x = to x = 500. 1 II- y-L = - Ay-L - By x y 3 y^O) = 1 y 2 = - Cy g y y 2 (0) = 1 y 3 = Ay x - By x y 3 - Cy 2 y 3 y^O) = A = .013 B = 1000. C = 2500. This set is also stiffly stable and integrated from x = to x = 50. III. y = y 2 y y 1 (0) = y 2 - - y-^ y 2 (°) = i y 3 = - .517^2 y 3 (0) = 1 This system is not stiffly stable and was integrated from x = to x = 10. IV. y-L = v 1 y x (o) = i y 2 = - y 2 y 2 (o) = i y 3 = y k y 3 (o) = o y^ = - y 3 y u (o) = l This system is not stiffly stable and was also integrated from x = to x = 10. APPENDIX III 900,000-- 810,000-- 720,000-- 630,000-- (/) Q Z o Lj 540,000-- (/) O K O 5 4 50,000- O 360,000-- O 270,000-- 180,000-- 90,000-- C0RRECT0R ITERATIONS ERROR 111 234 234 234 10 10 -6 10" 10 ERROR BD. 10 10 MAT. EVAL. 3. I R D k' E R 5' 6< 13 10 12 10 22 11 25 18 13 FN. EVAL. MAT. EVAL. U23 kkl 318 362 389 3^0 360 390 U57 U65 U8T II 6 6 17 23 11 9 21 11 51 18 26 234 234 234 FN. EVAL. 30UU 3138 31U6 1177 13U2 1367 818 890 987 787 850 1919 15^2 2193 10' 10 — r "6 lo- MAT. EVAL. 13 lU 10 18 13 10 ho 17 35 FN. EVAL 15166 15919 15923 3539 38U9 3857 1893 2000 2019 1389 1562 1632 2095 20U7 3123 2 3 4 2 3 4 2 3 4 10" 10" 6 10' e EQUATION SET I FUNCTION EVALUATION SUBROUTINE" 200 MICROSECONDS MATRIX EVALUATION SUBROUTINE-366.6 MICROSECONDS (PARTIAL DERIVATIVES) APPENDIX IV a g ro ii o UJ Q co 2 O 00 < Ll. U_ 00 i CO 5 z> o UJ < DC Ld cr o o Ld LT or O o cr Ld m u. u. »- en « ■ • e z o z u. < D O * cn ae UJ a o «z CvJ <■ (0 I- o tr. cr rr 1% UJ UJ s O Ul Q o »- -J tr (E (T 3 UJ o O o < V) SN0llVn~IVA3 NOLtONnd JO d38WnN IO Ul " Jt GO < < or H UJ co H >- LJ cr Li. o U. H 1 O UJ cc 1 cc o (J h- Ul U- ^s CO o ^ 2 ~ o o CC UJ & IMBE LOW ID o => _i UJ z < u. li- n- en 4 ■ • ® z o z u. u. <3 D o * «« ■ • e z o z Ll u. < D o * cr UJ 2z CM «• CO £2 cr UJ cr UJ cr UJ |a /-\ n Q »-_i |8 cr o cr o 3UJ SN0iivnnvA3 NOiiONnj do u38wnN APPENDIX V LIST OF SIGNIFICANT CHANGES IN THE AUTOMATIC ORDER SELECTION INTEGRATION PROGRAM Old Version New Version The step size can only be changed by a factor of 2 or .5 at any one time. The step size can be changed by a factor of 2 n (n is any positive or negative integer). This is limited only by a control involving the printing points for results. Old error coefficients are used in order-changing criteria only. New error coefficients are used for order-changing criteria and for the test for successful integration. The test for successful integration (error bound) is simply the requested error. The error bound on the integration depends on requested error, error coefficients, the order of method being used, and step size. Error bound for integration is used as follows : Error bound is used as follows Error Bound > predicted corrected y i - y i n( Error Bound) >L i=l P c y. - y. for all i in order for step to be taken. where n is the number of equations Order can be changed after step doubling, halving (non-convergence), partial derivative matrix evaluation, and it may be changed upward or downward. Order can be changed only after a certain number of steps (depending on the order), and if non-convergent, the order can only be decreased if changed. If the corrector has not converged, step halving or matrix re-evaluation is performed in the main program and order changing is allowed immediately. If the corrector does not converge, alternate matrix re-evaluation and step size changing (by a factor of 2~j) take place in the integration subroutine', until convergence is accomplished or the step size becomes too small. APPENDIX VI COMPARISON OF MATRIX EVALUATION USING LAST ACCEPTED VALUES OF THE SOLUTION AND, USING PREDICTED VALUES OF THE SOLUTION (ACTUAL ERROR/0 OF FN. EVALUATIONS) REQ. ERROR EQUATION SET I EQUATION SET II 1 EQUATION SETni EQUATION SET IV A ID" 3 p i 7.5xlO -2 315 2.1+xlO -3 Ikk UxlO~ 3 260 l.UxlO° 27h 7.1xl0" 2 308 2.7xlO _3 1*5 2.6xl0~ 3 219 l.UxlO 2 r jh A ID" 5 P 2.8xl0~ k 620 3.8x10 p 137 6.6xl0" 6 396 6.3xlO" 2 U80 1.3xlO -3 551 ■ 1.7x10 1U8 2.ixio" 5 1+01 5.3xlO -2 U80 A ID" 8 p 1 8.9xlO _T 1008 o 1.7x10" ! 360 1 2.5x10 987 3.H-X10" 5 1191 5.3xl0" T IN- 1 lxlO _T 309 Q 2.5x10" 986 l , 3-h-xIO" 5 1191 = Accepted = Predicted COMPARISON OF OLD AND NEW VERSIONS OF THE AUTOMATIC ORDER SELECTION PROGRAM IN TERMS OF: (ACTUAL ERROR/ # OF FN. EVALUATIONS) o &. BUT N o o CM m in I NO 1 O O » • Eh '-» BUOl > CQ U 1 . r-l X -3- ■a t- CM o "3 ON ON r-l O H X m o CO J- o <-\ X CM On O en o <-\ X -3/ ON r-l o -A CM m o rH •• • CC -C _, 3 OC 2 1 E 2 o ~^ * -3 NO rn CO Z W O CM e- w =0= Q CJ H CO U 1 .. . X o 2 IKK -H S cc o O =5 W O o o rH X in CM in CM CM o o rH X -a- no l/N rH CM CM 1 O •A NO CM CO 1 o rH X NO CM H O m in O co in m NO 1 o i m CO m t— CO * r. < 55 K W o o __* 3 qGw O i CM CM J- rH o s CO CO CM rH 1 O rH X ON ON 1 o 3 o l/\ J- 1 o -A CM m r- t-i -3- 1 o i m H O rH CO QKK 1 rH rH rH CM CM JKO- O H CJ CM u Z X w CO NO l/N | CO 1 CO 1 I I EW: RR.BD.O ORRECTO (-K-2), CO 1 o 3 O NO CM 1 o rH X Or CM H O •A NO NO ON O m CM m o •A in t— CO ON o r-t X m -3- CO -3 rH NO rH CM r-l M Z W O CM E- w z cc o o en • Eh <-h l/N U"N NO CO CO CO 1 o en O O o o -3" o O C_> iH rH NO rH in rH rH ON r-l O m CQ W 1 X rH X X CM X o X o X -3" .. . cc o t— r-i fr- I— -3 O t- t- in m CO 3 CC CC -H Z WKO 2KO m P-1 ■» On CM A H Eh V> Z CC W CO -3 l/N ■X, ON < CO* en 1 1 1 1 NO • Eh *—• 1 o O o o CO o t— CD O U CM o ON NO rH NO •-i CM rH CM CQ W 1 -3/ X CM X -3" X CM X rH X CM « •• • DC St! X CM m CM On NO NO f— -3/ Q CC CC 1 o\ O W O CM m ri z cc w en NO l/N t— CO CO o o * i 1 1 1 II I in • Eh — o o Lf\ o o l/N o o CM QUO] rH -3" rH t— t— t— o r-l O i-h CQ M 1 X -3 X l/N X m X NO X NO X •• ■ tc « rH CM 00 t-\ rH t— m CO in M 3 CC DC 1 l/N ON CM CM CT\ m CM -3/ ZHO CM Eh w z X in co o o 1 C— t— CO CO CO in • Eh rH O 1 IA I 1 l/N 1 1 t- Q O rH r-l rH o CM o CM o CM o ON o PQ W 1 X O r-i NO H rH H -3" rH •• • CC O -H/ CM X -3/ X CM X NO X m X ON Z 3 DC CC rH w cc o CO _3/ NO -3" r-l rH o z w u Eh u < z cc w CO l/N o o o * -3 1 NO l/N ON ao in Z> • Eh — - 1 o 1 o t— 1 1 ON Q O CM o CO —1 O o m H CO o t— o cy cq w i rH ON X LfN rH CO X H in r-l -3/ • • • cc *: X m LTN X r~i m CO X -3 X m W O CC DC 1 m m t— m rH j- JKO- rH O W O CM u z cc ta CM l/N j- C— O o o * i 1 O 1 t— o I r-l o . fH -^ o o o o 1 UA o CO 1 o Q O CM r-l in r~\ o <-\ o O t— •A o o o CQ Id 1 X rH X X CM H o •A .. . CC i>C in en CO t— CO NO X J* ON NO m a cc cc i ■a m NO -3- H c— NO CM CO r-l Z M O CM m Eh w Z DC CO NO o CO O O in 1 o -3 1 o r— rH o • Eh rH 1 -3/ o l/N 1 NO O o 1 t- 1 o Q O rH o CO rH t— O CM •A LTN o CO o o CQ Id 1 •x 1 r-t X H H -A o H .. . 0= o c— rH X r-\ if\ ON rH X m S 3 CC CC rH S cc o -3- t— t— t— CO CM CM r-{ rH o z w o NO Eh h < z cc w l/N o o o * 1 O NO 1 O ON rH O CD • Eh ^ 1 o o 1 o m 1 ao o Q O CM o MD rH l/N o ON rH CM o ON O o Of CO Id 1 rH MO X rH 1*1 X 1 H in rH .. . CC bd X m X \o ON On X r-l X l/N id q cc cc i NO CO m ON NO NO -3- rH rH H O W O CM rH aoHjav; ac CO Eh m Cvte. CO Eh rHfx. &h CO EH>-fc Cx. CO Eh rHCz.pL, CO EHrHfc.fc. CO H rH (k b. SJ,N3I0IAi30: Z O Z z o z Z O Z 1 i a 1 - aoaaa "baa O o i 6 rH rH