From: glgesm@engr.psu.edu (Gary L. Gray) In article , horras@aries.scs.uiuc.edu (Greg Horras) wrote: > I am a graduate student in physical chemistry studying models of > chemical reactions in stirred tank reactors. These reactions are > typically modeled by systems of ordinary differential equations > which must be numerically integrated. Some of the published models > exhibit exotic, irregular behavior, including what has been diagnosed > as being chaos. > > I have attempted to reproduce some of these models using a GEAR routine > > but have not met with success (but I have many more to try). > > Often, crucial infomation is not provided in the publication such that > the models canot otn be reproduced. Also, numerical integration routines are very > subtle things, and many parameters must be set in order to get a specific > result. > > I was wondering; has anyone out > there come across a system of ODE's which gives chaotic trajectories > under some integration conditions but gives non-chaotic trajectories > under other conditions? If so, I would love to show this to some > people.....I would not be surprised if there is a system such as this > somewhere. > We are currently working with a system that we have modeled by four coupled and highly nonlinear ODEs. We use a routine due to Petzold and Hindmarsh called lsode that uses either a variable order, variable step-size Adams method for non-stiff systems or a Gear routine for stiff systems. We can do two runs using exactly the same physical parameters and initial conditions in our equations of motion and by just changing how often we output the results, we can get vastly different answers (lsode uses interpolation to output results at the user-specified time). Integrating exactly the same equations in a fixed step size Runge-Kutta routine gives yet a different set of answers. If the equations are nonlinear, it often doesn't take much for the type of behavior you describe to occur. From: mhosea@shavano.math.niu.edu (Michael Hosea) In article <1995Nov6.120505.15@jarvis.cs.toronto.edu>, Wayne Hayes wrote: >Michael Hosea wrote: >>I'm more familiar with VODE (Hindmarsh, Brown, and Byrne) than the >>older LSODE, but I suspect that the frequency of output requested is >>used to help determine the scale of the problem and hence has an effect >>on step size selection. This would explain the different results. > >Not just different, but "vastly different". It's the chaotic nature >of the problem to exponentially magnify small changes. Not necessarily, because the magnification of the perturbation may occur only over a small interval, after which the integrations may be stable again. Think of the Lorenz model. Be that as it may, I think you missed the point. Gary indicated that he got different results simply by requesting output more or less frequently. This was (quite naturally) surprising to him because LSODE uses interpolation for dense output, so one would expect that the underlying integration would be precisely the same when performed on the same computing equipment (using the same object code for the integrator), irrespective of how chaotic or otherwise unstable the system might be. I was explaining that LSODE does in fact (I've since checked) make use of the frequency of output in a way that would change the underlying integration. In particular, if you let LSODE compute an initial step size, it makes use of the distance from T0 to TOUT to get a rough estimate of the scale of the problem. Because a different initial TOUT generally changes the initial step size, a small perturbation is introduced which can explain the different results. The code SDRIV2 (most should use the double precision version DDRIV2 on workstations and PCs) does the same thing. Shampine has a nice example in his most recent book (example 3.1) where DDRIV2 gives very different results on a stable problem when the frequency of output requested was the only difference in the input data. It isn't a bug, it's a feature. ;-) >>RKSUITE, but in this case I doubt it would help. > >That probably won't help. If the global error estimate is anything >other than an exponentially growing number, I'd suspect the error >estimation routines. It won't make longer, more reliable integrations possible, but it can help the researcher. A global error estimate, like that in RKSUITE, can indicate when the integration results are no longer reliable. The fact that you can get vastly different results with LSODE, SDRIV2 and RKSUITE by varying the tolerance means that if you just run them one time you might not know that the results were unreliable. If you use RKSUITE with the global error estimate turned on, it will stop when it no longer "believes" in the global error estimate. Shampine (and I guess Brankin and Gladwell, though Shampine is the one who made the comment to me) had chaotic systems in mind when he designed it that way. Naturally, if you're prudent and you don't have a global error estimate available, you'll do a few runs with different tolerances (at least) and compare the results. >(b) it decides on it's own of the problem is stiff or non-stiff, so you >don't have to worry about it (too much). FYI, there is a variant of LSODE (LSODA) that detects stiffness and switches from AM to BDF methods automatically. Dynamical approach study of spurious steady-state numerical solutions of nonlinear differential equations. 1. The dynamics of time discretization and its implications for algorithm development in computational fluid dynamics. H.C. Yee, P.K. Sweby and D.F. Griffiths, Journal of Computational Physics, Vol. 97, p.249-310 (1991). Lorenz, E. (1989) ``Computational Chaos -- A Prelude to Computational Instability,'' Physica D, 35, pp. 299-317. -------------------------------------------------------------------- Title: On the reliability of gravitational N-body integrations Authors: QUINLAN, GERALD D.; TREMAINE, SCOTT Affiliation: AA(Canadian Inst. for Theoretical Astrophysics, Toronto, Canada) AB(California Inst. of Technology, Pasadena) Journal: Royal Astronomical Society, Monthly Notices (ISSN 0035-8711), vol. 259, no. 3, p. 505-518. Publication Date: 12/1992 Origin: STI Category: Astrophysics NASA/STI Keywords: CELESTIAL MECHANICS, MANY BODY PROBLEM, NUMERICAL INTEGRATION, STELLAR ORBITS, STELLAR SYSTEMS, STAR CLUSTERS, STELLAR GRAVITATION Bibliographic Code: 1992MNRAS.259..505Q Abstract In a self-gravitating system of point particles such as a spherical star cluster, small disturbances to an orbit grow exponentially on a time-scale comparable with the crossing time. The results of N-body integrations are therefore extremely sensitive to numerical errors: in practice it is almost impossible to follow orbits of individual particles accurately for more than a few crossing times. We demonstrate that numerical orbits in the gravitational N-body problem are often shadowed by true orbits for many crossing times. This result enhances our confidence in the use of N-body integrations to study the evolution of stellar systems. ---------------------------------------------------------------------- Title: Roundoff error in long-term planetary orbit integrations Authors: QUINN, THOMAS; TREMAINE, SCOTT Affiliation: AB(Toronto, University, Canada) Journal: Astronomical Journal (ISSN 0004-6256), vol. 99, March 1990, p. 1016-1023. Research supported by the Connaught Fund. Publication Date: 03/1990 Origin: STI Category: Astronomy NASA/STI Keywords: ORBITAL MECHANICS, PLANETARY ORBITS, SOLAR SYSTEM, CARTESIAN COORDINATES, EQUATIONS OF MOTION, ERROR ANALYSIS, FLOATING POINT ARITHMETIC, JUPITER (PLANET) Bibliographic Code: 1990AJ.....99.1016Q Abstract Possible sources of roundoff error in solar system integrations are studied. It is suggested than when floating point arithmetic is optimal, the majority of roundoff error arises from two sources: the approximate representation of the numerical coefficients used in multistep integration formulas and the additions required to evaluate these formulas. An algorithm to remove these two sources of error in computers with optimal arithmetic is presented. It is shown that the corrections result in a substantial reduction in the energy error at the cost of less than a factor of 2 increase in computing time.