UNIVERSITY OF 
 
 ILLINOIS LIBRARY 
 
 AT URBANA-CHAMPAIG.M 
 
U»»v.r.|^ ,,CMO,, — "*y re,„,» , n °dl,r T^— * 
 tcucu, -'. '" "«ni«*ol from 
 
 "••ft, 
 
 * or Ascipiino, 
 
 ^S^-r- 
 
 Feb 2 
 
 1 
 
 
 KC0 72010 
 
 ^en renewing by Dhnn. 
 
 P<**ous due ^^^"^^due date below 
 
 L162 
 
UIUCDCS-R-79-983 
 
 3 10, S^ 
 JLJtCr 
 
 no. 9X3 
 
 MccVVt 
 
 UILU-ENG 79 1731 
 COO-2383-0061 
 
 October 1979 
 
 ODEs: Is there anything left to do? 
 
 by 
 
 C. W. Gear 
 
 • •; 
 
 z*i 
 
 1*2*1 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 
 
 UNIVERSITY Or H 
 
 MPQANA-^Hy-- 
 
UIUCDCS-R-79-983 
 
 ODEs: Is there anything left to do? 
 
 by 
 
 C. W. Gear 
 
 October 1979 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 61801 
 
 Supported in part by the U.S. Department of Energy under grant ENERGY/EY-76-S-02-2383 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/odesisthereanyth983gear 
 
ODEs : Is there anything left to do ? 
 
 C. W. Gear 
 
 Department of Computer Science 
 
 University of Illinois 
 
 Urbana, IL 61801 
 
 J.. INTRODUCTION 
 
 In 1962 Gene Golub questioned why I was starting to work on the 
 numerical ODE problem. "Haven't all the important problems been solved?" 
 he asked. Indeed, it seemed at the time that they had. Dahlquist's 
 important paper [8] had appeared 6 years earlier and solved what seemed to 
 be all of the important and interesting theoretical questions; most 
 numerical analysts were unaware of any major classes of ordinary 
 differential equations that could not be solved by computer at speeds that 
 seemed to be reasonable for the difficulty of the problem; and the 
 available computer codes seemed to perform fairly well. At the most it 
 appeared that there were a few "clean-up" questions to take care of. 
 However, at that time, the problem of stiffness was already known to some 
 people. It appears to have first appeared in 1952 (Curtiss and 
 Hirschfelder [7]) and Dahlquist was about to publish his classic 1963 paper 
 [9] on it. Ten years later, reasonably robust software began to appear, 
 but even today there are still several important theoretical questions 
 waiting to be answered that may have significant practical implications for 
 stiff problems. 
 
- 2 
 
 Apart from these questions, is there anything else to do in the field? 
 This paper examines the steps in the development of a typical problem from 
 the time it is first identified by the "user" (that is, by the 
 scientist/engineer who runs into an intractable computational problem) to 
 the time that robust software is developed for classes of problems. This 
 process may initially take 20 years; during the early years very few will 
 be aware of the problem, and, even at the end after robust software has 
 appeared, development may continue. After examination of the development 
 process, the paper looks at a number of problems that are in various stages 
 of development and reviews recent results in these areas . We will see that 
 there are a number of problems of great interest at all stages of 
 development and that not only have there been a number of important new 
 results recently, but that there are many open problems. 
 
 _1 ._1. The Development Staircase 
 
 The progress from problem identification to robust software can be 
 broken into six levels. These are: 
 
 1 Identification of problem by a user. 
 
 2 Development of a method and preliminary codes. 
 
 3 Mathematical analysis of the problem and method. 
 
 4 Preparation of initial production codes. 
 
 5 Analysis of software for classes of problems. 
 
 6 General purpose, robust software development. 
 
- 3 - 
 
 In the case of stiff ODEs, these developments took two decades from the 
 early 1950's. Today, the solution of most stiff problems is relatively 
 routine unless additional problems (such as very large systems or many 
 discontinuities) are present. However, the development has not stopped. 
 New methods and new analyses continue to appear, giving rise to the 
 "staircase" picture shown in Figure 1. Each level of the staircase 
 continues to spawn new developments at higher levels. 
 
 The stiff problem has already reached the top of this staircase, so, 
 once again, we might ask "Is there anything left to do?" besides a few 
 clean-up questions. This paper examines a number of problems, some of 
 which are at or near the top of the staircase, but others of which are only 
 on the first step. Variable-order, variable-step methods, for example (or, 
 more precisely, problems for which variable-order, variable-step methods 
 are essential for efficient integration) are also on the top level, but 
 problems which require low accuracy methods are at the first step: they 
 have been identified but little has been done for their solution. The 
 problem areas we will examine are : 
 
 Stiff problems 
 
 Variable-order, variable-step methods 
 
 Problems with eigenvalues near the imaginary axis 
 
 Oscillating problems 
 
 Very large stiff problems 
 
- 4 - 
 
 Problems with frequent discontinuities 
 
 Multirate methods 
 
 Low accuracy methods 
 
 The approximate status of these problems is shown in Figure 2. This figure 
 does not include any of the problems in ODEs arising from PDEs such as 
 variable-sized systems that occur when the method of lines is used with a 
 varying space discretization, and the special structure of the Jacobian 
 that arises from the method of lines which may make very large problems 
 more tractable than in the general case. If we had a crystal ball, we 
 could extend Figure 2 to see what problems we will be facing in the 
 future — but that would automatically put those problems on step 1 by virtue 
 of their having been identified! However, the fact that there are problems 
 at all levels suggests that there will no shortage of problems in the 
 future. 
 
 _2. Problems showing recent progress 
 
 This section examines some results of the last few years that have 
 contributed to our understanding of some of the problems and methods for 
 solving them. The results mentioned are frequently applicable to several 
 of the problems; in the space available we will discuss their application 
 only to one problem. 
 
 2«_1. Stiff problems 
 
 The major outstanding problem of the last few years bas been the error 
 and stability analysis of stiff, non-linear systems. Burrage and Butcher 
 
- 5 - 
 
 [4] and Burrage [3] have examined the stability of non-linear, non- 
 autonomous systems for Runge-Kutta methods. They have derived a sufficient 
 algebraic condition (i.e., a condition on the coefficients of the method) 
 for stability of the numerical solution to problems which satisfy 
 
 (2.1) Re<f(y,t) - f(z,t),y - z> < 
 
 Note that this condition implies the stability of the differential system. 
 Stability of the method is defined to be the condition 
 
 (2.2) Hy - z || < ||y - z || 
 
 n n n— 1 n-l 
 
 This condition implies that the difference between numerical solutions 
 
 reduces at each step. A similar approach is now being used for multistep 
 
 method analysis by Dahlquist [10], [11], and Nevalinna and Liniger [19] who 
 
 use eq. (2.2) with y and z interpreted as vectors of all past values used 
 
 in computing the forward solution. This is a much stronger condition than 
 
 absolute stability which only requires that | |y - z | be bounded as n 
 
 n n 
 
 goes to infinity. The former condition is called contractivity . The 
 advantage of working with contractivity is that it is possible to examine 
 each step independently. The difficulty with the analysis of stability of 
 multistep methods is that adjacent steps interact, so that non-linearities 
 in the problem or variable stepsizes give rise to products of changing 
 operators for which there may not be a single fixed norm in which they are 
 all bounded by one although the product may be bounded . When we have 
 contractivity in some fixed norm for the range of problems and methods 
 used, it automatically becomes a global property. 
 
 The analysis has been greatly aided by the introduction of one-leg 
 
- 6 - 
 
 methods by Dahlquist (see [10]). These are intimately related to multistep 
 methods, taking the form 
 
 k k 
 
 (2.3) I ay + hf( I By ) - 
 1-0 x n ' i i=0 x n_i 
 
 In fact, for fixed-step methods, these methods are equivalent to multistep 
 
 methods with the same coefficients. However, they differ for variable-step 
 
 methods. The one-leg methods are easier to analyse and probably have 
 
 better stability properties. The ease of analysis arises from the fact 
 
 that each step is governed by a difference formula in which only one value 
 
 of the function f is computed. Contractivity can be studied by local 
 
 linearization of the problem and there is only one Jacobian to deal with in 
 
 a linearization of a one-leg method. For example, if we apply the 
 
 trapezoidal rule to stepsize h , . /0 = t , . - t , namely, 
 
 n+l/z n+1 n 
 
 y .i = y +h ,, ,„(f ,. + f )/2, we find that contractivity requires 
 n+1 nn+l/2n+ln 
 
 II" -k + l/2 J n + l 1 "' tI+ i t, n + l/2 J n ] ll ' l 
 
 where J and J ., are the Jacobians at T and t ,, respectively. Since J 
 n n+1 n n+1 
 
 may change over the interval, this bound may be very difficult to 
 
 establish; indeed it may not be true. On the other hand, the equivalent 
 
 one-leg method, namely y ,, = y + h ,, y „f((y ., + y x /0 \ leads to a 
 
 n+1 n n+1/2 n+1 n)/2) 
 
 requirement that 
 
 i |[I --k + i/2 jrl(I+ i\, + i/2 J! " * > 
 
 where J is evaluated at t . . ,_ in both instances in the last equation. 
 
 n+1/2 
 
 Dahlquist [10] has shown that if the problem satisfies eq. (2.1), a 
 
- 7 - 
 
 fixed step, constant coefficient, A-stahle multistep method is such that 
 tbere exists a norm in which the method is contractive and conversely, 
 ''ence, A-stahility is an attractive property for multistep methods. 
 Unfortunately, A stability limits the order to two so there is considerable 
 interest in generating similar results for less stringent stability 
 conditions such as A(a) and stiff stability. Somewhat surprisingly, the 
 algebraic stability condition of Burrage and Butcher [4] is stronger than 
 A-stabilitv, meaning that an A-stable PK method may not be contractive for 
 a problem satisfying eq. (2.1). Since contractivity and stability of one- 
 step methods are essentially identical (a method that is not contractive in 
 any steps is hound to he unstable), an A-stable RK method may not be 
 suitable for a non-linear stiff equation. 
 
 _?._?. Problems with imaginary eigenvalues 
 
 If the system Jacohian has eigenvalues near the imaginary axis and the 
 corresponding components are known to be absent from the solution, one is 
 interested in a method which damps such components numerically, 
 independently of the stepsize used. Since absolute stability is desired 
 for eigenvalues with more negative real parts, one is forced to ask for A- 
 stability (unless a special nethod that does exponential fitting is used). 
 'Hie limitation on order of multistep methods that are A-stahle forces 
 consideration of other methods. Two sets of recent results are 
 particularly significant here. A charminp paper by banner, w a irer, and 
 Norsett f?5] answers many open questions about the maximum order of a 
 variety of A-stable methods and gives simple proofs for many Vnown results. 
 In narticular, it is shown that the maximum order of an A-stable method 
 
- 8 - 
 
 cannot exceed 2q, where q is the number of derivatives or stages used. 
 This means that the choice for higher-order A-stable methods is restricted 
 to highly implicit methods such as implicit Runge-Kutta methods and multi- 
 derivative methods. The latter are not very attractive because the 
 Jacobian must be computed exactly. (They may be suitable for problems in 
 which the Jacobian is known explicitly or is fixed.) Implicit Runge-Kutta 
 methods are potential candidates, but the size of the non-linear problem 
 that has to be solved is a problem. For systems of s equations using a q- 
 stage method, a simultaneous system of sq equations must be solved. Work 
 by Bickart [2] and Butcher [6] shows how to implement implicit Runge-Kutta 
 methods in a relatively efficient manner but a more promising approach is 
 to use semi-implicit Runge-Kutta methods with equal diagonal elements so 
 that only one LU decomposition is needed along with q back substitutions 
 for a q-stage method. Unfortunately, these q-stage methods are limited to 
 order q+1. Another approach being explored by Butcher [5] is based on 
 iterated defect correction, and also requires the solution of a set of s 
 non-linear equations. Yet another scheme is the blended multistep method 
 of Skeel and Kong [23]. This is a linear combination of two multistep 
 methods, where the combination depends on the Jacobian. The accuracy is 
 derived from the methods and is not dependent on the accuracy of the 
 Jacobian. On the other hand, if the stability for linear problems is 
 analyzed using the true Jacobian, the combined method can have the 
 stability properties of a multiderivative method. The multiderivative and 
 equivalent blended methods require the solution of a linear equation 
 involving a polynomial in the Jacobian. If this polynomial is chosen to 
 have equal roots, an LU factorization of a linear polynomial can be 
 performed, followed by multiple back substitutions. This makes maximum use 
 
- 9 - 
 
 of the sparsity by avoiding powers of the Jacobian. 
 
 2.3. Oscillating Problems 
 
 We distinguish between problems with almost purely imaginary 
 eigenvalues corresponding to components not present in the solution and 
 problems with oscillating solutions. In the former, the problem is to 
 avoid introducing the oscillating components. Since difficulties arise in 
 the former case only when the imaginary parts are very large, these have 
 been viewed in the same category as stiff problems, and this seems to be a 
 reasonable viewpoint. (See Miranker and Wahba [18].) If, however, the true 
 solution exhibits very high frequency oscillations, the problem can be very 
 different. Only if the effect of the oscillations is linear is it 
 reasonable to damp them and seek the average solution. In this restricted 
 case, it is still reasonable to view the problem as a stiff one because the 
 objective of the numerical method is the same, _to ignore components due to 
 large eigenvalues . However, if the effect of the oscillation on the rest 
 of the system is non-linear, it is not possible to damp the oscillation 
 without error. In this more general case, it is not reasonable to view it 
 as a stiff problem; rather we must ask what it is we want to determine 
 about the solution. Some recent work by Petzold and Gear [22] and Petzold 
 [21] addresses this issue. 
 
 Suppose we want to know the waveform after a very large number of 
 cycles, or the envelope of the solution (see Figure 3). Unfortunately, the 
 waveform is not defined for an initial value problem because the initial 
 value defines a unique solution. However, the "envelope" in Figure 3 is 
 apparent to any engineer! To make the problem mathematically (and 
 
- 10 - 
 
 numerically) tractable, we define a quasi envelope . This is a smooth curve 
 which can be "integrated" using stepsizes which may be very much larger 
 than the period of the oscillation T. The definition of the quasi envelope 
 is illustrated in Figure A. If the period T is specified for an autonomous 
 differential equation, the quasi-envelope z(t) is defined in terms of an 
 arbitrary initial segment for t in [0,T) by integration over single 
 periods. Thus, the value of z at s+t is defined in terms of z(s) by the 
 formula 
 
 z(s+T) = z(s) + (v(s+T) - v(s)) 
 
 where v is the solution of the differential equation with initial values 
 (s,z(s)). Clearly, if z(0) = y(0), z agrees with the solution of the 
 initial problem at all integral multiples of T. We will write the above 
 equation in the form 
 
 (2. A) z(t+T) = z(t) + Tg(t,z(t)) 
 
 where 
 
 (2.5) g(t,z) = (v(t+T) - v(t))/T 
 
 Equation (2. A) looks like the Euler method applied to the differential 
 equation z' = g(t,z), and this idea is exploited in developing 
 interpolation methods. After some manipulation we can devise 
 approximations that look like integration methods in which the coefficients 
 are functions of T/H where H is the stepsize used. As this ratio tends to 
 zero, they tend to standard integration techniques. (Note that equation 
 (2. A) tends to a the N differential equation z' = g(t,z) as T tends to zero.) 
 For example, the multistep methods take the form 
 
- 11 - 
 
 k 
 
 (2.6) I (a z + H0 g ) = 
 j_n n— 1 1 n— i 
 
 If this is exact for polynomials z of degree k in t, the coefficients a 
 and 3 are polynomials of degree at most k in T/H. Note that each step of 
 this "multistep method" requires the evaluation of g(z), and this requires 
 the integration of the original differential equation through one cycle. 
 
 If the period T is unknown, or if it varies slowly, it can be 
 estimated by minimizing 
 
 t+T 
 / ll(y(s+T) - y(s)irds 
 t 
 
 with respect to T. In fact, this serves to define T(t) since it is not 
 
 otherwise defined. The previous formulas were stated for a fixed period T. 
 
 The same formulas can be used for a variable period by using a change of 
 
 independent variable to one in which the period is constant . When the 
 
 equations are expressed in the new independent variable, we get a 
 
 relationship similar to eq. (2.4) but with a fixed period. In addition, we 
 
 must add an equation for the real time. It takes the form 
 
 t(x + u) = t(x) + T(t) 
 
 where jj is the fixed period in the new independent variable t. 
 
 The general form of the method consists of a pair of integrators, the 
 inner integrator called from a program which determines g and estimates T, 
 plus the outer integrator which computes z. Each step of the outer 
 integrator requires one or more evaluations of g in a typical predictor- 
 
- 12 - 
 
 corrector operation. In fact, the outer integration can be done by minor 
 modifications of a conventional integrator, and the step and order 
 selection algorithms can be carried over. Each time the the function g is 
 evaluated for the outer integration, an integration over a full cycle by a 
 conventional integrator is performed. If the period can vary, it is 
 necessary to integrate over two full cycles to compute the term which has 
 to be minimized. 
 
 Non-autonomous problems in which the t-dependency of the differential 
 equation contains frequencies which determine the period can be handled by 
 the same techniques, but the large stepsizes H must be synchronized with 
 the period. If the t-dependencies are very slow compared to the period, 
 this synchronization is not important. 
 
 The method is not efficient if the problem also happens to be stiff 
 (meaning that there are rapidly decaying components present also) because 
 it is then necessary to use an implicit outer integration scheme which 
 means that a Jacobian must also be calculated. In this case, the Jacobian 
 is a state transition matrix over one period, and this matrix is never 
 sparse even if the system is sparse. The method is also not effective if 
 the solution is composed of several almost periodic functions with nearby 
 periods. These types of problems remain almost intractable. 
 
 2i'J±. Discontinuous problems 
 
 After each discontinuity, earlier information cannot be used to 
 extrapolate forward so either a single-step method must be used or a 
 restart of a multistep method is necessary. Most current automatic codes 
 for multistep methods use as a first step a first- or second-order one-step 
 
- 13 - 
 
 method that is also one of their multistep methods, and then use the 
 variable-order mechanism to build up the order. This is moderately 
 efficient, and, if it occurs infrequently is an adequate technique. 
 However, frequent restarting is costly. Therefore, we have developed a 
 Runge-Kutta starter for multistep methods that uses a single Runge-Kutta 
 like step to generate all the information needed to switch to a moderately 
 high order multistep method in the second step [16]. The method consists 
 of a single step of the form 
 
 k = hf(y + I 3k) , 1-1, ..., q, 
 
 y m = y n + l Y.k , m = 1, ..., s. 
 
 This is a q-stage R-K method that attempts to generate s new values in one 
 step. These s values can be used to start a multistep method of order s+1 
 or s+2. It has been shown that for implicit RK starters, it is possible to 
 generate a set of s points of order s using s stages, and for the more 
 promising explicit methods, we can generate a set of fourth-order points 
 using 6 stages. Hence, a non-stiff method can be started at fourth order 
 using six evaluations. In fact, there is a nine-parameter family of such 
 methods. Preliminary tests indicate that it reduces the number of function 
 evaluations to get started by about half. 
 
 A problem that still remains to be solved is that of detecting the 
 presence of a discontinuity efficiently. If nothing is known about its 
 nature, a code has great difficulty reducing the step until it lands right 
 on the discontinuity with a mesh point. 
 
- 14 - 
 
 3^, Problems presenting major difficulties 
 
 This section examines some of the problems which appear to be more 
 difficult to solve. Perhaps the most important of these is the large, 
 stiff problem, both because it is a practical problem, and its solution 
 will have an important impact on the application of ODE methods to PDE 
 problems via the method of lines. The major remaining problems with 
 variable-step, variable-order codes arise from the desire to put them on a 
 firm theoretical basis. In general, they work very well and further 
 advances are unlikely to improve their performance significantly. The last 
 two problems, multirate methods and low accuracy methods, are in their 
 infancy. The latter is very important in real-time simulation in which the 
 integration is time-critical , that is, it is of no interest unless one 
 second of simulated time can be computed in no more than one second of real 
 time! This problem may also be important in the method of lines for PDEs 
 where typical accuracies (two digits) are very low compared to typical 
 accuracies requested in ODE solution (5 to 15 digits) . 
 
 _3._1. Large stiff problems 
 
 We can distinguish two categories of "largeness" in stiff problems: 
 those in which an LU decomposition of the Jacobian is possible but a time 
 consuming part of the calculation, and those in which it is impossible to 
 contemplate factoring the Jacobian because of space considerations. In the 
 
 latter case it may be impossible to store the Jacobian explicitly, even if 
 
 3 
 it is sparse. The former consist of order of 10 equations, while the 
 
 latter may have as high as 10 equations or larger. The latter usually 
 
 arise from two- and three-dimensional PDEs treated by the method of lines. 
 
- 15 - 
 
 If the number of large eigenvalues causing stiffness is small, it may 
 be possible to treat the problem efficiently using ideas such as those in 
 the Correction in the Dominant Subspace (CDS) method described in Alfeld 
 and Lambert [1]. The basic idea behind this class of methods can be 
 explained as follows. If a stiff method, such as backward Euler, is 
 applied to the problem 
 
 y' = Ay 
 we get a difference equation of the form 
 
 Vl " R(hA)y n 
 
 where R(hA) = [I - hA] . The effect of large eigenvalues in hA is to 
 
 remove the corresponding eigenvectors in y ,.. Small eigenvalues cause the 
 
 n+l 
 
 corresponding eigenvectors to be multiplied by an order one approximation 
 
 hA 
 to e . If an explicit method such as Euler is used, we find that R(hA) = 
 
 [I + hA]. This correctly treats the eigenvectors corresponding to the small 
 
 eigenvalues, but amplifies those corresponding to large eigenvalues. The 
 
 CDS method consists of applying a simple (probably explicit) method, and 
 
 then damping the eigenvectors corresponding to the large eigenvalues (the 
 
 dominant subspace). This requires a knowledge of the dominant subspace 
 
 but, if its dimension is small, does not involve much storage or 
 
 calculation. The damping is applied by removing from the solution most of 
 
 the component in the dominant subspace. Details can be found in the cited 
 
 paper. The method requires a few inner products and vector operations so 
 
 that the calculation is proportional to the dimension of the system and the 
 
 dimension of the dominant subspace; much better than the n or n^ 
 
 operations required with a full LU handling of the Jacobian. The potential 
 
- 16 - 
 
 difficulties with this sort of approach are that a time varying Jacobian 
 means a relatively costly recalculation of the dominant supspace, and that 
 the vectors spanning a dominant subspace are not sparse even when the 
 matrix is very sparse (except in unusual circumstances) . 
 
 If a conventional direct method is applied to problems in the 
 
 intermediate size range, the linear equation subproblem has two important 
 
 3 
 components, the LU decomposition which may be as slow as n operations, and 
 
 2 
 the back substitution which may be as slow as n operations. A new 
 
 decomposition must be performed under a number of conditions. These are 
 
 when the actual Jacobian changes and when the order or stepsize changes. 
 
 The latter occurs because the matrix which must be decomposed has the form 
 
 [I - h3J] where J is the system Jacobian. Enright [12] suggested one way 
 
 of avoiding a complete refactorization when only h or 3 changes (the most 
 
 common case) . It involves factoring the matrix into SHS where S is a 
 
 suitable similarity or unitary transformation, and H is an easy matrix to 
 
 work with (for example, an upper Hessenberg matrix). Changes in h and 3 
 
 can then be accomodated by adding a multiple of the identity to H. Enright 
 
 and Kamel [13] have since suggested another scheme in which S is a unitary 
 
 transformation (a product of Householder transformations for example), and 
 
 a permutation matrix is also folded into S so that all of the large 
 
 elements of H appear in the first rows. The reduction to this form can be 
 
 stopped at an intermediate point when H consists of large elements in the 
 
 upper rows, zeros in the bottom left-hand part, and small elements in the 
 
 bottom right hand part. At this point, the small elements can be discarded 
 
 and the resulting matrix used as an approximation in the iteration. It 
 
 achieves an effect similar to correction in the dominant subspace since it 
 
- 17 - 
 
 damps the dominant space by applying an implicit method with a quasi-Newton 
 iteration to them. (Of course, the underlying method in this scheme must 
 have appropriate stability properties for stiff problems.) 
 
 In the case of extremely large systems, it appears that the only way 
 to handle the problem is by iterative methods so that at most the original 
 Jacobian has to be stored. (If the equations are derived by the method of 
 lines, the Jacobian elements often can be generated as needed, so even they 
 need not be stored explicitly.) Examination of the ideas used in the CDS 
 method suggests that the important property of any iterative method is that 
 it must be reasonably accurate for components in the dominant subspace. 
 Warming and Beam [26] explore approximate factorizations by operator 
 splittings that are related to ADI methods. In these, the equivalent of 
 the matrix [I - h(3J] is approximated by the product [I - h3J,] [I - h$J_] 
 where J = J. + J„. The matrices J. and J„ contain the components due to 
 derivatives in the x and y spatial directions, respectively. For 
 rectangular regions, the eigenvalues of J are the sums of eigenvalues of J. 
 and J-. Further, the eigenvectors corresponding to large eigenvalues of 
 either J. or J„ are in the null space of the other. As long as this 
 property is approximately true, the factorization is a reasonably accurate 
 approximation. Since the factors correspond to one-dimensional operators, 
 they are effectively tridiagonal for second-order differences, and suitably 
 
 banded for higher-order differences. Unfortunately, terms involving 
 
 2 
 3 /9x3y do not permit this factorization, but frequently these terms do not 
 
 contribute to the stiffness which arises from parabolic terms. 
 
 Consequently, they can be ignored in the quasi-Newton step. 
 
 The underlying theme of all of these methods is similar, and can be 
 
- 18 - 
 
 related to the use of one or more predictor-corrector schemes where the 
 corrector iteration is performed using quasi-Newton in which the Jacobian 
 is approximated by another matrix chosen to get a "simpler" factorization. 
 The important property of the approximation is that it must lead to fast 
 convergence, and this is determined by its ability to represent the 
 Jacobian reasonably accurately in the dominant subspace. 
 
 _3._2. Variable-step , variable-order methods 
 
 An adequate theory for automatic methods is still lacking. A theory 
 is needed to justify the error estimators used to control the step and 
 order. Most current estimators are based on the asymptotic analysis of 
 fixed-step, fixed-order methods in which a power series in the stepsize h 
 is developed for the error. Higher-order terms are dropped and differences 
 of the solution are used to estimates the derivatives that occur in the 
 leading terms. A variable-order, variable-step method is controlled by a 
 user parameter, say, £, which specifies the error tolerance. If the global 
 error and local error were a smooth function of £, one could consider an 
 asymptotic expansion in £ instead of h, but even if a code attempts to 
 adjust parameters so that the local error is exactly equal to the desired 
 multiple of £ at each step, the global error would still not be a 
 reasonable function of £ because the switch from one order to another is a 
 discrete process. In practice, codes do not achieve the desired local 
 error exactly, but accept any step that appears to give an error no larger 
 than desired. Furthermore, they tend to avoid too many changes to stepsize 
 and order in the interest of efficiency, and this tends to make the global 
 error an erratic function of £. As £ is reduced, the global error 
 generated by the code should reduce. Ideally, the reduction should be 
 
- 19 - 
 
 proportional, but this seems to be difficult if not impossible to achieve. 
 See Stetter [24]. An alternative approach could be based on the 
 effectiveness theory of Hull [17]. That is based on showing that a code 
 satisfies the error requirements for any problem from a specified class of 
 problems. However, this is not a trivial task, and it is not usually 
 possible to verify that the problem lies within the given class. It is 
 hard to see any promising directions for this problem. 
 
 _3._3. Multirate methods 
 
 A multirate method is a technique to handle each variable with a 
 stepsize appropriate to its behavior. Thus, if the system of equations is 
 
 dy i 
 
 the i-th variable should be integrated using a stepsize function h (t). 
 Each f must be capable of being evaluated separately because the t values 
 at which one f must be evaluated may be different from the t values used 
 for another. Since the f.'s depend on the y for other values of j, it is 
 necessary to interpolate for the values of y at the values of t used as 
 mesh points for y. but not for y . The cost of interpolation is non- 
 negligible compared to the cost of an integration step (in fact, it is 
 about the cost of a predictor step) . Hence, multirate methods can only be 
 competitive if one of two conditions is true: either the cost of function 
 evaluations must be very high or the dependency of one function f on other 
 values of y must be sparse. Then there can be considerable savings. In 
 the former case the savings can arise because the total number of f's 
 evaluated will be reduced; in the latter case the savings can arise because 
 
- 20 - 
 
 very few interpolations are needed. The latter case occurs when the ODE's 
 arise from a method of lines treatment of PDE's. 
 
 Even if the necessary conditions are met, there are many problems to 
 be solved before the method is practical. The first concerns the errors 
 produced by the interpolation process. While it is simple to show that a 
 straightforward scheme has the appropriate order of convergence (see Gear 
 [14]), it is necessary to derive estimation formulas for the errors created 
 by each integration process so that the stepsizes can be controlled 
 automatically. Perhaps the usual estimates of the derivatives will 
 suffice, but there is a danger that interpolation errors from other 
 components may interfere because these errors will have varying signs even 
 over the space of a few steps since they depend on the relationship between 
 the mesh used in one equation and that used in another. 
 
 When we have error estimation formulas, we have to determine how to 
 structure the software to overcome another serious problem, that of backup. 
 Orailoglu [20], in his masters thesis, created an experimental 
 implementation of a multirate method using a separate integrator for each 
 equation. Each integrator was actually a form of DIFSUB. An "event 
 driven" approach was used. The basic idea is to examine the estimated 
 stepsizes for all equations, and to choose that equation (and corresponding 
 variable) for integration which will progress to the earliest point in 
 time. The variable is then integrated, and the estimated stepsize for the 
 next step in that variable calculated. Each mesh point in each equation 
 can be thought of as an "event" from a simulation viewpoint. The 
 simulation proceeds by advancing to the next event to occur. The problem 
 with this approach is that the stepsizes for subsequent steps can only be 
 
- 21 - 
 
 estimated. If when the step is actually taken it turns out to be too 
 large, it must be reduced. The difficulty arises when one variable is 
 being integrated with a very small stepsize and another with a very large 
 one that turns out too be much too large. For example, suppose that y and 
 z are being integrated with stepsizes h and lOOh repectively. If they are 
 both currently known at time 0, 99 or 100 integrations of y will be 
 performed before one integration of z is attempted. When we attempt to 
 integrate z and find that the estimated stepsize is much too large, we will 
 reduce the stepsize and find ourselves trying to integrate z to a much 
 earlier point in time than the current t value of y. Consequently, we will 
 not have any values of y on which to base an interpolation for the 
 evaluation of the derivative of z. This problem may not be as serious as 
 the fact that y has been integrated using extrapolated values of z over an 
 interval that is too large for an integration step. Since the error bound 
 for interpolation or extrapolation has the same general form as the error 
 bound for integration, there is a danger that the integration of y will 
 have been contaminated with interpolation errors from z. Orailoglu tried 
 various techniques to avoid this problem. One method was to use "pre- 
 integration" where every step was done twice, with the first step used to 
 get a better estimate of the stepsize. Even though this was expensive, the 
 total number of function evaluations was still less than in a conventional 
 method. (Here, a function evaluation is an evaluation of a single f .) 
 Another mechanism was tried with some success . We monitored the coupling 
 of errors from variable to variable by examining the Jacobian of the 
 system. Unfortunately, this approach required knowledge of the Jacobian 
 even for non-stiff systems. We were able to "predict" when an estimated 
 stepsize had to be reduced by this technique, but the cost was high. We 
 
- 22 - 
 
 have started to investigate a different approach based on the observation 
 that the coupling from a rapidly varying component to a slowly varying 
 component must be weak (or the latter would not be slowly varying) . 
 Therefore, it may be possible to to integrate the slowly changing 
 components over their large stepsizes first, their values at the end of the 
 interval being relatively independent of the fast components. Then, the 
 fast components can be integrated, using more accurate interpolation for 
 the intermediate values of the slow components rather than extrapolation. 
 The process looks like a top-down, successive refinement process because it 
 can be implemented in the following way: 
 
 (i) Initially, assume all variables are known at the same time value, 
 say t = 0. 
 
 (ii) Select the set of variables known at the earliest time. (Initially, 
 this is the complete set.) Integrate these equations with the 
 largest of the estimated stepsizes of any of the variables. 
 
 (iii) Check the error of each variable individually. 
 
 (iv) Reject the steps on the variables with large errors, and halve their 
 suggested stepsizes. 
 
 (v) Repeat the whole process from step (ii) • 
 
 This technique performs a compound step whose length is the largest 
 stepsize for which one equation can be integrated successfully. It does 
 successive halving of other stepsizes, and reintegrates them until the 
 stepsize is small enough. Although this may lead to many wasted steps, the 
 fact that the binary subdivision reaches the stepsize in a logarithmic 
 
- 23 - 
 
 number of subdivisions means that the number of failed steps never exceeds 
 the number of successful steps. While 50% efficiency may seem poor, the 
 name of the game is to reduce the number of integrations in the slow 
 components to such an extent that we are ahead. When applied to equations 
 arising from the method of lines, this approach automatically refines the 
 mesh in regions of rapid transition. 
 
 There are, however, many unanswered questions with this technique. 
 For example, the choice of the outermost large step is critical, so a 
 theory to guide its selection is needed. Error estimates are also 
 critical. If an estimate is small for a particular variable, its step will 
 not be subdivided further, even though subdivision of steps for other 
 variables may cause so much change that the integration of the first 
 variable was invalid. 
 
 ^.^_. Low accuracy methods 
 
 Current ODE methods are based on the concept of polynomial order. 
 This matches the numerical solution to the true solution in an asymptotic 
 sense. Unfortunately, for large stepsizes several problems arise. The 
 first is that error estimates are unreliable because they also are based on 
 asymptotic theory. The second is that the actual methods are inefficient 
 and display poor stability properties. Low accuracy methods are needed in 
 a number of applications, in particular, real-time simulation [15] and the 
 time integration of PDEs. In real-time integration, one of the most 
 important properties is the stability correspondence between the actual 
 system and the numerical model because the simulation is frequently being 
 done to test the "handling properties" of a proposed device. The method 
 
- 24 - 
 
 should be stable or unstable as the actual device Is. The same 
 characteristics are probably desirable for applications to PDEs . Since we 
 are talking about a large stepsize, polynomial accuracy is probably not 
 important. Some form of minimax approximation is more appropriate (as it 
 is in approximation over an interval rather than in the neighborhood of a 
 point). For example, if we consider the class of differential equations 
 
 y' = x y 
 
 where X is known to lie in some region D of the complex plane, we can find 
 the Runge-Kutta or multistep method with a given number of stages or steps 
 that gives the best approximation to the class of equations in the minimax 
 sense. The approximation also depends on the stepsize h, so the actual 
 maximum error depends on the stepsize. In fact, it can be seen that this 
 dependency is monotonic, so the step can be chosen when the region D is 
 given and the error tolerance is specified. The general approach would be 
 to assume a form for the class of problems with some parameters. Formally, 
 we would have a class of problems C(p) where the parameters p are assumed 
 to lie in a region D of the parameter space. A class of methods M(h,a) 
 would be considered, where a is a set of method coefficients, and h is the 
 stepsize. For a given stepsize h, the coefficients should be picked to 
 minimize the worst case error of the class of problems. Thus, if the error 
 for problem C(p) using method M(h,a) is E(p,h,a), we would pick a for a 
 given h to minimize the worst case error to get 
 
 E(h) = min max E(p,h,a) 
 a p in D 
 
 The stepsize would be picked so that E(h) is sufficiently small, that is, 
 
- 25 - 
 
 so that E(h) < e. An automatic program might function by estimating the 
 regions in which the parameters lie, and using that information to reduce 
 the region over which the minimax approximation is taken. The utility of 
 an approach of this form is dependent on determining some simple 
 approximations to the minimax solutions. Clearly one cannot hope to solve 
 a complex minimization problem at each integration step! 
 
 4^. Conclusion 
 
 This paper has surveyed a number of recent advances, many of which are 
 very important, either offering major advances in the theory or suggesting 
 ways of handling currently very difficult problems. It has also examined 
 some problems that remain very challenging, and these are problems whose 
 solution can have a significant impact on the solution of PDE problems. 
 Some of these problems have appeared very recently, which suggests that 
 there will continue to be new problems to solve in this field for some 
 time, and that continued effort can reap worthwhile rewards. 
 
 BIBLIOGRAPHY 
 
 [1] Alfeld, P., Lambert, J. D. , Correction in the Dominant Space: A 
 Numerical Technique for a Certain Class of Stiff Initial Value 
 Problems, Math. Comp . , _31, # 140, pp 922-938, October, 1977. 
 
 [2] Bickart, T., An Efficient Solution Process for Implicit Runge Kutta 
 Methods, SIAM J. of Num. Anal., _14, pp 1022 - 1027, Dec 1977. 
 
 [3] Burrage, K. , High Order Algebraically Stable Punge Kutta Methods, BIT 
 J_8, pp 373 - 383, 1978. 
 
- 26 - 
 
 [4] Burrage, K., and Butcher, J., Stability Criteria for Implicit Runge 
 
 Kutta Methods, SIAM J. of Num. Anal., _16_, #1, pp 30 - 45, 1979. SIAM 
 1979 
 
 [5] Butcher, J. , Some Implementation Schemes for Implicit Runge Kutta 
 
 Methods, to appear, Proceedings of the 1979 Dundee Biennial Numerical 
 Analysis Conference, Springer Verlag, Berlin. 
 
 [6] Butcher, J., On the Implementation of Implicit Runge Kutta Methods, 
 BIT _16, pp 237 - 240, 1976. 
 
 [7] Curtiss, C. F., Hirschfelder, J. 0., Integration of Stiff Equations, 
 Proc. Nat. Acad. Science, IT. S., _38, pp 235 - 243, 1952. 
 
 [8] Dahlquist, G. , Numerical Integration of Ordinary Differential 
 Equations, Math. Scandinavica, _4, pp 33 - 50, 1956. 
 
 [9] Dahlquist, G. , A Special Stability Problem for Linear Multistep 
 Methods, BIT, 3, pp 27 - 43, 1963. 
 
 [10] Dahlquist, G. , G-Stability is Equivalent to A-Stability, BIT, _18, pp 
 384 - 401, 1978. 
 
 [11] Dahlquist, G. , Some Properties of Linear Multistep and One-leg Methods 
 for Ordinary Differential Equations, in Working Papers for the 1979 
 SIGNUM Meeting on Numerical ODEs , ed. R. Skeel, Dept. of Computer 
 Science Report # 963 University of Illinois at Urbana-Champaign. pp 
 1-1 to 1-4. 
 
 [12] Enright, W. H. , Improving the Efficiency of Matrix Operations in the 
 Numerical Solution of Stiff Ordinary Differential Equations, ACM 
 
- 27 - 
 
 Trans, on Math. Software, 4^ #2, pp 127 - 136, June 1978. 
 
 [13] Enright, W. H., Kamel, M. S., Automatic Partitioning of Stiff Systems 
 and Exploiting the Resulting Structure, in Working Papers for the 1979 
 SIGNUM Meeting on Numerical ODEs , ed. R. Skeel, Dept. of Computer 
 Science Report # 963 University of Illinois at Urbana-Champaign, pp 
 12-1 to 12-4. 
 
 [14] Gear, C. W. , Multirate Methods for Ordinary Differential Equations, 
 
 File # 880, Department of Computer Science, University of Illinois at 
 Urbana-Champaign, 1974. 
 
 [15] Gear, C. W., Simulation: Conflicts between Real-Time and Software, in 
 Mathematical Software III, ed J. R. Rice, Academic Press, New York, 
 1977. 
 
 [16] Gear, C. W. , Runge-Kutta Starters for Multistep Methods, Dept. of 
 computer Science Report # 938, University of Illinois at Urbana- 
 Champaign, Sept 1978. 
 
 [17] Hull, T. E., The Effectiveness of Numerical Methods for Ordinary 
 
 Differential Equations, in Studies in Numerical Analysis 2, ed. J. M. 
 Ortega and W. C. Reinboldt, SIAM, Philadelphia, 1970, 
 
 [18] Miranker, W. L., Wahba, G., An Averaging Method for the Stiff Highly 
 Oscillatory Problem, Mathematics of Computation, J30, pp 383 - 399, 
 1976. 
 
 [19] Nevalinna, 0., Liniger, W., Contractive Methods for Stiff Differential 
 Equations, BIT, J_8, pp 457 - 474, 1978. 
 
- 28 - 
 
 [20] Orailoglu, A., A Multirate Ordinary Differential Equation Integrator, 
 MS Thesis, Department of Computer Science Report # 959, University of 
 Illinois at Urbana-Champaign, 1979. 
 
 [21] Petzold, L., An Efficient Method for Highly Oscillatory Ordinary 
 
 Differential Equations, Ph. D. Thesis, Dept of Computer Science Report 
 # 933, University of Illinois at Urbana-Champaign, August 1978. 
 
 [22] Petzold, L., Gear, C. W., Methods for Oscillating Problems, Dept. of 
 Computer Science File # 889, University of Illinois at Urbana- 
 Champaign, October, 197 7. 
 
 [23] Skeel, R. D. , Kong, A. K-Y. , Blended Linear Multistep Methods, ACM 
 Transactions on Mathematical Software, ^, pp 326 - 345, 1977. 
 
 [24] Stetter, H. J., Tolerance Proportionality in ODE Codes, in Working 
 Papers for the 1979 SIGNUM Meeting on Numerical ODEs, ed. R. Skeel, 
 Dept. of Computer Science Report # 963 University of Illinois at 
 Urbana-Champaign, pp 10-1 to 10-6. 
 
 [25] Wanner, G. , Hairer, E., and Norsett, S. P., Order Stars and Stability 
 Theorems, BIT, J_8, pp 475 - 489, 1978. 
 
 [26] Warming, B. F. , Beam, R. M. , Factored, A-Stable, Linear Multistep 
 Methods - An Alternative to the Method of Lines for Mult id imens ions , 
 in Working Papers for the 1979 SIGNUM Meeting on Numerical ODEs, ed. 
 R. Skeel, Dept. of Computer Science Report # 963 University of 
 Illinois at Urbana-Champaign, pp 9-1 to 9-3. 
 
cu 
 
 CO 
 
 « 
 o 
 u 
 
 •H 
 n) 
 ■U 
 W 
 
 ■u 
 C 
 ai 
 
 e 
 a, 
 o 
 
 H 
 CU 
 > 
 
 cu 
 Q 
 
 cu 
 u 
 
 3 
 M 
 ■H 
 Pn 
 
(O 
 
 UJ 
 
 -I 
 
 CO 
 O 
 CC 
 
 Q. 
 
 o_ 
 ui 
 
 h- 
 
 CO 
 
 i 
 
 UJ 
 
 _J 
 m 
 < 
 
 CC Q 
 < Z 
 > < 
 
 CO 
 O 
 O 
 
 I 
 h- 
 UJ 
 
 tr 
 
 UJ 
 Q 
 
 cc 
 o 
 
 X 
 
 * Hi 
 
 01 
 
 5 
 
 UJ 
 
 _l 
 
 CD 
 
 o 
 
 CC 
 Q. 
 
 CO 
 
 UJ 
 
 _J 
 < 
 > 
 
 Z 
 UJ 
 
 o 
 
 UJ 
 
 > 
 
 < 
 
 z 
 
 o 
 < 
 
 5 
 
 01 
 
 
 2 
 
 
 UJ 
 
 
 _J 
 
 
 m 
 
 
 o 
 
 
 CC 
 
 CO 
 
 CL 
 
 2 
 
 
 UJ UJ 
 
 o 
 
 O H 
 
 z 
 
 cc en 
 
 < 
 
 < >■ 
 
 _J 
 _l 
 
 O 
 
 C/1 
 
 ERY 
 TIFF 
 
 o 
 
 > to 
 
 CO 
 
 2 
 
 UJ 
 
 _J 
 m 
 o 
 tr 
 a. 
 
 CO 
 O 
 
 Z 
 
 Z 
 
 o 
 a 
 
 CO 
 
 CO 
 
 a 
 O 
 
 x 
 »- 
 
 UJ 
 UJ 
 
 I- 
 < 
 
 CC 
 
 I- 
 
 _l 
 
 D 
 
 2 
 
 CO 
 
 a 
 o 
 
 x 
 
 H 
 UJ 
 
 S 
 
 >- 
 o 
 < 
 
 CC 
 
 o 
 o 
 
 < 
 
 o 
 
 o 
 
 CO 
 
 en 
 
 o 
 
 CO 
 
 w 
 
 Q 
 O 
 
 en 
 
 B 
 
 rH 
 O 
 
 cu 
 •u 
 CO 
 
 o 
 n 
 
 a 
 
 CN 
 
 <U 
 M 
 
 3 
 
 00 
 •H 
 
 O 
 cx» 
 
c 
 o 
 
 •H 
 
 ■u 
 
 iH 
 O 
 
 1-1 
 O 
 
 •H 
 O 
 CO 
 
 o 
 
 m 
 
 3 
 •H 
 
M 
 
 P 
 
 + 
 
 in 
 
 ^^^ \ 
 
 ^^ / 
 
 / 
 
 2,' 
 
 > \ 
 \ 
 
 ■> -^ -^ 
 
 
 
 </> 
 
 01 
 
 ex 
 o 
 
 rH 
 
 > 
 
 c 
 w 
 
 CO 
 
 OJ 
 
 >-i 
 
 3 
 60 
 •H 
 
BIBLIOGRAPHIC DATA 
 SHEET 
 
 1. Report No. 
 
 UIUCDCS-R-79-983 
 
 3. Recipient's Accession No. 
 
 4. Title and Subtitle 
 
 ODEs: Is there anything left to do? 
 
 5. Report Date 
 
 October 1979 
 
 6. 
 
 7. Author(s) 
 
 C. W. Gear 
 
 8. Performing Organization Rept. 
 
 No. 
 
 R-79-983 
 
 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. 
 
 US ENERGY/ 
 
 EY-76-S-02-2383 
 
 12. Sponsoring Organization Name and Address 
 
 U. S. Department of Energy 
 Chicago Operations Office 
 9800 South Cass Avenue 
 
 13. Type of Report & Period 
 Covered 
 
 presented at SIAM 
 Spring Mtg., June 1979 
 14. 
 
 Argonna, IL 6Q/i39 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 Two decades ago many people thought that there was very little left to occupy 
 the researcher in the numerical initial value problem for ordinary differential 
 equations but there has since been an enormous increase in the effort devoted 
 to this problem and related software. Is everything finally done so we can 
 turn to more difficult areas? This paper examines the development of numerical 
 solution techniques from the identification of a problem to the never-final 
 preparation of automatic codes for the solution of classes of similar problems. 
 It discusses a number of ODE problems that are only a part of the way along 
 this path of development. These problems include ones with highly oscillatory 
 solutions, ones with frequent discontinuities, and ones in which very little 
 accuracy is needed (and which should be integrated at low cost). 
 
 17. Key Words and Document Analysis. 17a. Descriptors 
 
 ordinary differential equations 
 numerical software 
 
 17b. Identifiers /Open-Ended Terms 
 
 17e. COSATI Field/Group 
 
 18. Availability Statement 
 
 unlimited 
 
 FORM NTIS-35 ( 10-70) 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIFIED 
 
 21. No. of Pages 
 
 32 
 
 22. Price 
 
 USCOMM-DC 40329-P7 1 
 
Form AEC-427 
 
 (6/68) 
 
 AECM 3201 
 
 U.S. ATOMIC ENERGY COMMISSION 
 
 UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT 
 
 ( See Instructions on Reverie Side ) 
 
 1. AEC REPORT NO. 
 
 COO-2383-0061 
 
 2. TITLE 
 
 ODEs: Is there anything left to do? 
 
 3. TYPE OF DOCUMENT (Check one): 
 
 r%] a. Scientific and technical report 
 
 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): 
 
 [">3 a. AEC's normal announcement and distribution procedures may be followed. 
 ~2 b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 ] 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 U-C 
 
 Urbana, IL 61801 
 
 ^ 
 
 ^>S i *-- 
 
 Date 
 
 October 1979 
 
 y FOR AEC USE ONLY 
 
 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 8. PATENT CLEARANCE: 
 
 Q a- AEC patent clearance has been granted by responsible AEC patent group. 
 I~) b. Report has been sent to responsible AEC patent group for clearance. 
 I I c. Patent clearance not required. 
 
JUN l 2 1980 
 
FEB 2