[SciPy-User] scipy odeint question

[Luke]

I am generating some ODE's that I integrate with scipy.odeint, and
everything is working well, but there is something that I would like
to do that I'm not sure if odeint can currently handle.

The function representing the right hand side of the ODE's I generate
are of the standard form that scipy needs, namely:
def f(t, x, params):
   ....
   return dxdt

In the ... there are many (sometimes thousands or tens of thousands of
lines) intermediate variables that are used for the final calculation
of dxdt.  The calculations generally involve the basic arithmetic
operators (+-*/**) and some trig functions (sin, cos, tan).
There are several other functions of t, x, and params that I
would like to evaluate *only at the integration step times specified
by the user*, and I would like to be able make use of the fact that
the computation of the intermediate variables has already been done
for the time step that is to be returned in the x array.

Here is an example of what I am trying to acheive:

def f(t, x, params):
   z1 = ...
   z2 = ...   (z1, z2, z3 would be in reality be more like z1, ...,
z1000, all functions of t, x, params)
   ...
   zn = ...
   dxdt = [z1+z2**2, z2*z3-z1**3, ...]
   # The following only to be evaluated at time points in the user
specified vector, e.g.,  t = arange(t_i, t_f, t_s)
   eqn_1 = f_1(z1,...., zn)
   eqn_2 = f_2(z1, ..., zn)
   ....
   eqn_n = f_n(z1, ..., zn)
   return dxdt

>From a practical point of view, I can easily integrate the equations,
determine the state trajectory, then call a function that evaluates
the other quantities eqn1, eqn2... *after* numerical integration.  But
for certain applications, such as realtime control applications, it is
useful to avoid repeated calculations (as would happen if I just
called a function that took that state trajectory and calculated the
eqn1, eqn2).

I can see a couple of approaches that might work:
1)  Make the time vector a global variable, and use an if statement
inside of f(t, x, params) like:
if t in t_span:
   eqn_1 = ...
   eqn_2 = ...
   ....
   eqn_n = ...
   accumulate_eqns(eqn1, ..., eqn_n)   # accumulate_eqns could be
either a function passed to f through pararms, or a global function

2) Make eqn_1, ..., eqn_n globals, and only integrate one time step at
a time so that they can be stored at the end of each integration step.
 In this way, odeint would be called in a for loop that would use a
time vector that had only two elements (their difference would be the
step time specified by the user)

I'm not 100% sure the above methods would work because I don't know if
globals in the namespace which calls odeint would be viewable to
globals declared in the function f.

What would be really awesome would be if odeint could hold on to any
variables in the f namespace when it reaches a time step that is
requested, and then call another function to evaluate the output
equations, but that function would have access to the state of all the
variables in the f namespace *at the time point* where the state x is
to be evaluated.  Ideally, the extra equations wouldn't be computed
except at the time points which are specified by the time vector
specified by the user.

Examples of use would be computing the linear/angular momentum of a
system of rigid bodies and particles, center of mass of the bodies,
total energy (kinetic and potential), and the value of constraints /
conserved quantities during the numerical integration of the
equations of motion.  The equations of motion can be integrated
without computing these quantities, but they are very useful for
checking that the numerical integration is well behaved.  For
complicated systems such as spacecraft, these equations can be
extremely long, and it can be very computationally advantageous to
make use of these intermediate quantities rather than having to be
recalculated after the numerical integration is complete.

Any ideas?

Sincerely,
~Luke