CENTRAL CIRCULATION BOOKSTACKS 
 
 The person charging this material is re- 
 sponsible for its renewal or its return to 
 the library from which it was borrowed 
 on or before the Latest Date stamped 
 below. You may be charged a minimum 
 fee of $75.00 for each lost book. 
 
 Theft, mutilation, and underlining of books aro reason* 
 for disciplinary action and may result In dismissal from 
 the University. 
 
 TO RENEW CALL TELEPHONE CENTER, 333-8400 
 UNIVERSITY Of ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 APR 5 1996 
 
 When renewing by phone, write new due date below 
 previous due date. L162 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/multirateordinar959orai 
 
UIUCDCS-R-79-959 
 
 UILU-ENG 79 1704 
 C00-2383-0056 
 
 A MULTIRATE ORDINARY DIFFERENTIAL EQUATION INTEGRATOR 
 
 by 
 A. Orailoglu 
 
 March 1979 
 
 IHE LIBRARY OF. THE 
 
 APR 1 7 1979 
 
UIUCDCS-R-79-959 
 
 A MULTIRATE ORDINARY DIFFERENTIAL EQUATION INTEGRATOR 
 
 by 
 A. Orailoglu 
 
 March 1979 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URB ANA- CHAMPAIGN 
 
 URBANA, ILLINOIS 61801 
 
 Supported in part by the U. S. Department of Energy, Grant US ENERGY/EY-76-S-02-2383 
 and submitted in partial fulfillment of the requirements of the Graduate College for 
 the degree of Master of Science in Computer Science. 
 

 in 
 
 ACKNOWLEDGMENT 
 
 The author wishes to express his gratitude to Prof. C.W. 
 Gear for his invaluable ideas, suggestions, criticisms and 
 supervision. Without them, this thesis could not have been 
 possible. The author would also like to thank every single member 
 of his research group and all his friends for their unfailing 
 help, suggestions and support. Special thanks are due, of 
 course, to his parents, Hristo and Atanasia Orailoglu, for 
 everything they have done, and for this the author would like to 
 acknowledge his great debt to them. 
 
I 
 
 II 
 
 III 
 
 IV 
 V 
 
 VI 
 
 iv 
 
 TABLE OF CONTENTS 
 
 INTRODUCTION , 1 
 
 IMPLEMENTATION 3 
 
 RESULTS 11 
 
 a - Efficiency Comparisons 11 
 
 b - Parametric Study 20 
 
 PROCRAM USACE 25 
 
 FURTHER IMPROVEMENTS 2 7 
 
 CONCLUSIONS 29 
 
 RIRLIOCRAPHY 30 
 
 APPENDIX 31 
 
 Program 31 
 
I. 
 
 INTRODUCTION 
 
 This thesis will concern itself with the development of a 
 p a c 1- a g e for solving systems of ordinary differential equations 
 using multlrate techniques. In a paper hy C.W. Cear [1] some 
 of the theoretical questions associated with multlrate methods 
 are discussed. Yet most of the work that has he en done up to now 
 in solvinp, systems of ordinary differential equations has 
 concerned itself with solvinp, the whole system simultaneously. A 
 disadvantage of solving systems of ode's simultaneously is clear 
 even intuitively. When some dependent variables change much more 
 rapidly than others, t he whole system will be forced to take time 
 steps to suit the fastest dependent variable. Coupled with the 
 fact that it is easy to construct examples in which different 
 equations will be slowly moving at different parts of the 
 integration interval, it becomes obvious that to integrate the 
 equations of the system of ode's separately would be 
 advantageous. The idea of integrating the equations of the system 
 of the ode's separately is thus the one that gives rise to 
 multlrate methods. 
 
 Attractive and easy though it sounds, multlrate methods 
 are not simple to implement. l/hen we are solving systems of 
 ordinary differential equations, we usually have systems in which 
 the equations are coupled . Therefore we cannot let each equation 
 proceed as far as its local error estimator would permit. The 
 reason for this becomes clear when we consider that in the system 
 
y ' = f ( y ) t h o behaviour of y depends not only on y but on the 
 whole vector y. Therefore when we are advancing alon;i the 
 solution curve nf y , our results, which depend on our estimates 
 of v, will only he as pood a s our estimates of y. Our estimate o f 
 
 y, will depend , on the other hand, on the solution curve of y 
 and on the length of the step taken. 
 
 Pec a use of the peculiarities of the mult irate method, we 
 cannot he content with the recommended step that our ode 
 integrator rot urns. V'e have to use more careful monitoring and 
 checking of acceptable next steps. An experimental package, 
 MR ATP , which uses mult irate methods, has been developed based on 
 DIFSUB [2] . In this package, a number of different schemes have 
 been tested. Descriptions of implementation techniques, different 
 sell ernes tested and experimental results follow. 
 
II. IMPLEMENTATION 
 
 The general organization of the scheme is to view the 
 
 system as n single equations, each being integrated separately. 
 
 At any intermediate stage , each equation has been integrated to 
 
 a different time called t . This is a vector with a component 
 
 cur r r 
 
 for each member of the system. A step size h has also been 
 
 estimated for the next step for each equation, leading to the 
 
 definition oft » t + h . 
 
 poss cur r 
 
 The next equation to be integrated one step is that with 
 
 the smallest t 
 
 poss 
 
 During the integration of a particular 
 
 equation, for example y , the function f,(y,;y.»j^i;t) must be 
 calculated. The values of y.(t) are calculated by extrapolation 
 from the last step. 
 
 Were it not for error control of stepsizes, the task 
 would be straight-forward. The difficulties arise because of the 
 possibility that the estimated step for an equation may be far 
 too large, so that that equation is not integrated but its 
 correct values are extrapolated too far into the future for use 
 in other equations. If this is not discovered in time, it is 
 virtually impossible to back up the equations which have been 
 integrated with the incorrect values. 
 
 The need to monitor time steps suggested by the ode 
 integrator for a single equation of the system of ordinary 
 differential equations has prompted the following additions to 
 
h 
 
 - t + h) 
 poss cur r 
 
 the basic organization. 
 
 1 - Since after acceptable projected times (t 
 
 have been calculated, other equations may be integrated, an 
 estimate of the effect of the changes on the acceptability of 
 the projected time steps has been included. If the effects 
 pass a certain acceptability level, the projected time steps 
 are recalculated. 
 
 2 - To determine the influence of each integration on the 
 
 projected time steps of the coupled equations, information on 
 the coupling of the systems is needed. The Jacobian of f is 
 needed since our estimate of the effect depends on it. Also, 
 to economize on the work done by the interpolating routine, 
 information on which system is being influenced by which 
 variable is needed. The reason for this is that usually only 
 a small number of variable values are needed in order to 
 evaluate the equation y ■ f.(y;t). All this information is 
 presently being requested of the user, but in a more 
 comprehensive software program, this information can probably 
 be deduced by the program itself. 
 
 With this information the following implementation was 
 
 chosen : 
 
 1 - As an ode integrator, DIFSUB was used. The integrator was 
 chosen for its ready accessibility and also for its particular 
 ease in extrapolating to different time values. Only the Adams 
 predictor-corrector method (of the three possible options of 
 the program) was used. 
 
Since DIFSUB does not use any multirate integration, 
 
 but integrates the system simultaneously, extensive 
 
 modifications occured. At present, the multirate integration 
 
 package, MRATE, handles only one equation per call instead of 
 
 the entire system. To facilitate correctness and ease of 
 
 debugging, the ode integrator was broken down into the 
 
 following subroutines, as it was being modified. 
 
 Scale (to scale information associated with y depending 
 on h ) 
 
 Predc t (Predictor) 
 
 Cor ec t (Corrector) 
 
 (Step and order change mechanism) 
 
 Soe 
 Sharp 
 
 (to reduce to a first order method if integration 
 does not seem to converge) 
 
 Update (to update variables associated with order each time 
 order used is changed) 
 
 E r r c h k (to check if converged corrector's error estimate is 
 within the tolerance specified) 
 
 Also MRATE now calls the interpolating routine in ter p 
 before a call to di f un 1 (which evaluates f (y;t) ) to supply 
 the values of y through interpolation. 
 2 - Since DIFSUB saves all known derivatives of each variable at 
 its current time, the natural choice for interpolation was 
 Taylor's expansion. 
 
 DIFSUB saves the j-th derivative of y scaled by the 
 quantity h-'/j! where h is the next step recommended. Since a 
 Taylor's expansion would be 
 
 h* - t t - t 
 
* - (h ) y (t ) 
 y (t ♦ h ) - Z J 
 
 J J m-0 
 
 k (h*) m yj rn) (t ) 
 = E 
 
 m=0 
 
 m! 
 
 it follows that all we need in order to interpolate to the t. 
 
 point is rescaling the information saved in DIFSUB. If s is 
 r m 
 
 the saved information for the m-th derivative in DIFSUB, our 
 interpolating procedure produces: 
 
 yj( tj + h ) 
 
 k 
 
 m= 
 
 m 
 
 Within this structure, it was easy to test several 
 schemes for choosing which variable to integrate one step and the 
 size of that step. 
 
 A simple scheme of implementing a multirate method was 
 tried first. The variable which would integrate to the smallest 
 point in the time axis was integrated first based on the 
 recommended time steps returned by the ode integrator. This 
 scheme helped keep all the equations moving within at most the 
 largest integration stepsize that can be taken (hmax). However, 
 this scheme proved inadequate since the effect of one rapidly- 
 changing variable on the projected integration step of another 
 equation was not taken into account. Therefore, integrations 
 frequently tended to cut down the recommended integration step, 
 thus upsetting the flow of the program. Equations with high 
 projected times did not cut down on the integration step until 
 
the actual integration of the equation, by which time all the 
 related variables had moved too far along the tine axis, thus 
 making interpolated results unreliable and inaccurate. The 
 result of this was that it was frequently impossible to integrate 
 the equation within the requested tolerance. 
 
 In a second scheme, the acceptability of the recommended 
 steps returned by the ode integrator was checked. This was 
 accomplished by actually calling the ode integrator once more at 
 the end of each integration with the recon mended stepsize. The 
 step size of the second call ( which could be equal to or less 
 than the recommended stepsize) was used as a possible integration 
 stepsize. This pseudo-call to the ode integrator changed none of 
 the values in the program. Also, to monitor the acceptability of 
 the projected integration step thus obtained, an estimator of the 
 changes introduced as the other equations were integrated was 
 added. As the changes that were introduced exceeded a certain 
 acceptability level, the equation whose recommended stepsize thus 
 became no longer reliable was integrated to bring it to the same 
 time as that of the last equation integrated. Calling the ode 
 integrator twice for each integration step ( once for determining 
 the recommended next step and once for actually integrating) was 
 quite expensive. Therefore, another parallel scheme was tried, 
 with basically the same structure. It used the recommended step 
 returned by the ode Integrator without rec becking It as a 
 possible integration step but also monitored the effect of 
 changes Introduced by other variables. 
 
P. oth of the above schemes producer! acceptable results. 
 
 Yet one of the c onp 1 i ca t ions that arose with the schenes above 
 
 was that they were forced to take snail steps. The reason was 
 
 that vhen one equation was forced to actual integration because 
 
 of the unreliability of its projected integration step, it was 
 
 integrated to the same time as the last equation integrated. 
 
 This ha (I the disadvantage of occasionally causing the forced 
 
 integration to t a V. c a very small integration step when the 
 
 difference between the time of the last integrated equation and 
 
 the equation whose recommended step had become unreliable was 
 
 snail. Thus the integrator was spending a large amount of tine 
 
 over a small interval when it was not necessary. Furthermore, 
 
 once having reduced the steps ize greatly DIFSUB can only increase 
 
 it slowly without introducing a lot of error. At the same time, 
 
 using a smaller steps! 7. e would not have helped the problem and 
 
 using a larger one might have conflicted with the scheme of 
 
 integrating first the equation with the smallest 't 
 
 poss 
 
 The problem with the schemes above was that of requiring 
 an actual integration step to be taken every time the projected 
 integration step was felt to be unreliable. A change in the 
 progran logic was therefore undertaken and now every time that 
 the projected stcpsiz e was felt to be unreliable it was only 
 r ec h e c 1: e d . Thus the final program logic was: 
 
1 - at the end of each integration step determine a projected 
 
 possible integration step. ( This can be done by doing a 
 pseudo-integration or by simply accepting the step returned by 
 the ode integrator as the possible integration step.) 
 
 2 - every time one of the projected integration steps is felt to 
 
 be unreliable do a pseudo-integration to determine a new 
 reliable projected integration step. 
 
 One question that has not been discussed up to now is how 
 
 the program checks that a projected integration step has become 
 
 unreliable. This is accomplished by having a mechanism for 
 
 accumulating the worst possible case errors that might be 
 
 introduced by a connected equation. The analysis is that we want 
 
 to find the errors introduced to the j-th equation as a result of 
 
 the i-th variable being integrated. Assuming that the i-th 
 
 variable predicted value for time (t ), was y. before 
 
 currj ' i 
 
 integration and is y after integration we want to estimate the 
 
 error introduced to the j-th equation as a function of the 
 
 difference between the Interpolated value to the (t ) . of y. 
 
 r poss j 'i 
 
 based on its value at the old (t ), and its interpolated value 
 
 curr i v 
 
 based on the value at the new (t ),. Let y. stand for the 
 
 curr 1 1 
 
 value interpolated from the old (t ). and y. stand for the 
 
 curr 1 y i 
 
 value interpolated from the new (t ).. Then we have for our 
 
 cur r 1 
 
 local error estimate: 
 
 h ji f j (y j ;y i } • W^i 
 
 9 ii 
 
 j ay. (y i" y i ) 
 
10 
 
 Now, (y - y") is our local error estimate, had we actually 
 
 inteerated to (t ). point. We extrapolate this difference fron 
 b poss j r 
 
 the predictor-corrector difference ( which is the estimate for 
 local error) at the point ^curr^i* Let R 8tand for 
 
 (t ) . - (t ), 
 poss ,1 cur r i 
 
 . Then the extrapolation produces: 
 
 (y< - y, ) - el R J 1 
 1 j-0 
 
 J.i 
 
 and therefore our local error estimate becomes 
 
 3f k 
 
 h,T— e, I R J 1. . 
 
 which is 
 
 6 J + 
 
 q.i 
 
 " * 
 
 id 
 
 . t ctj 
 
 3 ii 
 
 j 9 Yj 
 
 i refers to the variable that is actually being integrated 
 
 j refers to the variable connected to the equation actually 
 
 being integrated 
 
 k is the order of the method being used to solve the i-th 
 equation 
 
 h refers to the actual integration step 
 
 e refers to the estimate of the local error 
 
 1 refers to the coefficients of the method being used. 
 
11 
 
 III - RESULTS 
 
 a - Efficiency Comparisons 
 
 The two schemes that were tested in MRATE as multirate 
 me t hods were : 
 
 1 - No p r e in t eg r a t i on 
 
 In this scheme, an actual integration step would be 
 
 followed by the setting of projected time to current time 
 
 plus recommended time step returned by the ode integrator. 
 
 A pseudo-integration step would only be taken In case the 
 
 error control mechanism indicated that a certain t had 
 
 poss 
 
 become unreliable. 
 
 2 - P relnteg ration 
 
 The actual integration step in this scheme would be 
 
 followed by a pseudo-integration step which would test the 
 
 recommended step returned by the actual integration. t 
 
 poss 
 
 would be set to time current plus step used in the pseudo- 
 integration. Again a pseudo-integration step would be taken 
 
 in case our error control mechanism indicated that t 
 
 poss 
 
 for a certain equation had become unreliable. 
 
 A number of different test problems were run on both 
 of the schemer. and a sizable amount of statistics was 
 generated for each integration. The statistics gave us 
 information on how many function evaluations were needed, how 
 many real or pseudo-calls to the ode integrator were made, 
 how many times corrector was called, how many of those times 
 
 < 
 
12 
 
 it could not converge and therefore had to cut down the 
 stepslze, how many times the converged corrector's local error 
 estimate failed to meet the user-specified local error 
 tolerance, how many times method had to readjust the projected 
 integration step, and also order change data, mean stepslze 
 and standard deviation on stepslze. Out of this abundance of 
 data, the number of function evaluations, defined as how many 
 times the expression y,(t) ■ f.(y 5 t) was evaluated, was 
 chosen as the criterion to use in evaluating the efficiency 
 performance of the method. The decision was made because the 
 number of function evaluations seemed to be connected to 
 execution times and thus it offered a way of combining all the 
 rest of the da ta . 
 
 The second scheme was expected to be more expensive 
 then the first in terms of accuracy per function evaluation 
 ratio, since the step recommended by our ode integrator should 
 be a reasonable one unless the solution curve of the equation 
 is changing very rapidly. We included the 9cheme nonetheless, 
 to see if a better estimate of the next stepslze was needed. 
 As is to be seen from the results that follow, the 
 pre in tegra t ion scheme did on the average 20 - 30 % worse in 
 terns of accuracy per function evaluation than the original 
 non-nu 1 1 i ra te DIFSUB. The no-p re in t egra t ion scheme, on the 
 other hand, consistently did two to three times better than 
 DIFSUB and three to four times better than the pre i n t egra t ion 
 scheme. For most problems, it seems that the stepslze 
 
13 
 
 recommended by our integrator is accurate enough and does not 
 Justify the cost of doing additional integration steps to test 
 the recommended stepsize. 
 
 The test problems that were tested were as follows: 
 
 1 - The system of ode's derived from the 
 
 3 u u 
 
 3x ' 8t 
 
 partial differential equation by forward differencing. This 
 was tested with u(0,t) ■ and different x-axis initial 
 conditions such as: 
 
 u(x,0) - j 
 
 u (x , ) ■ sin(x) 
 
 u(x.O) - -Uf^- 
 This produced the following scheme : 
 
 y x (t) - 
 
 y 2 (t) - 
 
 (y 2 (t)- yi (t)) 
 
 dx 
 
 (y 3 (t)-y 2 (t)) 
 dx 
 
 y 9 (t) 
 
 (y 10 (t)-y 9 (t)) 
 
 dx 
 
 'm - (y " (t '- y io (t)) 
 
 i u aic 
 
 yio (t) 
 
 dx 
 
 Systems with 10 and 20 ode's were tested. 
 
 2 - The system of ode's derived from the same partial 
 differential equation with the same initial conditions. 
 
14 
 
 This time centered differencing instead of forward 
 differencing was used except for the first equation. 
 This produced the following equations : 
 
 y L (t) - 
 y 2 (t) - 
 
 (y 2 (t)-y l <t)) 
 
 dx 
 
 (y 3 (t)-y 1 (t)) 
 
 Fdx 
 
 y 9 (t) 
 
 yio (t) 
 
 (y 10 (t)-y 8 (t) ) 
 
 2dx 
 
 y 9 (t) 
 
 2dx 
 
 Again systems with 10 and 20 equations were tested 
 
 3 - The system of ode's made up of n of the followinj 
 equations : 
 
 My, + l)(x ,-t) 
 - [e X " ( yi + l) + e 
 
 * (y 1 - 1 )(x 1 _ 1 +x ) 
 
 (y^DJ 
 
 - i j ^i'y^ 
 
 (x 0" x n ) 
 
 M ». 5 
 
 k«=20. 
 T-. 5 
 
 and initial conditions of (-1,1,1,-1,-1,1) for the system 
 
 T 
 of six equations and (-1,1,1,-1,-1,1,-1,1,1,-1,-1,1) for 
 
 the system of 12 equations. 
 
15 
 
 Results for all of these systems were obtained and 
 plotted. Plotting was done in terras of natural log of 
 I Iglobal error| | vs. natural log of function evaluations. The 
 best fitting line In the least-squares sense was passed 
 through all the schemes being compared, namely, no- 
 pre In tegra t Ion and pre in tegra t ion schemes of MRATE, and 
 DIFSUB. The L_ error was found by obtaining analytical 
 solutions where possible, and otherwise by getting a true 
 answer estimate from the values of DIFSUB employing a strict 
 local error criterion. Graphs uniformly show superiority of 
 the no-preintegration scheme vs. the non-multlrate method and 
 possible limited merit for the pre in tegra t ion scheme vs. the 
 non-nu 1 t i r a te method at the high accuracy region of the graph. 
 Some of the graphs are presented on the following pages. (In 
 the graphs, MRATEl refers to the no pr Integra t ion scheme and 
 MRATE2 refers to the pre in tegra t ion scheme of MRATE.) 
 
r gd 
 
 / 
 
 */ 
 
 H / 
 
 111 
 
 to 
 
 in 
 m 
 
 CD 
 
 in 
 
 in 
 
 03 
 
 <D 
 
 in 
 
 * 
 
 in 
 
 
 Q 
 
 in 
 
 o 
 
 in 
 
 in 
 
 rH 
 
 rl 
 
 OJ 
 
 OJ 
 
 I 
 
 1 
 
 1 
 
 I 
 
 i 
 
 JZ DU. _jcu luactDo: 
 
TO 
 
 in 
 
 in 
 
 CD 
 
 CD 
 
 • 
 
 W 
 
 / 
 
 
 in 
 
 (0 
 
 J 
 < 
 > 
 ui 
 
 z 
 a 
 H 
 
 m o 
 . z 
 
 • £ 
 
 IL 
 
 a 
 
 z 
 J 
 
 in 
 I 
 
 o 
 
 H 
 
 I 
 
 in 
 
 I 
 
 o 
 
 OJ 
 
 I 
 
 in 
 I 
 
 <D 
 
 jz an -joj ujccceoa: 
 
in 
 
 o 
 H 
 
 H 
 111 
 
 f / 
 
 / 
 
 o 
 
 H 
 
 
 OP 
 
 in 
 m 
 
 GD 
 
 in 
 
 in 
 
 < 
 > 
 
 UJ 
 
 z 
 a 
 
 H 
 O 
 
 • z 
 
 IL 
 U. 
 
 a 
 
 z 
 
 in 
 i 
 
 in 
 
 I 
 
 in 
 
 i 
 
 o 
 
 rl 
 
 I 
 
 in 
 
 fti 
 
 H 
 
 I 
 
 in 
 
 H 
 
 I 
 
 in 
 
 * 
 
 H 
 
 i 
 
 JZ DLL JA1 LUQlKOa: 
 
in 
 
 CD 
 
 / 
 
 / 
 
 / 
 
 rl 
 Li 
 
 / 
 
 / 
 
 CD 
 
 in 
 
 • 
 
 in 
 
 CD 
 
 CO 
 
 -J 
 < 
 > 
 LI 
 
 Z 
 
 a 
 H 
 
 H 
 U 
 
 z 
 
 L. 
 L. 
 
 a 
 
 z 
 -J 
 
 in 
 
 cu 
 
 in 
 l 
 
 in 
 
 i 
 
 o 
 
 rl 
 
 I 
 
 in 
 
 w 
 
 H 
 I 
 
 in 
 
 H 
 
 l 
 
 m 
 
 rl 
 
 I 
 
 CD 
 
 _IZ OU_ _J(U ucuc^occ 
 
20 
 b - Parametric Study 
 
 The error in a single integration step is assumed to 
 arise from two sources : 
 
 1 - Errors in the integrator applied to a particular component 
 
 when the values of all other components are calculated by 
 extrapolat ion , 
 
 2 - Errors introduced when coupled components are Integrated. 
 
 One of the ad hoc parameters in the above was the 
 ratio of errors "allowed" in each of those components. The 
 original integrator recommended a stepsize which would have 
 kept the local error at the next step just within the local 
 error tolerance (with a certain safety margin), if the 
 equation's solution curve had not changed radically. It was 
 obvious that recommending as big a step as this in the 
 multirate method would have been unwise, since that meant 
 that our error control mechanism would report that a step had 
 become unreliable any time a coupled equation had been 
 integrated. Therefore, the program was modified to return a 
 recommended stepsize which would produce an error for the next 
 step within [ 0.7 * local error tolerance ] . Thus our error 
 control mechanism did not declare a stepsize unreliable until 
 the local error estimate surpassed [ 0.3 * local error 
 tolerance J . The results thus obtained were quite successful, 
 as presented above, yet there was still room for possible 
 improvement by a more careful testing of the above factors 
 
21 
 
 which were implemented in an ad hoc manner. Therefore, a 
 parametric study was done on the above factors to find a 
 near-optimum set of these factors. Problem (a) of the 
 preceding section was used with different initial conditions. 
 Unreliability thresholds were set from [ 0.1 * local error 
 tolerance] to [ 0.9 * local error tolerance ] with 
 corresponding changes made to integrator-returned recommended 
 stepsizes . 
 
-~ CD 
 
 -urn 
 
 .. r» 
 
 _.o 
 
 .. m 
 
 -.<*• 
 
 CO 
 <L 
 UJ 
 
 » Q 
 
 ' -J 
 O 
 
 X 
 W 
 
 III 
 
 K 
 
 CU I 
 
 H 
 
 >• 
 
 H 
 
 _l 
 
 i- 
 o. 
 
 UJ 
 
 o 
 
 o 
 
 , • < 
 
 o 
 
 H 
 
 in 
 
 CO 
 
 in 
 
 CD 
 
 m 
 
 in 
 
 -JZ QU. U.DZOHHQZ UJX-ICD 
 
-r- CD 
 
 J- CD 
 
 .- i> 
 
 .-© 
 
 -- in 
 
 .-*■ 
 
 (0 
 
 UJ 
 
 .. CO 
 
 o 
 
 H 
 
 in 
 
 CD 
 
 CO 
 
 in 
 
 
 in 
 
 • 
 
 • 
 
 • 
 
 m 
 
 CD 
 
 r* 
 
 -IZ QU_ L.DZUI-HOZ UJX-JCO 
 
 .-CU 
 
 -- ri 
 
 in 
 
 Q 
 
 o 
 
 X 
 (0 
 UJ 
 
 tr 
 
 i 
 
 > 
 I- 
 H 
 -J 
 H 
 
 m 
 
 < 
 
 Q. 
 Ul 
 U 
 
 o 
 
 < 
 
24 
 
 The preceding graphs were plotted by passing the 
 least-squares fit for each parameter to the equation 
 
 In error - a (In function evals) + 
 and then Interpolating the results for three different 
 values of global error ( actual values used were l.E-3, 
 i.E-6, l.E-9 ) . Then for each different global error, 
 points were plotted and the least-squares fit for each 
 
 global error was passed through the quadratic equation: 
 
 2 
 function evals - a (threshold pararo) + (.1 (threshold param) + y 
 
 (for constant global error ) 
 
 and threshold parameters which minimized number of function 
 
 evaluations were obtained. These were: 
 
 problem ( a ) problem ( b ) 
 
 global error ■ l.E-3 
 global error - l.E-6 
 global error ■ l.E-9 
 
 . 4414 
 .3624 
 .3439 
 
 .4075 
 .3833 
 . 3656 
 
 The high accuracy results are probably more 
 significant since results in the low accuracy region are 
 constrained by the maximum stepsize parameter. As can be 
 seen from this study, then, the best breakdown seems to be 
 around (0.60 - 0.65 ) * eps for returned stepsize by 
 integrator and (0.35 - 0.40 ) * eps for acceptability 
 threshold. 
 
25 
 IV - Prop ran Usage 
 
 The implemented program, MRATE, can be ti9ed for any 
 system of ordinary differential equations. It is most useful 
 nonetheless, for systems of coupled ordinary differential 
 equations which have differing speeds of integration. As such, it 
 can be particularly suited to a variety of systems of coupled 
 ode's obtained from discretizations of partial differential 
 equations. 
 
 The package as it presently stands needs three pieces of 
 information from the user. 
 
 1 - A main program to supply minimum stepsize, maximum stepsize, 
 
 local error tolerance, initial time, final time, report mode, 
 a y array containing initial values, four integer arrays to 
 implement two lists of which variable influences which 
 equation and vice versa, the dimension parameters for the 
 preceding arrays, number of equations, a double precision and 
 an integer work array. The main program need only call the 
 package once. 
 
 2 - A function evaluator d i f un 1 . The arguments to this 
 subroutine are : 
 
 a - The current values of the variables 
 
 b - The derivative ( output parameter ) 
 
 c - The index of the equation being integrated 
 
 d - The number of equations 
 
26 
 
 The user need only set the derivative. 
 3 - A double precision function ,1 a cob to evaluate Jacobian. 
 Arguments are : 
 a - Partial with respect to this variable 
 b - Equation number whose partial derivative we are 
 
 calculating 
 c - Current values of all variables 
 d - Number of equations 
 
V - Further Improvements 
 
 There appear to be two areas of fur titer research and 
 possible improvement for the multirate p a c k a p e developed. 
 
 As far as the software is concerned, numerical evaluation 
 of the Jacohian should eliminate the need for the tedious 
 evaluation and coding of the exact Jacobian. This would be 
 useful, since the symbolic Jacobian evaluation by hand is 
 particularly error-prone and tedious for real problems. 
 
 Of course another possibility would be implementing this 
 pad. age as a part or an option of a more general differential 
 equation integrator such as the MODEL [3] package , currently 
 under research at the University of Illinois. besides 
 alleviating the necessity of supplying a Jacobian evaluation 
 routine, this would eliminate the need to supply anything but a 
 specification of the problem in a more natural language. If this 
 gets implemented, it will also be beneficial to have code for 
 discrctiz ing partial differential equations by the method of 
 lines available in the pad age [4] . The multirate method can be 
 especially suited to such systems. 
 
 In terms of numerical analysis, there seems to be a lot 
 of room for improvement on three different issues which warrant 
 further research ■ 
 
 a - r. valuation of t he J acobian . It seems it would be sufficient 
 to use an approximate Jacobian, until the program finds that 
 
28 
 
 the approximate Jacobian is not useful any longer (i.e. when 
 ace u m u late d errors do not give a good estimate of the 
 reliahility of the next integration step) . It should be 
 sufficient to update Jacobian only then. 
 
 b - E rror e stimates f or c oupled variables when an e quation 1 s 
 he in p integrated . Our present scheme relies on a worst 
 possible case analysis and assumes that the errors arc the 
 larpcst possible. Since our control mechanism is a 
 pessimistic one, a possible scheme might be to keep some 
 history of the usefulness of our acceptability criterion and 
 change it accordingly. 
 
 c - llaxinuni a .1 1 own h 1 e s t e p s i z e f or an a ctual integration . It 
 seems that a large step for a certain equation can put a 
 connected equation so far behind in time from the first that 
 it can no longer converge. This would necessitate backtracking 
 and undoing the last integration which would be costly in 
 terms of storage and execution tine. Even though this has 
 been encountered only rarely (twice out of more than 100 
 Integrations completed successfully) , it would seem that a 
 mechanism for deternininp maximum possible stepsize that can 
 be talen during each Integration, without disturbing connected 
 equations irreversibly, would add to the efficiency of the 
 program* 
 
29 
 
 VI - Conclusion 
 
 The present study was undertaken to test the pronise and 
 potential of a nultirate nethod for ordinary differential 
 equations. As can be seen, the program in its present stage 
 offers a reasonably large ( 2-3 tines) gain in accuracy per 
 function evaluation. This is not strictly a gain in execution 
 tines because of the overhead in inplenenting the nultirate 
 nethod ( extra calls to subroutines, Jacobian evaluations, calls 
 to Interpolating routines before function evaluations), yet even 
 though it is difficult to quantify , there still seen to be gains 
 in accuracy per execution time. As cited in the preceding 
 section, there seens to be a nun her of questions to be studied, 
 as far as efficiently inplenenting a multirate method. 
 Therefore, a multirate method has a lot of potential and seems 
 to pronise efficiency gains in solving a variety of systems of 
 ordinary differential equations. 
 
31 
 APPENDIX 
 I - Program 
 
 The program that has been implemented is made up of 
 
 DIFSUB with subroutines that call DIFSUB and monitor the 
 
 projected integration step and report integration information* 
 
 The structure of the multirate ode integrator, MRATE , thus 
 developed is as follows : 
 
 INTERP 
 
 user-supplied 
 JACOB 
 
 user-supplied information 
 and work space 
 
 MRATE 
 
 MRATE1 
 
 INTALL 
 
 INTONE 
 
 SMALL 
 
 DFDRVR 
 
 DIFSUB 
 
32 
 
 The modified subroutine DIFSUB has the following 
 
 s t rue t ure : 
 
 PREDCT 
 
 CORECT 
 
 user-supplied 
 JACOB 
 
 DIFSUB 
 
 SHARP 
 
 SCALE 
 
 SOE 
 
 DIFFUN 
 
 UPDATE 
 
 ERRCHK 
 
 INTERP 
 
 user- 
 supplied 
 DIFUN1 
 
 The basic chnnp.es that were made to the original ode 
 Integrator ( HIFSUB ) were: 
 
 1 - As can be seen above, the ode integrator has been broken down 
 to a number of subroutines. 
 
33 
 
 2 - The ode Integrator now works on a single equation instead of 
 a system. 
 
 3 - The ode integrator works on temporary arrays with its 
 temporary results passed back to df Jrvr . Dfdrvr actually 
 decides whether to store the temporary results (if this is an 
 actual integration step) or to discard them ( if this is a 
 psGudn-intcpration) . 
 
 The subroutines of the outer package accomplish : 
 
 1 - Prate : Assigns the work space provided by the user to 
 different arrays. 
 
 2 - ll ratf 1 : Initializes standard integration parameters. Calls 
 i n t a ] 1 to do the integration over the rangn desired and prints 
 final values by interpolation. 
 
 3 - I n t a 1 1 : Main loop to do all of integration. Each loop 
 
 consists of calling intone to do one integration step, then 
 
 determining t for the variable just integrated and finally 
 p o s s J 
 
 determining t for the variables whose projected integration 
 poss r -* ° 
 
 step lias become unreliable because of the effect of the last 
 actual integration. 
 
 k - I ntone : Does one step of actual integration by a call to 
 dfdrvr and then accumulates effect of the last integration on 
 the coupled variables to the equation by updating delta array. 
 
 5 - Df'lrvr : Calls d 1 f s u b for one step of integration and 
 depending on flags supplied, stores temporary values produced 
 b y d i f s ub (if it is an actual into g ration step) or discards 
 
34 
 
 them ( If It is a pseudo- integration) . 
 The user is requested to provide : 
 
 I - A suhroutine d i f u n 1 to return the ^(y) for the i-th 
 
 variable . 
 2 - A double precision function J acob to return 
 
 8 f 
 
 3y 
 
 J 
 
 3 - Double precision work space of 15 times the number of 
 equations and integer work space of 5 times the number of 
 equations in the system. 
 
 U - A connectivity matrix of which variable influences which 
 equation. Since systems of ode's usually have sparse 
 connectivity matrices this has been implemented in list 
 structure with two integer arrays. 
 
 5 - Integration parameters like minimum step, maximum step, local 
 error tolerance, initial time and final time of integration. 
 
 Through the use of a multitude of subroutines of distinct 
 functions, the package has become easy to change and maintain. 
 Also since nothing on the package depends on the particular 
 parameters of a specific problem (i.e. variable working space) it 
 does not require any changes to suit particular problems. Finally 
 since the program has been developed and tested on an IBM 360/75, 
 CYBER 175 and INTERDATA 7/32, it has been written in standard, 
 transportable FORTRAN, and should be able to run on any machine 
 without any major modifications. 
 
ofmAEC-427 U s AT0(v1 | C ENERGY COMMISSION 
 
 (6/68) 
 AECM 3201 
 
 UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT 
 
 I See Instruction* on Rtvrm Sid* ) 
 
 AEC REPORT NO 
 
 COO-2383-0056 
 
 2. TITLE 
 
 A MULTIRATE ORDINARY DIFFERENTIAL EQUATION INTEGRATOR 
 
 l TYPE OF DOCUMENT (Check one): 
 
 K] a. Scientific and technical report 
 
 O b. Conference paper not to be published in a journal: 
 
 Title of conference 
 
 Date of conference 
 
 Exact location of conference 
 
 Sponsoring organization 
 
 □ c. Other (Specify) 
 
 I. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): 
 
 IS »• AEC's normal announcement and distribution procedures may be followed. 
 1 b. Make available only within AEC and to AEC contractors and other US. Government agencies and their contractors. 
 I I c. Make no announcement or distribution 
 
 REASON FOR RECOMMENDED RESTRICTIONS: 
 
 SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. W. Gear 
 
 Professor and Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 University of Illinois U-C 
 Urbana, IL 61801 USA 
 
 .^C&tp^ 
 
 Signature 
 
 Date 
 
 March 1979 
 
 FOR AEC USE ONLY 
 
 AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY. ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 
 RECOMMENDATION: 
 
 PATENT CLEARANCE: 
 
 J a. AEC patent clearance has been granted by responsible AEC patent group. 
 □ b. Report has been sent to responsible AEC patent group for clearance. 
 LJ c. Patent clearance not required. 
 
IBLIOGRAPHIC DATA 
 HEET 
 
 1. Kcport No. 
 
 UIUCDCS-R-79-959 
 
 3 Tit le and Subc itle 
 
 , A MULTIRATE ORDINARY DIFFERENTIAL EQUATION INTEGRATOR 
 
 3. Recipient's Accession No. 
 
 5- Report Date 
 
 March 1979 
 
 ] Author(s ) 
 
 i A. Oralloglu 
 
 8. Performing Organization Rept. 
 
 No - COO-2383-0056 
 
 1 Performing Organization Name and Address 
 
 Department of Computer Science 
 University of Illinois U-C 
 Urbana, IL 61801 USA 
 
 10. Project/Task/Work Unit No. 
 
 11. Contract /Grant No. 
 
 US ENERGY/EY-76-S-02- 
 
 2333 
 
 L Sponsoring Organization Name and Address 
 
 U. S. Department of Energy 
 Chicago Operations Office 
 9800 South Cass Avenue 
 Argonne. IL 60439 USA 
 
 13. Type of Report & Period 
 Covered 
 
 M.S. Thesis 
 
 14. 
 
 Supplementary Notes 
 
 . Abstracts 
 
 Most of the work on the numerical integration of systems of ordinary differential 
 equations concerns itself with integrating all the variables simultaneously. 
 There has been little theoretical work on integrating the variables separately 
 (multirate methods) and even less in implementing a multirate ordinary 
 differential equation integrator. In this thesis, a multirate integrator is 
 presented and different schemes are tested. The implemented multirate 
 integrator is applied to a variety of test problems and experimental results 
 are plotted. These show the possibility of large efficiency gains in using 
 multirate methods for a variety of systems of ordinary differential equations. 
 
 I. Key Words and Document Analysis. 17o. Descriptors 
 
 multirate integrator 
 
 ordinary differential equations 
 
 1»- Identifiers /Open-Ended Terms 
 
 ■- 
 
 COSATI Field/Group 
 
 Availability Statement 
 
 unlimited 
 
 rq M NTlS-38 I io-: 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIFIED 
 
 21. No. of Pages 
 36 
 
 22. Price 
 
 U3COMM-DC 40329-P7I 
 
m 3 
 
 :'H 
 
Aiir. 1 1H79