Я хотел бы адаптировать проблему начальных значений ( IVP ) к проблеме граничных значений ( BVP ), используя scipy.integrate.solve_bvp
. Подобный вопрос был задан здесь , но я не следую всему объясненному в ответе. Приведенный ниже пример модели SIR был взят с с этого сайта . Здесь начальное условие y0
принимается равным начальному значению S
, I
и R
в момент времени x0[0]
. Эта система ODE задается нижеприведенной функцией SIR
, которая возвращает [dS/dt, dI/dt, dR/dt]
в интервале от x[0]
до x[-1]
.
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp, solve_bvp
def SIR(t, y, prms):
S = y[0]
I = y[1]
R = y[2]
beta, gamma = prms
# return [dS/dt, dI/dt, dR/dt]
return np.array([-beta * S * I, beta * S * I - gamma * I, gamma * I])
infected = np.array([1, 3, 6, 25, 73, 222, 294, 258, 237, 191, 125, 69, 27, 11, 4])
xt = np.arange(infected.size).astype(int)
xw = 0.2 # spacing between points on x-axis (elapsed time)
t_eval = np.arange(xt[0], xt[-1]+xw, xw)
x0 = [xt[0], xt[-1]]
y0 = [762, 1, 0] # S0, I0, R0, beginning of outbreak
N = sum(y0) # population total
prms = [0.01,0.1] # beta, gamma
sol = solve_ivp(SIR, x0, y0, method='LSODA', t_eval=t_eval, args=(prms,))
fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0], label='S')
ax.plot(sol.t, sol.y[1], label='I')
ax.plot(sol.t, sol.y[2], label='R')
ax.plot(xt, infected, label='OG data')
ax.grid(color='k', linestyle=':', alpha=0.3)
fig.legend(loc='lower center', ncol=4, mode='expand')
plt.show()
plt.close(fig)
В качестве проверки работоспособности, запустив код выше создаст рисунок ниже:
Теперь предположим, что я хотел бы добавить еще одно граничное условие - скажем, x1
и y1
- для оценки в x0[-1]
.
y0 = [0, 200, N-200] # S1, I1, R1, end of graph of outbreak; values from eye-ing the graph # N-200 \approx 550
Из документации solve_bvp
видно, что bc
должны быть вызываемыми граничными условиями. Другие параметры solve_ivp
и solve_bvp
также отличаются. Как я могу использовать этот игрушечный пример для решения BVP таким образом?