scipy.integrate.odeint: how can odeint access a parameter set that is evolving independently of it?

Question

I might have some non-linear ODEs that are being solved by scipy.integrate.odeint. However, a parameter at each time step might have to be updated by using a non-DE rule, which uses the results of the ODE solver at each timestep, together with the variable along which integration is being performed, in order to update the parameter. For example, the non-DE rule might be a coarse discretization of a "background" PDE system.

---

Here's a toy example, related to what I am working on right now: We might have a particle, which has some properties ($x$, $y$ and $z$) associated with it that evolve with time -- let's say that $x$ and $y$ represent the 2D position of the particle, and $z$ its "colour". The particle moves over a 2D space. The 2D space is discretized, and each discretized point has a property $p$ associated with it - let's say that $p$ evolves via the following rules:

1. certain pre-selected spatial points have a $p$ source
2. if a particle has recently visited a spatial point, $p$ is set to zero
3. at every time step $t$, 10 random sets (of size 4) of neighbouring spatial points are selected, and their $p$ is set to be the average $p$ of the points in the set

When calculating the RHS ($dx/dt$, $dy/dt$ and $dz/dt$), we need to know what the average of value of $p$ is on the spatial points closest to the particle's current location. However, it's not straightforward to embed the evolution of $p$ on the spatial points as a bunch of DE rules...otherwise we could have simply passed an ODE system like $dx/dt$, $dy/dt$, $dp_1/dt$, ... , $dp_n/dt$) to odeint.

---

Let's say I want ODE results at t = [t1, ..., t2]. Then, I could setup the program so that each time I call odeint, I only make it give me results between one timestep, then use those results to update the non-DE rule, before calling odeint again. However, this would be quite inefficient, since I'd be breaking odeint's internal flow?

Is there any other way to get around this issue? Could it be that I can inherit from the object-oriented version of 'odeint', and then have my class execute the non-DE step in between each step? Would this be any more efficient than the naive solution suggested above?

Answer 1

EDIT: No, that process will not be compatible with any adaptive time-stepper. It's not deterministic, and even for non-adaptive ODE solvers, your proposed process will degrade the accuracy of your solution.

---

You'll have to describe what you want to do in more detail.

In general, changing a parameter changes the equations you want to solve. If your parameter changes in a smooth (infinitely continuously differentiable, all derivatives bounded over your interval of integration) fashion, depending on this "non-differential equation rule", it's probably better to figure out how to query an update as part of your right-hand side function evaluation (and Jacobian function evaluation).

If it does not change in a smooth fashion, it depends on how many bounded, continuous derivatives you have. In general, differential equation solvers use interpolating polynomials to solve approximately differential equations, and the error term for a method of order $k$ is usually related to the derivative of order $k+1$. If you have bounded, continuous derivatives up to order $k + 1$, then you can use a method of up to order $k$. If you have no continuous derivatives, you have to use special methods that have low orders of convergence.

What will probably happen if you try to kludge in something like you're proposing is that the integrator will restart itself repeatedly, which will be inefficient, and probably inaccurate. I don't think that subclassing the object-oriented version of odeint will rectify the problem because as I've described above, the problem is a fundamental numerical analysis issue.

Answer 2

Of course you can hack something together that might work.

logger = [] # visible in odefunc

def odefunc(P,t):

    if logger: # if the list is not empty
        if logger[-1]: # then read the last values 
            pass # and do something based on them

    **your complex math here***
    some_calculated_value = 0.123    

    logger.append([t, some_calculated_value]) 

    return something

print(logger)