[SciPy-User] Odeint Integration Method and its time step

[Athanasios Anastasiou]

Hello everyone

I am trying to run some simple examples with odeint in order to get a
better understanding of it but amongst other things i am a little bit
confused by the time parameter (t) that is passed to the function (f).

The specific questions that i have are:

1) What exactly does t represent in the derivatives function f? Is it
absolute time?
2) Which method is odeint using to integrate the derivatives function
and why does it call it more than one time for the same time instance?
(t)
3) Considering the following odeint scipy examples:
http://www.scipy.org/Cookbook/CoupledSpringMassSystem
and
http://www.scipy.org/LoktaVolterraTutorial
Why isn't time (or the time step) taken into account within the
function that returns the derivatives at each time instant?

For more information about these please see below.


Looking forward to hearing from you
Athanasios Anastasiou




Here are some more pieces of information regarding what i have done
and the results that generated those questions above. Please note that
i have not previously used the Fortran code that the scipy functions
are based on.

The first thing i tried was a spring-mass system with something like this:

def dS(y,t,K,M,dT,c):
  dY=[0.0,0.0]
  F = (-K*y[1])-(c*y[0]) #Forces acting on the spring... (linear
returning force, friction damping)
  a = F/M                #...cause an acceleration...
  dY[0] = a * dT         #...which sets the body in motion at some
velocity over the specific time step...
  dY[1] = y[0] * dT      #...which displaces the body by dS at every dT step
  return dY

Fs=512
dT = 1.0/Fs
y0 = [0.0,2.0]
tSim = 4 #This is 4 seconds...
nSim = tSim*Fs #...which translates to so many discrete samples

theTime=range(0,nSim)
S = SI.odeint(dS,y0,theTime, args=(32.0,4.0,dT,0.0))
pylab.plot(theTime,S[:,1]); #This plot is of course in "samples"
rather than "seconds"

This works as expected.

The next step was to "drive" the spring-mass system with a sinusoid like this:
(Only the dS function is given here the rest of the code from above
remains almost the same)
def dS(y,t,K,M,dT,c):
  dY=[0.0,0.0,0.0]
  F = (-K*y[1])-(c*y[0]) + scipy.sin(y[2]) * 8 #Forces acting on the
spring... (linear returning force, friction damping)
  a = F/M                #...cause an acceleration...
  dY[0] = a * dT         #...which sets the body in motion at some
velocity over the specific time step...
  dY[1] = y[0] * dT      #...which displaces the body at every dT step
  dY[2] = 2.0*scipy.pi*4.0*dT #A completely independent integration of
the phase component
  return dY

(In this case the spring mass system is excited with a sinusoidal
force running at a frequency of 4 Hz having an amplitude of 8 N)

This again, seems to be working as expected but the mean amplitude
does not have ONE resonant frequency. It seems to be "oscillating"
erratically with a change in the driving frequency. If i drive the
system with really low frequencies (0.5 Hz) the mean amplitude is
greater than when i am driving the system at its resonating frequency
of (1/2*pi)*sqrt(k/m) :-(

At this point i started suspecting that there is something i am
ignoring or not doing right in the whole simulation and decided to
examine the values of t more closely.

This is when i discovered that t(n),t(n+1) is not constant as i am
assuming with that dT over there. The first thing i did was to obtain
the dT from the t vector simply by dT=t(n)-t(n-1). Unfortunately this
turns out to be 0.0 for every other t(n) which of course produces no
change for dY which leads to ever decreasing time steps and eventually
a fatal odeint error.

The time variable (t) seems to be *absolute time*. Therefore i
modified the code so that odeint runs for only one step and used the
time variable directly for the calculation of dY as follows:

def dS(y,t,K,M,dT,c):
  dY=[0.0, 0.0, 0.0]
  F = (-K*y[1]) - (c*y[0]) + (scipy.sin(y[2]) * 1.0) #Forces acting on
the spring...
  a = F/M #- 9.81                #...cause an acceleration...
  dY[0] = a * t         #...which sets the body in motion at some velocity..
  dY[1] = y[0] * t      #...which displaces the body at every dT step
  dY[2] = 2.0*scipy.pi*4.0*t #Simple integration of a phase component
  return dY

Fs=512
dT = 1.0/Fs
y0 = [0.0,0.0,0.0]
#y0 = [0.0,0.0,0.0,0.0]
tSim = 4
nSim = tSim*Fs

theTime=range(0,nSim)
S = scipy.zeros((nSim,len(y0)))
k = 1
S[0,:] = y0
while k<nSim:
  S[k,:] = SI.odeint(dSWithGravWithDriveAndTime,S[k-1,:],[0,1],
args=(32.0,4.0,dT,0.0))[1,:]
  k+=1
pylab.plot(theTime,S[:,1]);

Now, obviously this is much slower, but i would not mind if the output
was correct....And it is not :-(

Furthermore i noticed that for one step t(n),t(n+1), the function dS
might be called 5-6 times with different values for t. In this case
where the function is called for every time step, direct
multiplication with the time value should not matter. It would just
return different slopes because the varying t is now effectively a
varying dT.

I understand that odeint might be obtaining intermediate values with a
varying integration step for its particular integration scheme to work
(something like the intermediate values obtained at dT/2 in a Runge
Kutta scheme (?) ) but what is this integration scheme and how can i
use (t) effectively?

(Please note that i have also tried to set h0, hmin, hmax to dT but
still, odeint dies a horrible death :-/ )