[SciPy-user] odeint for calculating trajectories

[Anne Archibald]

2009/2/12 Rob Clewley <rob.clewley at gmail.com>:
>> Is there a good way to use scipy.integrate.odeint to calculate trajectories
>> from an observed velocity field? I know you can do this when you have an
>> analytic expression for dx/dt, but in this case I have a spatial grid of
>> values for dx/dt.  The only way I've come up with is to make the function
>> passed to odeint something that will interpolate fromt the grid to the given
>> point.
>
>
> I don't think odeint is the right tool for this job - there is no ODE
> integration to do if you do not have an explicit function for the
> vector field. You should think of it purely as an interpolation
> problem. You have (t,x) values and (t, dx/dt) values, so this defines
> a piecewise quadratic function which has continuous *second*
> derivative everywhere (i.e. the trajectory smoothly agrees at your
> mesh points). I would use the polynomial interpolation classes that
> were recently added to scipy by Anne Archibald (search this list for
> details about it). You pass it your arrays of values and you get back
> a function that smoothly interpolates through your points. This is the
> most accurate trajectory that you can derive from this finite mesh
> vector-field.

Put another way, odeint is full of cleverness to figure out how fast
the derivative is changing, but that is of no use to  you here (unless
you have an extremely high-resolution, slowly-changing vector field).
So an ode solver that walks along the grid is about as good as you can
do.

There is actually cython code to do exactly this in the scikit
vectorplot, which uses it to implement line integral convolution. If
you want fast and dirty trajectories, you may be interested in
modifying that code. Using a piecewise polynomial will certainly give
you a smoother trajectory, though.

Anne