[SciPy-User] making use of args in odeint

Luke hazelnusse@gmail....
Mon Oct 19 16:58:57 CDT 2009


I am integrating some ODE's and want to calculate some extra quantities
during numerical integration.  I was hoping to use the args option of odeint
as a way to parameters which could be 'output' at each time step, but it
isn't working the way I would expect.  Here is a short example of
integrating the equation:
dx/dt = -a*x

where a is a parameter, and we also want to calculate at every time step:
b = t / 2.

This is just a silly example, the point is to illustrate that during
numerical integration, other quantities besides the state trajectory x(t)
are of interest.

If we wanted to integrate this from the initial condition x = 1, with a=1,
on the interval t=arange(0, 1.2, 0.2), we would expect that the trajectory
of b to simply be arange(0, 0.6, 0.1).  But instead, b ends up being not on
these exact times.

from numpy import sin, cos, arange, zeros, array
from scipy.integrate import odeint

def f(x, t, params):
    a = params[0]
    params[1] = t / 2.
    return -a*x

# Intial time, step time, final time
t_i = 0.0
t_s = 0.2
t_f = 1.0

# Initial condition
xi = 1.0
# Time array
t = arange(t_i, t_f + t_s, t_s)
n = len(t)

a = 1.0
# array for storing b
b = zeros(n)
params = [a, b[0]]

# State trajectory
xs = zeros(n)

for i, ti in enumerate(t[1:]):
    x = odeint(f, xi, [t[i], t[i+1]], args=(params,))
    xs[i+1] = xi = x[-1]
    b[i+1] = params[1]     #store the output parameters

print t
print xs
print b


The output is:
t = [ 0.   0.2  0.4  0.6  0.8  1. ]
x = [ 0.          0.81873077  0.67032008  0.54881168  0.44932901  0.3678795
]
b = [ 0.          0.10232638  0.20562135  0.30934571  0.41352631
0.51818916]

#  b is close to the times of [0, 0.1, 0.2, 0.3, 0.4, 0.5], not exact though

I would then expect that b would just have the same values as the array t /
2, but instead the time values are not on the points I specified.  My guess
is that this has to do with the internal step size selection and that the
last time that the right hand sides is evaluated is at a time point
dictacted by the internal time step, rather than the time step that I would
like it to be, namely, the same points as in the array of time points.

One solution I can see is to just store different time array to use with the
plotting of all the output parameters, but I was wondering if there was an
alternative way that would just let me ensure that the output params gets
assigned values corresponding with time values at the final time of the
interval I pass to odeint.

Does anybody know if there is a way to do this?

The reason for even taking this approach is that there are lot of
repetitious calculations that are done in more complicated systems that I
would like to avoid doing in a separate step after numerical integration.
These calculations are being done regardless In this example this isn't
apparent, but in other systems, the calculations are more significant and so
there is some computational savings to be had from this approach.

~Luke
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20091019/2da0b379/attachment.html 


More information about the SciPy-User mailing list