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.