[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