Эквивалент "OutputFcn" для Matlab в Scipy?

Question

Как и у решателей MATLAB ODE есть Outputfcn, где решатель ODE вызывает функцию после каждого успешного шага по времени

options=odeset('OutputFcn',@odeprint)
    [T,Y]=ode15s(@(T,X)sys(T,X,vin),[t0 .0005],X(:,1),options)

Есть ли какая-либо эквивалентная функция вывода в resol_ivp от scipy? который вызывается после каждого временного шага для сохранения вектора решения? Или есть ли какой-нибудь способ создать эту выходную функцию в execute_ivp Сципи?

Для моего кода: определение A, R как 4*4 и B, S как 4*1 Z as 4 * 1`

def func(t,Z):
  if  X[0]>5 or X[1]>0:
      Zdot=A*Z+B*U
  else:
      Zdot=R*Z+S*U
return Zdot

sol = solve_ivp(func, tspan,Z0)
aa=sol.t
bb=sol.y
X=v1.dot(bb)[:][i]
plt.plot(aa,X[0,:],'r')

У меня проблема с определением X=v1.dot(bb)[:][i], потому что я заранее не знаю длину bb. Это то, что я хочу сделать. Хотя этот код не работает. На каждом шаге решатель, я хочу X как массив 4 * 1, используя X=v1.dot(bb), где v1 - другая матрица 4 * 4

Answer 1

Если вы не зададите опцию t_eval, структура решения содержит в полях .t и .y внутренние точки шага. Если вы передадите опцию dense_output=True, то структура решения содержит в .sol кусочно-полиномиальную функцию интерполяции, которую вы можете использовать для получения решения в любое другое время или по списку раз.

Если вы хотите Для большего контроля вам, вероятно, потребуется напрямую использовать базовые классы степпера для метода интеграции.

---

Если вы хотите переключать модели y'=f_1(y) и y'=f_2(y) вдоль некоторой поверхности g(y)=0, тогда вы можете в принципе используйте

def odefun(t,y):
    return f_1(y) if g(y)>0 else f_2(y);

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

def event(t,y): x=v.dot(y); return max(x[0]-5,x[1]) 

event.terminal=True
event.direction = -1 if event(t0,y0)>0 else 1

, а затем внутри l oop после каждого события поменяйте знак направления и выберите правильную функцию

event.direction *=-1
odefun = f_1 if event.direction<0 else f_2

Затем добавьте код в l oop, чтобы добавить последнее решение к ранее собранным данным решения, см. { ссылка } для возможного решения.