[SciPy-User] Odeint: not following input time

[Warren Weckesser]

On Fri, Jun 10, 2011 at 8:02 AM, Jeroen Meidam <j.meidam at gmail.com> wrote:

> 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.
>


(Look again: the times printed from within the function are not the same as
the values in 't', and if you change 'wx' and 'wy', it will likely print
different times. )



> 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?
>
>

Jeroen,

This is normal.  The odeint function (which is a wrapper for the Fortran
library LSODA) uses adaptive step sizes.  Internally, it tries to take as
large a step as it can while keeping the estimate of the error within some
bound.  For well-behaved, slowly varying functions, it is able to take very
large steps.

When you print the time from within your function, you are seeing the times
at which the algorithm is evaluating the right-hand side of your equation.
In general, these times will not be at the times requested in the call to
odeint.  Interpolation is used to fill in the data at the times requested.
The error associated with interpolation is no worse than the error
associated with the underlying solver algorithm, so this is not a problem.

If you look at the solution 'sol', you'll see that it has shape (500, 4).
That's one row for each time requested, and all the rows are the same; i.e.
the solution is constant, as expected.

Warren



>
> ***** 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
>
>
> ******************
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110610/904185ea/attachment.html>