[SciPy-user] odeint
Ryan Gutenkunst
rng7 at cornell.edu
Tue Oct 4 11:00:49 CDT 2005
LOPEZ GARCIA DE LOMANA, ADRIAN wrote:
> Hi people,
>
> I'm using odeint from integrate to solve a system of ODEs. Although the output seems correct
> a strange error pop up at the terminal:
>
> [alopez at thymus scipy]$ python single_global.integ.py
> Traceback (most recent call last):
> File "single_global.integ.py", line 23, in func
> xdot[7] = + (v_NACT * x[5]**m_NACT) / (x[5]**m_NACT + k_NACT**m_NACT) - k_deg_h * x[7]
> ValueError: negative number cannot be raised to a fractional power
> odepack.error: Error occured while calling the Python function named func
<snip>
> And the output file, looks OK, there is no single negative number.
>
> Can anyone reproduce the error? Can anyone tell me what's going on?
I get the same error running your code, although it stops the
integration in my case.
I ran into the same error running some of our own biochemical
simulations. My hunch (someone correct me if this is nonsense) is that
this error comes from the extrapolating steps the integrator takes.
odeint is a variable step size integrator, so the routine tries a step,
checks tolerances, adjusts step as necessary, and so on. If the value of
your variable is getting close to zero, the extrapolated step may take
it to a negative value. When odeint evaluates xdot at that point, you
see those errors.
The solution we came up with was to integrate in terms of the logarithm
of the variable values. This ensures that they are always positive, and
for our systems the speed penalty is pretty small.
It's pretty easy to do since:
d log(x)/dt = dx/dt * 1/x
So you might want to try something like:
def func_log(log_x, t, *args):
x = exp(log_x)
dx_dt = func(x, t, *args)
return dx_dt/x
Of course, your initial condition must be log(x_IC). And a zero IC will
kill you, but maybe picking something like 1e-16 is reasonable for your
problem.
Cheers,
Ryan
> Thanks in advance,
>
> Adrián.
>
--
Ryan Gutenkunst |
Cornell LASSP | "It is not the mountain
| we conquer but ourselves."
Clark 535 / (607)227-7914 | -- Sir Edmund Hillary
AIM: JepettoRNG |
http://www.physics.cornell.edu/~rgutenkunst/
More information about the SciPy-user
mailing list