Использование адаптивных размеров шага с помощью scipy.integrate.ode

(Краткая) документация для scipy.integrate.ode говорит, что два метода (dopri5 и dop853) имеют управление шагами и плотный вывод. Глядя на примеры и сам код, я могу видеть только простой способ получить результат от интегратора. А именно, похоже, что вы просто надавите интегратора вперед на какой-то фиксированный dt, получите значение функции в это время и повторите.

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

EDIT: плотный выход

Связанное с ним понятие (почти противоположное) - это "плотный выход", при котором сделанные шаги имеют значение, которое должен выполнять шаговый, но значения функции интерполируются (как правило, с точностью, сопоставимой с точностью шагателя), к чему-либо ты хочешь. Fortran лежащий в основе scipy.integrate.ode, по- видимому способен к этому, но ode не имеет интерфейса. odeint, с другой стороны, основан на другом коде и, очевидно, делает плотный вывод. (Вы можете выводить каждый раз, когда вы вызываете свою правую сторону, чтобы видеть, когда это произойдет, и посмотрите, что это не имеет никакого отношения к времени вывода.)

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

EDIT 2: Плотный выход снова

По-видимому, scipy 1.0.0 представил поддержку для плотного вывода через новый интерфейс. В частности, они рекомендуют отходить от scipy.integrate.odeint и к scipy.integrate.solve_ivp, который как ключевое слово dense_output. Если установлено значение True, возвращаемый объект имеет атрибут sol который вы можете вызвать с помощью массива раз, который затем возвращает значения интегрированных функций в это время. Это еще не решает проблему для этого вопроса, но во многих случаях это полезно.

Ответ 1

Поскольку SciPy 0.13.0,

Промежуточные результаты из семейства dopri семейств решателей ODE теперь можно получить с помощью solout обратного вызова solout.

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def logistic(t, y, r):
    return r * y * (1.0 - y)

r = .01
t0 = 0
y0 = 1e-5
t1 = 5000.0

backend = 'dopri5'
# backend = 'dop853'
solver = ode(logistic).set_integrator(backend)

sol = []
def solout(t, y):
    sol.append([t, *y])
solver.set_solout(solout)
solver.set_initial_value(y0, t0).set_f_params(r)
solver.integrate(t1)

sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()

Результат:

Результат, похоже, немного отличается от Тима Д., хотя оба они используют один и тот же бэкэнд. Я подозреваю, что это связано с свойством dopri5. В подходе Тима я думаю, что результат k7 с седьмого этапа отбрасывается, поэтому k1 вычисляется заново.

Примечание. Известная ошибка с set_solout не работает, если вы установите ее после установки начальных значений. Он был установлен как SciPy 0.17.0.

Ответ 2

Я смотрел на это, чтобы попытаться получить тот же результат. Оказывается, вы можете использовать хак для получения пошаговых результатов, установив nsteps = 1 в экземпляр ode. Он будет генерировать UserWarning на каждом шаге (это можно поймать и подавить).

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
import warnings


def logistic(t, y, r):
    return r * y * (1.0 - y)

r = .01
t0 = 0
y0 = 1e-5
t1 = 5000.0

#backend = 'vode'
backend = 'dopri5'
#backend = 'dop853'

solver = ode(logistic).set_integrator(backend, nsteps=1)
solver.set_initial_value(y0, t0).set_f_params(r)
# suppress Fortran-printed warning
solver._integrator.iwork[2] = -1

sol = []
warnings.filterwarnings("ignore", category=UserWarning)
while solver.t < t1:
    solver.integrate(t1, step=True)
    sol.append([solver.t, solver.y])
warnings.resetwarnings()
sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()

Результат:

Ответ 3

Метод integrate принимает логический аргумент step, который сообщает методу вернуть один внутренний шаг. Однако оказалось, что решатели 'dopri5' и 'dop853' не поддерживают его.

В следующем коде показано, как вы можете получить внутренние шаги, решаемые решателем при использовании решателя 'vode':

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def logistic(t, y, r):
    return r * y * (1.0 - y)

r = .01

t0 = 0
y0 = 1e-5
t1 = 5000.0

backend = 'vode'
#backend = 'dopri5'
#backend = 'dop853'
solver = ode(logistic).set_integrator(backend)
solver.set_initial_value(y0, t0).set_f_params(r)

sol = []
while solver.successful() and solver.t < t1:
    solver.integrate(t1, step=True)
    sol.append([solver.t, solver.y])

sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()

Результат:

Ответ 4

FYI, хотя ответ уже принят, я должен указать на историческую запись, что плотный вывод и произвольная выборка из любой точки по вычисленной траектории поддерживаются в PyDSTool. Это также включает в себя запись всех адаптивно определяемых временных шагов, используемых внутренне решателем. Это взаимодействует с dopri853 и radau5 и автоматически генерирует код C, необходимый для взаимодействия с ними, вместо того, чтобы полагаться на (намного медленнее) обратные вызовы функции python для определение правой части. Насколько мне известно, ни одна из этих функций изначально или эффективно не предоставляется ни в одном другом решателе, ориентированном на питон.

Ответ 5

Вот еще один вариант, который также должен работать с dopri5 и dop853. В принципе, решатель будет вызывать функцию logistic() так часто, как это необходимо для вычисления промежуточных значений, так что где мы сохраняем результаты:

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt

sol = []
def logistic(t, y, r):
    sol.append([t, y])
    return r * y * (1.0 - y)

r = .01

t0 = 0
y0 = 1e-5
t1 = 5000.0
# Maximum number of steps that the integrator is allowed 
# to do along the whole interval [t0, t1].
N = 10000

#backend = 'vode'
backend = 'dopri5'
#backend = 'dop853'
solver = ode(logistic).set_integrator(backend, nsteps=N)
solver.set_initial_value(y0, t0).set_f_params(r)

# Single call to solver.integrate()
solver.integrate(t1)
sol = np.array(sol)

plt.plot(sol[:,0], sol[:,1], 'b.-')
plt.show()