UNIVERSITY OF 
 
 ILLINOIS LIBRARY 
 
 AT URBANA-CHAMPAIGN 
 
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. 
 To renew call Telephone Center, 333-8400 
 
 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 JMR22 
 JAM 2 5 
 
 Ikii 
 
 JAN 2 1W 
 
 DEC 1 1 1996 
 
 L161— O-1096 
 
j M/f 
 
 f-2. 
 
 //-X/Cf^V • 
 
 UIUCDCS-R-80-1019 
 
 UILU-ENG 80 1721 
 C00-2383-0068 
 
 AUTOMATIC DETECTION AND TREATMENT OF OSCILLATORY 
 AND/OR STIFF ORDINARY DIFFERENTIAL EQUATIONS 
 
 by 
 
 June 1980 
 
 
 ''Mil 
 [Y*| 
 
 C. W. Gear 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/automaticdetecti1019gear 
 
UIUCDCS-R-80-1019 
 
 AUTOMATIC DETECTION AND TREATMENT OF OSCILLATORY 
 AND/OR STIFF ORDINARY DIFFERENTIAL EQUATIONS 
 
 by 
 
 C. W. Gear 
 
 June 1980 
 
 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 ENERGY/EY-76-S-02-2383, 
 Presented at the Differential Equation Workshop, Center for Interdisciplinary 
 Research (Zif ) , University of Bielefeld, West Germany, April 21, 1980. 
 
1. Introduction 
 
 A truly automatic code for ordinary differential equations must not 
 only handle the most general case reasonably efficiently, but must also 
 automatically detect those classes of problems which are unreasonably 
 expensive by general methods and switch to methods which are more efficient 
 for those problems. This paper will consider two classes of problems, 
 stiff and oscillatory, for which special methods can be far more efficient 
 than general methods. This is not to say that there are not other classes 
 of problems that are worthy of special treatment. For example, linear 
 equations probably can be solved more efficiently if this fact is known 
 [1] , but at this time there are no methods that are sufficiently more 
 efficient for linear problems than the general methods that it seems worth 
 the effort to detect linear problems automatically. (Furthermore, most 
 users can tell if a problem is linear, while it may be difficult for them 
 to tell when it becomes stiff or oscillatory.) 
 
 Although it is common to talk about "stiff differential equations," an 
 equation per se is not stiff, a particular initial value problem for that 
 equation may be stiff, in some regions, but the sizes of these regions 
 depend on the initial values and the error tolerance. For most problems 
 the solution is initially in a transient and an accurate solution demands a 
 stepsize sufficiently small that the truncation error of that transient is 
 small. For such stepsizes stability is not a problem. When the transient 
 has decayed below the error tolerance, the problem may be stiff. At this 
 time a stiff method must be used. Many techniques and programs are 
 available for stiff equations ([2], [3], [4], [5], [6]) so we will not 
 repeat that material. 
 
 Until recently stiff methods have also been used in the transient 
 region, but the fact that they are generally less efficient than nonstiff 
 methods (both because of smaller error coefficients and the linear algebra 
 involved) has encouraged several people to investigate automatic detection 
 of stiffness. 
 
 The problem of highly oscillatory solutions has some parallels to the 
 stiff problem. Again, the solution may not be nearly periodic initially, 
 but after a transient starting phase, may tend towards a periodic solution 
 
- 2 - 
 
 or have a nearly periodic behavior. There are some methods that are 
 applicable in the latter phase, for example, Mace and Thomas [9], Graff 
 [6] , Graff and Bettis [7] and Petzold [10]. However, these methods cannot 
 be used in the transient phase so it is essential to detect the time when 
 nearly periodic behavior has begun and to estimate the period reasonably 
 accurately. There are, of course, problems for which the user knows that 
 the solution is nearly periodic throughout. Satellite orbits are a case in 
 point. (Most of the early methods were developed for these problems.) In 
 such cases the period is known reasonably accurately so there is no 
 detection problem. This paper is particularly concerned with those 
 problems which may become nearly periodic in later stages and methods for 
 detecting this behavior in order to switch to an appropriate scheme. 
 
 Methods for nearly periodic problems are generally known as 
 mult Revolutionary from their celestial orbit background. The idea of such 
 methods is to calculate, by some conventional integrator, the change in the 
 solution over one orbit. If the period of an orbit is T (for a moment 
 assumed fixed), then a conventional integrator is used to compute the value 
 of 
 
 d(t) = y(t + T) - y(t) 
 by integrating the initial value problem y' = f(y) over one period T. If we 
 consider the sequence of times t = mT, m integral, we have a sequence of 
 values y(mT) which are slowly changing if y is nearly periodic. The 
 conventional integrator allows us to compute the first differences d(mT) of 
 this sequence at any time mT. Under appropriate "smoothness" conditions 
 (whatever that means for a sequence) we can interpolate or extrapolate for 
 values of d(mT) from a subset of all values of d, for example from d(kqT), 
 k = 1, 2, 3,..., where q is an integer > 1, and thus estimate y(mT) by 
 integrating only over occasional orbits. 
 
 In a satellite orbit problem it is fairly easy to define the meaning 
 of "one period." For example, one could use a zero crossing of a particular 
 coordinate, or even a fixed period based on a first order theory. In her 
 thesis, Petzold considered problems for which it is difficult to find 
 physical definitions of the period and examined a method for determining 
 the approximate period by minimizing a function of the form 
 
- 3 - 
 
 I(t, T) = | |y(T + T) - y(T) | | t 
 
 where the norm measures the values of y(r + T) - y(t) approximately over 
 the range Tf:(t, t + T). The actual norm she used was 
 
 I(t, T) = / y(T + T) - y(T) L dx 
 t z 
 
 where T was the last estimate of the period. The use of T was for 
 pragmatic reasons. Ignoring that detail, the value of T which minimizes 
 I(t, T) is a function of t, and T(t) was said to be the period of the 
 solution. This enabled d(t) = y(t + T(t)) - y(t) to be calculated and 
 multirevolutionary methods to be used. The variable period was handled 
 easily by a change of independent variables to s in which the period is 
 constant, say 1. The equation 
 
 t(s + 1) - t(s) = T(t) 
 was appended to the system 
 
 z(s + 1) - z(s) = g(s) 
 where z(s) = y(t(s)) and g(s) = d(t/s) for integer values of s. (When T is 
 constant, this is the analog of the old device for converting a non- 
 autonomous system to an autonomous system by appending the differential 
 equation t' = 1.) 
 
 The scheme for period calculation used by Petzold suffers from three 
 drawbacks. The first drawback is that it is fairly expensive, involving a 
 numerical approximation to the first two derivatives of I(t, T) by 
 quadrature which requires integration for y(x) over two periods. In the 
 experimental implementation, integration was repeated for every iteration 
 of a Newton method to minimize I(t, T) by solving 8I/8T = 0. This could 
 have been eliminated by saving all values and interpolating, but the 
 storage cost becomes high and the quadrature/interpolation cost remains 
 non-negligible. The second drawback is that a reasonably accurate period 
 estimate is needed for the Newton iteration to converge. Outside the 
 region of convergence for Newton's method a search scheme for a minimum 
 could be used but this would be very expensive because of the computation 
 involved in each quadrature even if all previously computed values could be 
 saved. This makes the approach very unattractive for initial period 
 
- 4 - 
 
 detection when there is no starting estimate. The third drawback is that 
 minimizing a function subject to several sources of error (including 
 truncation errors in the integration and quadrature, and roundoff errors 
 revealed by considerable cancellation in || |y(T + T) - y(f)||) is likely to 
 yield a fairly inaccurate answer. Since the value of d(t) = g(s) is quite 
 sensitive to small absolute changes in the period T which may be large 
 relative to the period, the function g(s) may not appear to be very smooth. 
 
 This paper discusses an alternate approach to the period 
 identification problem. It overcomes the cost and convergence problems, 
 and also seems to help with the sensitivity problem. This is discussed in 
 the next section along with stiffness detection. The third section reviews 
 multirevolutionary multistep integrators and a technique for handling 
 variable periods based on the period identification algorithm. The fourth 
 section discusses a numerical example while the final section discusses 
 unsolved problems. 
 
- 5 - 
 
 2. Periodic and Stiffness Detection 
 
 A fully automatic method should be able to detect problems with 
 special properties that can be solved more efficiently, but the cost of 
 detection should be low compared to the integration cost so that problems 
 without those properties do not cost appreciably more. Since stiffness and 
 nearly periodic behavior are properties that may appear at any point in the 
 solution, the detection process must operate continuously. If it is to 
 have a low cost, it must not take more than a few operations on available 
 intermediate or final results of a standard integrator. This section first 
 discusses periodic behavior detection, then stiffness detection. 
 
 We have been deliberately imprecise about the meaning of "nearly 
 periodic," and will continue that way with the working definition In our 
 minds of "the type of problem that can be handled efficiently by 
 multirevolutionary methods." The types of problems that have the required 
 properties are differential equations for which there exist functions 
 F(t, t) and T(t) such that 
 
 y(t) = F(t, t) 
 
 is a solution of the differential equation dy/dt = f(y, t) , F is periodic 
 in t with period T(t), that is, 
 
 F(t + T(t), t) = F(t, t), 
 
 for all t and t, 9F/9t is very small compared to 3F/8t, and T(t) is slowly 
 varying. Here, t and t are the "fast" and "slow" times as shown in Figure 
 1. The "meaning" of this representation is that 
 P(t) = {F(t, t), te[0, T(t)]} is the local periodic behavior of the 
 solution, and the change of P(t) with respect to t represents the way this 
 behavior slowly changes. (This representation is only valid for problems 
 which are not phase locked to a periodic driving function.) F(0, t) was 
 called a quasi-envelope by Petzold. It is the function z(t) defined 
 
- 6 - 
 
 earlier for a discrete set of points only 
 
 t 
 
 olution 
 
 T(t) 
 
 periodic 
 direction 
 
 Figure 1. Nearly Periodic Solution Family 
 
 This representation is not unique, but depends on the choice of the 
 period T(t) and the values of F(0, t) over the initial period. It is 
 convenient to consider a change of variables to (s, t) with s = x/T(t). In 
 the new coordinate system, F(s, t) has period 1 in s and a unique quasi- 
 envelope is defined for any fixed s in terms of F(0, t). 
 
 The "period" of a nearly periodic function has not yet been defined. 
 We could use some intuitively reasonable mathematical description, in which" 
 case we would have to seek computational algorithms for its approximation. 
 However, the period is most easily defined in terms of the algorithm used 
 to calculate it. It should, of course, yield the exact period for periodic 
 functions and be close for small perturbations of periodic functions. This 
 replaces an analysis of the accuracy of period calculation with an analysis 
 of the efficiency of the multirevolutionary method with respect to 
 different period definitions. This latter may be an easier task. 
 
 Petzold's period definition, based on minimizing the norm in eq. (1), 
 is very expensive to apply and cannot be considered as a technique for 
 determining if an arbitrary output of an integrator is nearly periodic. 
 
- 7 - 
 
 Therefore, we look for alternate definitions of the period. First, note 
 that if the oscillation is due to a periodic driving function, we probably 
 know its period or can examine the system which generates the driving 
 function directly. Hence, we can restrict ourselves to autonomous systems 
 or systems which when made autonomous are nearly periodic. (This means 
 that the partials of the system with respect to time are small. ) The 
 solution of an autonomous system is completely determined by the 
 specification of the value of the solution vector y_ at one time. That is 
 to say, if we identify two times on the solution such that y_(t^) = x( t 2^' 
 we know that the solution is periodic with period to - ti. This suggests 
 determining the period by looking for minimum of | y_(t ^ ) ~ 7.^2^11* ^he 
 cost of this is not particularly low and it requires a clever adaptive 
 program with a lot of heuristics. The form of the program is to choose an 
 initial point t, and, as each new value y^(t) is calculated, to see if 
 |y_(ti) - y_(t)| has passed a local minimum. If so, interpolation is used 
 to locate the minimum. If the value of the norm at the minimum is small 
 (first heuristic parameter), we have possible periodic behavior. If not, 
 we must continue to examine more y(t). However, since the periodic 
 behavior of y may not have started by t-., we must also advance t-, 
 occasionally (next heuristic parameter). After some experiments, we 
 abandoned this approach. 
 
 Another way of defining the period is to identify certain points on 
 the solution at which a simple characterization is repeated, such as zero 
 crossing. The solution itself may not have zero crossings and if it 
 consists of a periodic function superimposed on a slowly growing function, 
 it may be difficult to choose any value which is crossed periodically. 
 However, its derivative will have periodic sign changes, so we have 
 
 experimented with a definition of period based on the zero crossings of 
 
 T *■ T 
 
 c_ y_ where c^ is the transpose of a vector of constants. The program 
 
 T 
 examines the integrator output for positive-going zero crossings of _c y_'. 
 
 (Currently, c_ is a vector of the weights provided by the user for error 
 
 norm calculations.) Anything but a simple periodic solution may lead to 
 
 more than one zero crossing in a single period, so the norm 
 
 |y'(tj_) - y' (t2 ) I is also examined, where ti and t^ are a pair of zero 
 
 crossings. If the norm is small, the possibility of a period is 
 
- 8 - 
 
 considered. The procedure used is as follows: 
 
 T 
 
 1. Identify a positive going sign change in c y . 
 
 2. Interpolate to find the t value t t of the zero crossing. Also 
 compute interpolated values of y_, y_'. 
 
 3. Save these values. (Up to ten prior values are saved.) 
 
 4. Compare current value with each prior value in turn until a small 
 
 y'^i j - y' .„...,„.. is found. 
 i ii old J- current 1 ' 
 
 5. Save period,,^, = t _„____«. - t^ 1 A . 
 
 r new current old 
 
 6. Compare period with period ■, . if one has previously been 
 calculated. If they are relatively close, accept the new period and 
 switch to multirevolutionary methods. 
 
 7. Examine several backward differences of recent periods. If they seem 
 to be smoothly varying, accept new period. 
 
 The last test was found to be necessary for some problems with variable 
 periods. 
 
 As can be seen, there are numerous heuristics, which implies that much 
 tuning is possible. However, it is important to note that tuning effects 
 efficiency only. If the tests for periodicity are too stringent, the 
 standard integrator will continue too long; if they are too lenient, 
 multirevolutionary methods will be invoked before they are efficient. 
 However, since they will then run using a stepsize of one period, they will 
 perform very little worse then the conventional integrator. (The only 
 losses are additional overhead. ) 
 
 The multirevolutionary method, described in the next section, is a 
 modification of a standard integrator. It calls on a subroutine to 
 evaluate g(z) = z(s + 1) - z(s) given z(s). It can suffer from stiffness 
 in exactly the same way that an ordinary integrator can suffer from 
 stiffness: if H9g/9z is large the method may be unstable and the corrector 
 iteration will not converge unless a Jacobian J = 3g/8z is used in a Newton 
 iteration. Shampine [11] has suggested monitoring the size of the Jacobian 
 
- 9 - 
 
 by estimating its norm when two or more function evaluations are used in a 
 single step. Essentially, 
 
 | | g(z-, ) - g(z 2 )| 
 
 L = max — 
 
 lf Z l ' Z 2|| 
 
 is calculated, where the max is taken over all steps and z^ and z 2 are two 
 different values of z for which g is evaluated in one step. (In practice, 
 it seems preferable to take an exponentially weighted max such as 
 
 ,r, n T I I sCz, ) - g(z 2 )| I 
 
 L = max (0.9 L -,., * 4 ) 
 
 new v old' z^ - z 2 
 
 but this is yet another tuning heuristic.) 
 
 This technique can be used in both the regular integrator used to 
 calculate g and the multirevolutionary integrator. However, there are a 
 number of tuning questions that pose some difficulties. One could decide 
 to switch to a stiff method the moment that J becomes large enough to 
 restrict the stepsize in which case there are no problems. However, that 
 is not efficient because nonstiff methods are both considerably lower in 
 cost per step and have a smaller error tolerance. The natural solution to 
 consider is to continue with the nonstiff method until the estimated 
 stepsize that could be used in a stiff method is sufficiently larger to 
 offset the increased cost per step. This requires estimating the error in 
 stiff methods by estimating various derivatives. This is done with 
 suitable difference formulas but if care is not used, the derivatives 
 estimated may appear to be such that little step increase is possible with 
 a stiff method. The reason for this difficulty is that most codes do not 
 directly restrict their stepsize on the basis of a Jacobian estimate and 
 stability needs. At the most they restrict h||j| so that the corrector 
 would converge if iterated. However, higher order methods may be on the 
 verge of instability. Once h||j| is too large for stability, small errors 
 grow rapidly. The automatic error control quickly reduces the step to keep 
 the error near the tolerance level and a well-engineered nonstiff solver 
 will produce a perfectly good solution, albeit slowly. However, the 
 solution contains errors of the order of the error tolerance due to these 
 marginally stable components, and these errors usually oscillate. When a 
 
- 10 - 
 
 difference formula is applied to them, large values result, and completely 
 obscure the derivatives we want to estimate. For example, the marginal 
 stability could introduce an error of (-l) n e into the numerical solution y 
 
 at the n-th step. If we form the k-th backward difference to estimate 
 
 k (k^ n k 
 
 h y v ' we have a component due to this error of (-1) 2 e. If we now ask 
 
 k fk") 
 what stepsize ah can be used in a stiff method whose error is C, h y v ' to 
 
 achieve an error of t, we find that a = (2 Ci ) ' independent of current 
 
 h. For BDF, C^ = 1/k so a is always less than one, falsely indicating that 
 
 the stepsize cannot be increased. To avoid this difficulty it is necessary 
 
 to keep h||j| small enough that components with the most negative 
 
 eigenvalues are at least moderately damped. 
 
 In addition to estimating |j||, it is possible to estimate the 
 largest eigenvalue (or eigenvalue pair if complex conjugate) using 
 evaluations of g(z) from more than one step, as suggested in Gear [12]. 
 However, experiments on that technique indicate that the eigenvalue 
 estimates are not too reliable. The reason for wanting to know the largest 
 eigenvalue is to know whether to use a technique efficient for real 
 eigenvalues, such as BDF, or a technique better suited to eigenvalues close 
 to the imaginary axis, such as BLEND [13] or implicit Runge-Kutta [3]. K. 
 Stewart* pointed out that it is sufficient to wait until a decision to use 
 stiff methods has been made. 
 
 At that time a Jacobian must be calculated, and power methods can be used 
 to determine the arguments of the large eigenvalues. (This poses an 
 interesting question for the numerical linear algebraist: how to calculate 
 cheaply the maximum argument of all eigenvalues exceeding a certain size in 
 a matrix. ) At the time the Jacobian is first calculated, it can also be 
 checked for other properties such as handedness and sparsity so that a 
 decision on which linear equation scheme to use can be made. 
 
 *K. Stewart, Jet fropulsion Lab, Pasadena, CA. Private communication. 
 
- 11 - 
 
 3. Variable Period Mult Revolutionary Methods 
 
 In the original coordinates we have z(t + T) - z(t) ■ d(t). For small 
 T this says z'(t) = d(t)/T. Hence, it is not surprising that the numerical 
 interpolation for z(t) given a technique for computing g(t) ■ d(t)/T is 
 very similar to a numerical integration technique. In the new coordinate 
 system, the basic structure of the program is an outer integrator which 
 solves the equations 
 
 z(s + 1) - z(s) - g(z) 
 t(s + 1) - t(s) = T(t(s)) 
 using an outer stepsize H. The method varies the order and stepsize just 
 as an ordinary integrator does. See Petzold [10] for details. It calls a 
 subroutine to evaluate g and T given z and t. This is done by integrating 
 the underlying ordinary differential equation y' = f(y) starting from y(t) 
 = z, determining when a period has elapsed and computing g(z) - y(t + T(t)) 
 y(t). Both methods for defining the period discussed in the pevious 
 section have been tried. The first, minimizing | |y(t + T(t)) - y(t)|j, is 
 now easier to implement because the left end, y(t), is fixed. The only 
 tuning difficulty is to ignore intermediate minima, and we have done this 
 by considering only values of T starting from 0.9 of the previous period 
 estimate. (If T changes more rapidly over H than this, either H is too 
 large or the nearly periodic assumption is questionable.) The norm actually 
 used has the form 
 
 ' i " to I A - ( * )): 
 
 where vi^ is the j-th derivative of the i-th component, and the A . are 
 weights. It appears to be best to use A . =0, s ^ 1, and A^., = weight of 
 i-th component of error weight vector. This allows for arbitrary 
 nonperiodic linear functions to be ignored. Higher derivatives can be 
 Included, but knowledge of them is subject to larger errors due to the 
 inner integration. 
 
 The second method, looking for a zero crossing, has a difficulty: the 
 
 T ~ 
 function c y_' will not necessarily be zero at y(t) ■ z. (It can be shown 
 
 that it will be zero except for roundoff error for linear problems.) This 
 
- 12 - 
 
 has been overcome by choosing a vector c^ separately for each period such 
 that p(t) = c j_' = at the start of the period. A future zero of p(t) 
 defines the length of the period. A unique c_ is defined by choosing c_ to 
 maximize p'(t) at the start of the period subject to |cj ■ 1. This value 
 is chosen because it minimizes the roundoff error effects in determining a 
 zero of p(t). The value of £, apart from scaling, is 
 
 £ = z " - ft i, 7 ,2* 1 
 
 This requires a knowledge of y" at the initial point y(t). This is 
 available because a Runge-Kutta starter which computes the first four 
 derivatives is used for the multistep inner integrator [14]. After c_ has 
 been calculated, future positive going zero crossings of p are examined, y_ 
 and its derivatives are calculated by interpolation at the zero crossing 
 point t + t, and |y_'(t + t) - y_'(t)| is evaluated. If it is small 
 enough, the period is set to t and g is calculated. 
 
 The variable period multirevolutionary integrator is based on a 
 modified Nordsieck scheme. Each component of z is represented by the 
 history vector 
 
 a = [z, Hg, H 2 g'/2, H 3 g"/6,..., H^" 1 ^! ] T 
 
 Petzold has shown that in this representation the predictor has the form 
 
 ln,(0) = ^-1 
 where A is the Pascal triangle matrix except for the first row which is 
 
 [1, 1, a L (r), a 2 (r),..., ot^Cr)] 
 
 where r = 1/H. She also showed that the corrector takes the form 
 
 £n = ^n,(0) + l u 
 
 where m is chosen so that a„ "satisfies" the relation 
 
 — n 
 
 z(s n + 1) - z(s n ) = g(s n ) and %_ is the conventional corrector vector except 
 in the first component which is a function of r = 1/H. Petzold gives these 
 functions for generalized Adams methods. (They are polynomials in r. ) The 
 corresponding functions for BDF methods are inverse polynomials in r. They 
 
- 13 - 
 
 are 
 
 Order First Coefficient of l_ 
 
 1 -1 
 
 2 -1/(3/2 + r) 
 
 3 -1/(11/6 + 2r + r 2 ) 
 
 4 -1/(25/12 + 35r/12 + 5r/2 + r 3 ) 
 
 5 -1/(137/60 + 15r/4 + 17r 2 /4 + 3r 3 + r 4 ) 
 
 6 -1/(147/60 + 203r/45 + 49r 2 /8 + 103r 3 /12 + 7r 4 /2 + r 5 ) 
 
 Petzold suggests a linear combination of the generalized Adams and BDF 
 coefficients, for example, r • Adams + (1 - r) • BDF so that the method has 
 the properties of BDF formula for large H (r -> 0), namely stiff stability, 
 and the property of Generalized Adams for r = 1, namely the outer 
 integrator is exact. Since it is not proposed to use BDF methods until 
 stiffness has set in (when H is large), it does not seem worth considering 
 this complication. 
 
- 14 - 
 
 4. A Numerical Test 
 
 Several example problems have been constructed using the Van der Pol 
 oscillator to give a nonlinear oscillation. Typical of these problems is 
 the following system of four equations 
 
 ii " Q L 
 
 "l " U' 
 
 = -(u, - Uo) + 2(uo - (u, - Uo) )u. 
 
 u 3 
 u 4 
 
 -10" 3 (u 3 - 1) 
 10" 3 sin 10" 3 t 
 
 All initial values were zero so uo = 1 - e~ and ua = 1 - cos .OOlt. 
 u^ and u£ are the solution and first derivative of a Van der Pol oscillator 
 oscillation about a level U3 and peak amplitude about 2uo. The period is 
 about 2tt for small Ui and steadily increases to about 7.63 for u<% = 1. 
 
 The matrix Q used was 
 -1 111 
 1-111 
 11-11 
 111-1 
 
 1/2 
 
 ,-1 
 
 (It is idempotent and Q = Q .) All components of y oscillate after an 
 initial period. The periodic detector located the oscillation at about t = 
 156. The integration was continued to t = 10,000 with an average outer 
 step of about 28 periods. At that stepsize the outer problem is quite 
 stiff. The oscillatory behavior is initially close to sinusoidal and 
 changes to the steep-edged behavior typical of the Van der Pol oscillator 
 by t = 1000. A local error tolerance of 10 in the inner integrator 
 required about 400 inner steps per period at first, increasing to about 
 
 1200 at the end. The outer integrator took 50 steps with local tolerance 
 
 _3 
 of 10 from t » 156 to t = 10,000, using 154 inner integrations over one 
 
 period including those for occasional Jacobian evaluations of g by 
 
 numerical differencing, an average speed up of ninefold over the standard 
 
- 15 - 
 
 inner integrator which would have used about 10 steps for the whole 
 problem. 
 
 Plots of the phases of the solution are shown in Figures 2 to 5. In 
 all of these figures the vertical scales for the four components have been 
 renormalized to put y^ between 2i - 1.9 and 21 - 0.1 for I = 1 to 4. 
 Figure 2 shows the first phase of the integration prior to detection of the 
 oscillation. 
 
 Figure 2. Initial Phase before Period Detected 
 (For extraneous reasons, only one integration point in ten has been plotted 
 here, hence the jagged curves.) The amplitude of the oscillation at this 
 point is 0.99. The shape of the oscillation at t = 156 is shown in Figure 
 3. This shape changes and grows in amplitude to 3.02 by t = 10,000, as is 
 shown in Figure 4. Figure 5 shows the smooth "quasi-envelope" z found by 
 the multirevolutionary integrator. It was generated using the 50 outer 
 integration steps so the actual solution y is found by superimposing the 
 oscillatory behavior of the form shown in Figures 3 and 4 and the 
 
- 16 - 
 
 appropriate t values in Figure 5. 
 
 i .i 
 
 iiHi 
 
 7.1 
 
 e. t 
 
 6. I 
 
 5. 1 
 
 Figure 3. First Period in Multirevolutionary Integration 
 
 Figure 4. Fiftieth (last) Period in Multirevolutionary Integration 
 
- 17 - 
 
 a. 2»m. 4«m. saaa. saaa. laaaa. 
 
 isaa. sa*a. ca«a. 7eee. ^aaa. iiaee. 
 
 Figure 5. Quasi-envelope z 
 
 0. aaaa. 4000. aaaa. saaa. 10000. 
 
 iaaa. sbbb caae . ?aaa. 9000. 11000. 
 
 Figure 6. Quasi-envelope Corresponding to u (Qz) 
 It should be noted that these are not small oscillations. For example, the 
 range of z, in Figure 5 (bottom line) is -1.13 to -0.6 (approx. ). The 
 
- 18 - 
 
 oscillation changes in amplitude from 0.99 to 3.02 (peak to peak) over the 
 interval. 
 
 In order to check the accuracy, the equivalent quasi-envelope for _u 
 was recovered from £ by the transsformation _u ■ Qz_. Since uo and ut are 
 
 not oscillatory, G3 and u^ should also be 1 - exp(-10 t) and 
 
 1 - cos(10 t), respectively. Since the cosine component is neutrally 
 
 stable, any integration errors will not be damped in later steps. The 
 
 relative error in the cosine component at t = 10,000 was .005 (.008 
 absolute) . 
 
 These results were without tuning. There are a number of 
 inefficiencies in the software that can be removed. For example, we did 
 not give the multirevolutionary integrator the information gathered during 
 the period detection phase about g and its differences. These can be used 
 to allow a high order start in the multirevolutionary integrator. We 
 believe that additional improvements are possible. 
 
- 19 - 
 
 5. Further Problems 
 
 There are a number of additional problems of concern. Here we will 
 discuss three problems: non-autonomous problems, detecting the end of 
 periodicity, and the multiple oscillator problem. Some problems require 
 only simple extensions; others, in particular the multiple oscillator 
 problem, pose serious difficulties. 
 
 There are two cases of the non-autonomous problem y" ■ f(t, y) to 
 consider: those in which f(t, y) is a slowly changing function of t and 
 those in which the t-dependence is responsible for the oscillation — we say 
 that the oscillation is driven . In the former case we can conceptually 
 convert to an autonomous system by appending the usual t" = 1. This term 
 is slowly varying so the solution of the enlarged problem remains "nearly 
 periodic. " In the latter case we can determine the period by examining the 
 driving term (that is, the t-dependent terms in f) and continue to use the 
 same method. 
 
 Some nonlinear oscilltors are such that a variable in the system 
 increases to a point that the oscillation is quenched. Then there is a 
 period of relaxation until it starts again. For such systems an automatic 
 program must detect the end of the oscillation and revert to a conventional 
 method. This is analogous to the problem of detecting a derivative 
 discontinuity in the solution of a conventional differential equation and 
 similar techniques can be used. When the multi-revolutionary integrator 
 calls for an evaluation of the local period and of g(z), the period 
 detection scheme will be unable to find a period anywhere close to prior 
 value, or even to find one at all. After it has looked a modest distance 
 beyond the expected value, it should report failure so that the 
 multirevolutionary integrator can reduce its stepsize to find the 
 "discontinuity" where the oscillator begins to be quenched. When the 
 stepsize has been reduced to one of only a few periods, the software can 
 revert to a conventional integration method. 
 
 The multiple oscillator problem poses difficulties unless there is a 
 large gap in frequencies between the two highest frequencies. In that 
 case, the lower frequencies can be viewed as the slowly changing components 
 
- 20 - 
 
 of F(t, t) and possibly the method can be used recursively for the second 
 highest frequency, and so on. If the two highest frequencies w, and o>2 are 
 of the same order of magnitude, the behavior will be far from nearly 
 periodic unless u^/wo ■ p/q for small integers p and q. In that case there 
 is a subharmonic w^/p of the two frequencies which can be used as the 
 period. If not, there does not seem to be much hope unless the oscillators 
 can be isolated and treated separately by the techniques discussed above. 
 Suppose we can visualize the system as consisting of two oscillators u' = 
 P(u, y) and v" ■ q(v, y) where y is a slowly varying term, and a slow part 
 described by y" = f(y, u, v). If f is linear in u and v, it is sufficient 
 to find the behavior of the average of u and v and this can be done for 
 each separately. However, if f is nonlinear in u and v, we must also keep 
 track of the relative phases of the oscillations of u and v so that each 
 time f is evaluated on a coarse mesh, the correct relative phases can be 
 used. 
 
 References 
 
 [1] Shampine, L.F., Linear equations in general purpose codes for stiff 
 ODEs, Report 80-0429 Sandia Laboratories, Albuquerque, NM, February 
 1980. 
 
 [2] Butcher, J.C., A transformed implicit Runge-Kutta method, Report 111, 
 Dept. Mathematics, Univ. Auckland, New Zealand, May 1977. 
 
 [3] Butcher, J.C., Burrage, K. and F.H. Chipman, STRIDE: Stable Runge- 
 Kutta Integrator for Differential Equations, Report 150, Dept. 
 Mathematics, Univ. Auckland, New Zealand, August 1979. 
 
 [4] Gear, C.W., The automatic integration of stiff ordinary differential 
 equations, Proceedings IFIP Congress 1968 , 1968 , 187-193 . 
 
 [5] Hindmarsh, A., GEAR: ordinary differential equation solver, Report 
 UCID-30001, Rev. 3, Lawrence Livermore Laboratory, CA, 1974. 
 
 [6] Byrne, G.D. and A.C. Hindmarsh, EPIS0DEB: An Experimental Package for 
 the Integration of Systems of Ordinary Differential Equations with 
 Banded Jacobians, Report UCID-30132, Lawrence Livermore Laboratory, 
 
- 21 - 
 
 CA, April 1976. 
 
 [7] Graff, O.F., Methods of orbit computation with multirevolution steps, 
 Report AMRL 1063, Applied Mechanics Research Laboratory, Univ. Texas 
 at Austin, TX, 1973. 
 
 [8] Graff, O.F. and D.G. Bettis, Modified multirevolution integration 
 methods for satellite orbit computation, Celestial Mechanics 11 , 1975, 
 443-448. 
 
 [9] Mace, D. and L.H. Thomas, An extrapolation method for stepping the 
 calculations of the orbit of an artificial satellite several 
 revolutions ahead at a time, Astronomical Journal 65 (5), June 1960. 
 
 [10] Petzold, L.R., An efficient numerical method for highly oscillatory 
 ordinary differential equations, Report UIUCDCS-R-78-933, Dept. Comp. 
 Sci., Univ. Illinois at Urbana-Champaign, IL, August 1978. 
 
 [11] Shampine, L.F., Lipschitz constants and robust ODE codes, Report 
 79-0458, Sandia Laboratories, Albuquerque, NM, March 1979. 
 
 [12] Gear, C.W., Method and initial stepsize selection in multistep ODE 
 solvers, Report UIUCDCS-R-80-1006, Dept. Comp. Sci., Univ. Illinois 
 at Urbana-Champaign, IL, February 1980. 
 
 [13] Skeel, R.D. and A.K. Kong, Blended linear multistep methods, ACM 
 Trans . Math . Software 3 (4), December 1977, 326-345. 
 
 [14] Gear, C.W., Runge-Kutta starters for multistep methods, Report 
 UIUCDCS-R-78-938, Dept. Comp. Sci., Univ. Illinois at Urbana- 
 Champaign, IL, September 1978, to appear ACM Transactions on 
 Mathematical Software. 
 
FormAEC-427 U.S. ATOMIC ENERGY COMMISSION 
 
 J?' 6 ?' UNIVERSITY -TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIF'C AND TECHNICAL DOCUMENT 
 
 I See Instruction* on Raver* Side ) 
 
 I. AEC REPORT NO. 
 
 COO-2383-0068 
 
 2. TITLE 
 
 AUTOMATIC DETECTION AND TREATMENT OF OSCILLATORY 
 AND/OR STIFF ORDINARY DIFFERENTIAL EQUATIONS 
 
 3. TYPE OF DOCUMENT (Check one): 
 
 CCI a. Scientific and technical report 
 
 I I 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) 
 
 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): 
 
 Kl a. AEC's normal announcement and distribution procedures may be followed. 
 
 I I b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 "2 c. Make no announcement or distrubution. 
 
 5. REASON FOR RECOMMENDED RESTRICTIONS: 
 
 6. SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. W. Gear 
 
 Professor and Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, IL 61801 
 
 Signature . ) 
 
 Date 
 
 June 1980 
 
 FOR AEC USE ONLY 
 
 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 8. PATENT CLEARANCE: 
 
 D a. AEC patent clearance has been granted by responsible AEC patent group. 
 O b. Report has been sent to responsible AEC patent group for clearance. 
 I~~l c. Patent clearance not required. 
 
BIBLIOGRAPHIC DATA 
 SHEET 
 
 1. Report No. 
 
 UIUCDCS-R-80-1019 
 
 4. Title and Subtitle 
 
 AUTOMATIC DETECTION AND TREATMENT OF OSCILLATORY 
 AND/OR STIFF ORDINARY DIFFERENTIAL EQUATIONS 
 
 3. Recipient's Accession No. 
 
 5. Report Date 
 
 June 1980 
 
 6. 
 
 7. Author(s) 
 
 C. W. Gear 
 
 8. Performing Organization Rept. 
 
 '' UIUCDCS-R-80-1019 
 
 9. Performing Organization Name and Address 
 
 Department of Computer Science 
 University of Illinois U-C 
 Urbana, IL 61801 
 
 10. Project/Task/Work Unit No. 
 
 11. Contract/Grant No. 
 
 ENERGY/EY-76-S-02-238: 
 
 12. Sponsoring Organization Name and Address 
 
 U.S. Department of Energy 
 Chicago Operations Office 
 9800 S. Cass Avenue 
 Argonne, Illinois 
 
 13. Type of Report & Period 
 Covered 
 
 technical 
 
 14. 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 The next generation of ODE osftware can be expected to detect special problems and to 
 adapt to their needs. This paper is principally concerned with the low-cost, automtic 
 detection of oscillatory behavior, the determination of its period, and methods for its 
 subsequent efficient integration. It also discusses stiffness detection. In the first 
 phase, the method for oscillatory problems discussed examines the output of any inte- 
 grator to determine if the output is nearly periodic. At the point this answer is 
 positive, the second phase is entered and an automatic, nonstiff, multirevolutionary 
 method is invoked. This requires the occasional solution of a nearly periodic initial 
 value problem over one period by a standard method and the re-determination of its 
 period. Because the multirevolutionary method uses a very large step, the problem has 
 a high probability of being stiff in this second phase. Hence, it is important to 
 detect if stiffness is present so an appropriate stiff, multirevolutionary method can 
 be selected. 
 
 17. Key Words and Document Analysis. 17o. Descriptors 
 
 multirate 
 integration 
 differential ODEs 
 
 17b. Identifiers/Open-Ended Terms 
 
 17e. COSATI Field/Group 
 
 18. Availability Statement 
 
 unlimited 
 
 FORM NTIS-38 (10-70) 
 
 19. Security Class (This 
 Report) 
 
 _„ UNCLASSIFIED 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIFIED 
 
 21. No. of Pages 
 25 
 
 22. Price 
 
 USCOMM-DC 40329-P71 
 
* 
 
FEB 1