Адаптация начальной задачи к краевой задаче с помощью scipy.integrate.solve_bvp? - PullRequest
Новая
       111

Адаптация начальной задачи к краевой задаче с помощью scipy.integrate.solve_bvp?

0 голосов
/ 30 марта 2020

Я хотел бы адаптировать проблему начальных значений ( 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)

В качестве проверки работоспособности, запустив код выше создаст рисунок ниже:

IVP example

Теперь предположим, что я хотел бы добавить еще одно граничное условие - скажем, 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 таким образом?

...