Можно ли scipy.integrate.odeint вывести внутренние расчеты

Question

Мне нужно передать внутренние вычисления от вызова к odeint.Обычно я заново вычисляю значения после завершения интеграции, но я бы предпочел выполнить все вычисления в вызываемой функции для odeint.

Мои проблемы не являются сложными в вычислительном отношении, поэтому приемлемо немного увеличить производительность при выполнении вычислений внутри решателя од.

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)

Этот код не запускается, но, надеюсь, покажетчто я послеЯ хочу получить значение «internal_calc» из функции eom для каждого прохода через интегратор.

Я искал варианты, но один из лучших известных мне программистов на Python сказал мне написать свой собственный интегратор, чтобы я мог делать то, что хочу.

До того, как я это сделал, яЯ хотел бы спросить, есть ли у кого-нибудь еще метод для получения значений из решателя odeint.

Answer 1

Возможно, вы просто не можете использовать возвращаемые значения вашей eom функции.Поэтому вам нужен какой-то другой способ контрабанды данных из eom.Есть много, много разных способов сделать это.Возможно, самое простое - просто использовать глобальную переменную:

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)

Вывод:

182

Также обратите внимание, что в приведенном выше коде я использовал solve_ivp вместо odeint.odeint документы говорят, что при написании нового кода вы должны теперь использовать solve_ivp вместо старого odeint.

Если бы это был мой собственный код, я бы, вероятно,выполнить задачу, передав объект-аккумулятор в частичную версию моей функции:

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)

Вывод:

182

Однако, если вы просто пишете короткий скрипт или что-то в этом роде,Вы, вероятно, должны просто использовать более простой подход global.Это довольно хороший пример использования.