0

I would like to adapt an initial-value-problem (IVP) to a boundary-value-problem (BVP) using scipy.integrate.solve_bvp. A similar question was asked here, but I do not follow everything explained in the answer. The example below regarding the SIR model was taken from this website. Here, the initial condition y0 is taken to be the initial value of S, I, and R at time x0[0]. This system of ODEs is given by the function SIR below, which returns [dS/dt, dI/dt, dR/dt] over the interval from x[0] to 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)

As a sanity-check, running the code above produces the figure below:

IVP example

Now, suppose I would like to add another boundary condition - say x1 and y1 - to be evaluated at 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

From the documentation of solve_bvp, it appears that bc must be callable boundary conditions. The other parameters of solve_ivp and solve_bvp also appear different. How can I use this toy-example to solve a BVP in this way?

  • 1
    #(boundary conditions) = #(ODE equations) + #(variable parameters). So for example, if you declare the two parameters as variable, then you get 5 slots for boundary conditions, that you would have to select from among the 6 values that you propose. – Lutz Lehmann Mar 30 '20 at 06:35
  • From what you present it seems more likely that you want to fit the model to your data. That you would do with one of the minimization or fitting functions from scipy.optimize, like presented in http://adventuresinpython.blogspot.com/2012/08/fitting-differential-equation-system-to.html. One could clean up and simplify that code a little, but not much. – Lutz Lehmann Mar 30 '20 at 06:37
  • @LutzLehmann Regarding the first comment, am I correct in thinking that the function `SIR(t, y, prms)` is my boundary conditions `bc`? Regarding the second comment, I see `ds = integrate.odeint(funct,initial_cond,t)` inside the for-loop - this looks like a shooting method. I would like to get there, but the math involved in shooting methods is a little above me right now. I would prefer to understand why the solution works instead of copy-pasting an algorithm. –  Mar 30 '20 at 06:41
  • Ie, the fitting part comes later; the procedure is similar to that presented in the example linked in the problem description. The idea is to use solve_bvp to get a callable function that yields y_est such that residuals can be calculated against y_obs; then I can use scipy minimize to get beta and gamma. –  Mar 30 '20 at 06:48
  • Making a detour over a BVP formulation makes little sense. The solution curve has 5 parameters, the 2 coefficients and 3 initial or boundary conditions. For the fitting task it does not matter where these are located, in theory one could even use other integration constants that are not directly values of the solution. A BVP solution is always more costly than and ODE integration. – Lutz Lehmann Mar 30 '20 at 09:14
  • See https://stackoverflow.com/q/59379143/3088138 for a similar question of converting an ODE IVP into a BVP. – Lutz Lehmann Mar 30 '20 at 09:24
  • One fitting experiment finds parameters `[beta, gamma] = [0.00219851, 0.45259793]` with an average deviation of 20. – Lutz Lehmann Apr 01 '20 at 18:46

0 Answers0