Как пройти через точные точки с помощьюsolve_ivp?

У меня есть система ОДУ, и я хочу изменить значение переменных, когда решатель достигает точного момента времени. То, что я пытаюсь сделать, похоже на этот пример в Джулии

То, что я пытался сделать, это использовать else и проверить, достигнуто ли время t. Однако это не работает, потому что массив, содержащий время, не проходит через эту точку. Также безуспешно пытался использовать t_eval.

Используя ту же проблему, показанную в примере с юлией, мы имеем:

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

def f(t, u):
    if t == 4 and u[0] < 4:
        u[0] += 10
    du = np.empty(1)
    du[0] = -u[0]
    return du

u0=[10.0]
V = 1
t = [0,10]
u0 = [10.0]
sol = solve_ivp(f,t,u0)

plt.plot(sol.t, sol.y[0])

Есть ли способ воспроизвести функцию обратного вызова julia?


person Eiem Lupus    schedule 19.07.2019    source источник


Ответы (1)


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

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


def f(t, u):
    du = -u
    return du


dose = 10.0

u0 = dose
V = 4

t0 = 0
t1 = 4
t2 = 10

t = np.linspace(t0, t1, 100)
sol1 = solve_ivp(f, [t0, t1], [u0], dense_output=True, t_eval=t)

u0 = sol1.y[0, -1]
if u0 / V < 4:
    u0 += dose

t = np.linspace(t1, t2, 150)
sol2 = solve_ivp(f, [t1, t2], [u0], dense_output=True, t_eval=t)

plt.plot(sol1.t, sol1.y[0], 'b')
plt.plot([sol1.t[-1], sol2.t[0]], [sol1.y[0, -1], sol2.y[0, 0]], 'b--')
plt.plot(sol2.t, sol2.y[0], 'b')
plt.grid()
plt.xlabel('t')
plt.show()

сюжет

person Warren Weckesser    schedule 20.07.2019