[SciPy-user] odeint

LOPEZ GARCIA DE LOMANA, ADRIAN alopez at imim.es
Tue Oct 4 12:52:59 CDT 2005


Thanks for your answer, I'll check if it works in my machine. I guess your're right, the errors appear when the values of the simulation are close to zero.

The problem is that I'm not just interested on solving this specific example but it is part of a bigger program where I search for the parameters, run the simulations, and evaluate a score given some data. I'll check how if I can integrate your suggestion into the program.

But, don't you think we are talking about a bug? Overall, if we want to compare scipy with Octave ... I mean, I definetly prefer SciPy, it's much faster, but I also want it to be reliable, and I don't want to be caring if the simulation is "too" close to zero.

Thanks for all, 

Adrián.

-----Original Message-----
From: scipy-user-bounces at scipy.net on behalf of Ryan Gutenkunst
Sent: Tue 04/10/2005 17:00
To: SciPy Users List
Subject: Re: [SciPy-user] odeint
 
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/

_______________________________________________
SciPy-user mailing list
SciPy-user at scipy.net
http://www.scipy.net/mailman/listinfo/scipy-user



More information about the SciPy-user mailing list