[SciPy-User] Odeint: not following input time

Jeroen Meidam j.meidam at gmail.com
Fri Jun 10 09:02:21 EDT 2011


Hello,

I have been searching for a solution on other forums, and this mailing 
list as well, but without luck. So I decided to post my problem.

When I run a simple odeint example, I get nice expected results. When I 
try my own odeint implementation I get strange results. What's happening 
in my own implementation of odeint is that the integrator seems to 
ignore my input time-line. I tested this by putting a print statement 
inside the function list.

Lets take the working example program (code1 at bottom of message):
I printed the time inside f(), to check what time it actually uses. It 
prints the expected 0 to 100 time sequence.

Now when I try this on my own code (code2, which I completely dumbed 
down to the most simple case imaginable, all derivatives are zero):
I again print the time, but now I get:
time: 0.0
time: 0.000293556243968
time: 0.000587112487935
time: 2.93614955216
time: 5.87171199184
time: 8.80727443152
time: 38.1628988283
time: 67.5185232251
time: 96.8741476218
time: 390.43039159
time: 683.986635557
time: 977.542879525
time: 3913.1053192
Which makes no sense at all to me.
First of all it exceeds my time limit and second, it only prints a few 
values.

Does anyone know what could cause this?


***** code1 ******
from scipy.integrate import odeint
from pylab import plot, axis, show

# Define the initial conditions for each of the four ODEs
inic = [1,0,0,1]

# Times to evaluate the ODEs. 800 times from 0 to 100 (inclusive).
t = linspace(0, 100, 800)

# The derivative function.
def f(z,time):
     """ Compute the derivate of 'z' at time 'time'.
         'z' is a list of four elements.
     """
     print time
     wx = sqrt(2.)
     wy = 1
     return [ z[2],
              z[3],
              -wx * z[0],
              -wy * z[1] ]

# Compute the ODE
res = odeint(f, inic, t)

# Plot the results
plot(res[:,0], res[:,1])
axis('equal')
show()
*****************

***** code2 ******
import numpy as np
from scipy.integrate import odeint
import matplotlib
from matplotlib.pyplot import figure, show, rc, grid
import pylab as pl

# Initial conditions
r0 = 0.0025**(-2./3) #omega0 is just a number
phi0 = 0.0
pphi0 = r0**(1./2)
pr = 0.0

T = 1200
N = 500

t = np.linspace(0,T,N)

p = [M,nu]
init_cond = [r0,phi0,pr0,pphi0]

def vectorfield(state,time,p):
     """
     Arguments:
         state :  vector of the state variables (reduced)
         t :  reduced time
         p :  vector of the parameters
     """
     print 'time:', time

     dr_dt    = 0.
     dphi_dt  = 0.
     dpr_dt   = 0.
     dpphi_dt = 0.

     return [ dr_dt,dphi_dt,dpr_dt,dpphi_dt ]

sol = odeint( vectorfield, init_cond, t, args=(p,) )

#followed by some plotting stuff, which is not important for now


******************



More information about the SciPy-User mailing list