2

I have a need to pass internal calculations from a call to odeint. I normally recalculate the values again after I am finished with the integration, but I would prefer to do all of the calculations in the called function for odeint.

My problems are not computationally heavy, so taking a little extra performance hit in doing calculations inside of the ode solver is acceptable.

from scipy.integrate import odeint
import numpy as np

def eom(y, t):
    internal_calc = y/(t+1)
    xdot = y

    return xdot, internal_calc

if __name__ == '__main__':

    t = np.linspace(0, 5, 100)
    y0 = 1.0  # the initial condition

    output, internal_calc = odeint(eom, y0, t)

This code doesn't run, but hopefully shows what I am after. I want to get the 'internal_calc' value out of the eom function for each pass through the integrator.

I've looked around for options, but one of the best python programmers I know told me to write my own integrator so I can do what I want.

Before I did that, I thought I would ask if anyone else has a method for getting values out of the odeint solver.

Chris
  • 23
  • 3

1 Answers1

0

It is possible, you just can't use the return values of your eom function. So you need some other way of smuggling data out from eom. There are many, many different ways to do this. The simplest is probably to just use a global variable:

import scipy.integrate as spi

count = 0

def pend(t, y):
    global count

    theta, omega = y
    dydt = [omega, -.25*omega - 5*np.sin(theta)]

    count += 1
    return dydt

sol = spi.solve_ivp(pend, [0, 10], [np.pi - 0.1, 0.0])
print(count)

Output:

182

Also, note that I used solve_ivp instead of odeint in the code above. The odeint docs say that when writing new code, you should now use solve_ivp instead of the older odeint.

If it were my own code, I'd probably accomplish the task by passing an accumulator object into a partial version of my function:

class Acc:
    def __init__(self):
        self.x = 0

    def __str__(self):
        return str(self.x)

def pend_partial(acc):
    def pend(t, y):
        theta, omega = y
        dydt = [omega, -.25*omega - 5*np.sin(theta)]

        acc.x += 1
        return dydt
    return pend

count = Acc()
sol = spi.solve_ivp(pend_partial(count), [0, 10], [np.pi - 0.1, 0.0])
print(count)

Output:

182

However, if you're just writing a short script or something, you should probably just use the simpler global approach. This is a pretty good use case for it.

tel
  • 13,005
  • 2
  • 44
  • 62
  • ahh. Interesting. I hadn't considered using a global variable in this case. I've generally avoided global variables in other languages, and I'm guessing that is still a good practice with python, but this is a case where I seems clean. – Chris Dec 10 '18 at 14:51
  • You mentioned other ways of doing it. Could you list a couple of thoughts on how to best do that? Maybe just having a direction would help me come up with the method for doing it without global variables. – Chris Dec 10 '18 at 14:53
  • Actually, now that I've learned a little more and some of the keywords for what I am trying to do it seems other folks have found solutions as well. I'm going to look at them some more, but if anyone else has the same issues they can take a look at these as well. [link](https://stackoverflow.com/a/46815621/10768739), [link](https://stackoverflow.com/a/16921932/10768739) – Chris Dec 10 '18 at 15:02
  • I added another example that gives a pattern that would be appropriate to use in a larger project. However, for simple short scripts, you should probably just use a `global`. You're right that this is a pretty good use case for one. – tel Dec 10 '18 at 15:13