[SciPy-user] integrate.odeint , stiff chemical equations and mass conservation -any hint?
Tue Jun 2 13:29:06 CDT 2009
2009/6/2 ms <email@example.com>:
> I am trying to integrate a large (50-200 equations) system of chemical
> kinetics ODEs using scipy.integrate.odeint.
> The qualitative behaviour of the system looks sensible, but mass is not
> conserved at all -it increases a lot during time, sometimes wildly,
> depending on parameters.
> This puzzles me, and I am *very new* to this kind of problems; however
> reading here and there it seems that stiff ODEs can show this kind of
> artefacts. I checked the equation values and looks stiff;
> infodict['mused'] is always 1... so odeint sees it as stiff too. This
> does not look good.
> I wonder what can one do to tame this beast. Any hint?
In terms of software, you could take a look at the pydstool package,
which has several more modern C-level solvers (which may particularly
help with stiff systems), as well as various other useful tools for
dealing with dynamical systems.
In terms of your specific problem, you could (using the above tool)
impose mass conservation as an algebraic constraint. This should, I
think, help the solver slow down in situations where mass conservation
would otherwise be violated by the (necessarily approximate) solution.
On the other hand, it may be more valuable to keep the total mass as a
free parameter so that you can judge the quality of your solutions by
looking at how much it varies. After all, the total mass is only one
direction in which your approximate solution can deviate from the true
Which of these two approaches is more useful probably depends on your
problem. If it is one where, like some problems in planetary dynamics,
even if the solution is somewhat wrong in detail, it's still useful as
long as the energy (and angular momentum) don't change, then it's
useful to impose the constraint (or in the case of planetary dyamics,
use a solver that has that constraint built in). On the other hand, if
total mass does not have a special status beyond the fact you know it
shouldn't change, so that errors in other directions are just as
important as errors in mass, it's probably more useful to keep it as
an error-monitoring mechanism.
It's also, of course, possible there's a bug in your derivative
function - it's worth checking that the derivative vector is always
orthogonal to the gradient of mass-as-a-function-of-your-parameters.
But you are probably right, and this is probably just a tough ODE.
> SciPy-user mailing list
More information about the SciPy-user