如何使用带有时间变量的scipy.integrate.odeint来解决ODE系统

3 回答

• 3

  据我从您的问题下面的评论中了解到，您尝试合并可能有噪音的测量数据 . 您可以使用这些数据来适应您的时间课程，而不是直接插入数据 . 在这里，我显示变量 S 的结果：

  

  green dots 是从您提供的ODE系统的解决方案中采样的 . 为了模拟测量误差，我在这些数据中添加了一些噪声（ blue dots ） . 然后，您可以适合您的ODE系统，以尽可能好地再现这些数据（ red line ） .

  对于这些任务，您可以使用lmfit . 重现绘图的代码如下所示（可以在内联注释中找到一些解释）：

  # zombie apocalypse modeling
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import odeint
  from lmfit import minimize, Parameters, Parameter, report_fit
  from scipy.integrate import odeint
  
  
  # solve the system dy/dt = f(y, t)
  def f(y, t, paras):
  
      Si = y[0]
      Zi = y[1]
      Ri = y[2]
  
      try:
          P = paras['P'].value
          d = paras['d'].value
          B = paras['B'].value
          G = paras['G'].value
          A = paras['A'].value
  
      except:
          P, d, B, G, A = paras
      # the model equations (see Munz et al. 2009)
      f0 = P - B * Si * Zi - d * Si
      f1 = B * Si * Zi + G * Ri - A * Si * Zi
      f2 = d * Si + A * Si * Zi - G * Ri
      return [f0, f1, f2]
  
  
  def g(t, x0, paras):
      """
      Solution to the ODE x'(t) = f(t,x,p) with initial condition x(0) = x0
      """
      x = odeint(f, x0, t, args=(paras,))
      return x
  
  
  def residual(paras, t, data):
      x0 = paras['S0'].value, paras['Z0'].value, paras['R0'].value
      model = g(t, x0, paras)
      s_model = model[:, 0]
      return (s_model - data).ravel()
  
  # just for reproducibility reasons
  np.random.seed(1)
  
  # initial conditions
  S0 = 500.              # initial population
  Z0 = 0                 # initial zombie population
  R0 = 0                 # initial death population
  y0 = [S0, Z0, R0]     # initial condition vector
  t = np.linspace(0, 5., 100)         # time grid
  
  P = 12      # birth rate
  d = 0.0001  # natural death percent (per day)
  B = 0.0095  # transmission percent  (per day)
  G = 0.0001  # resurect percent (per day)
  A = 0.0001  # destroy percent  (per day)
  
  # solve the DEs
  soln = odeint(f, y0, t, args=((P, d, B, G, A), ))
  S = soln[:, 0]
  Z = soln[:, 1]
  R = soln[:, 2]
  
  # plot results
  plt.figure()
  plt.plot(t, S, label='Living')
  plt.plot(t, Z, label='Zombies')
  plt.xlabel('Days from outbreak')
  plt.ylabel('Population')
  plt.title('Zombie Apocalypse - No Init. Dead Pop.; No New Births.')
  plt.legend(loc=0)
  plt.show()
  
  # generate fake data
  S_real = S[0::8]
  S_measured = S_real + np.random.randn(len(S_real)) * 100
  t_measured = t[0::8]
  
  plt.figure()
  plt.plot(t_measured, S_real, 'o', color='g', label='real data')
  
  # add some noise to your data to mimic measurement erros
  plt.plot(t_measured, S_measured, 'o', color='b', label='noisy data')
  
  # set parameters including bounds; you can also fix parameters (use vary=False)
  params = Parameters()
  params.add('S0', value=S0, min=490., max=510.)
  params.add('Z0', value=Z0, vary=False)
  params.add('R0', value=R0, vary=False)
  params.add('P', value=10, min=8., max=12.)
  params.add('d', value=0.0005, min=0.00001, max=0.005)
  params.add('B', value=0.01, min=0.00001, max=0.01)
  params.add('G', value=G, vary=False)
  params.add('A', value=0.0005, min=0.00001, max=0.001)
  
  # fit model
  result = minimize(residual, params, args=(t_measured, S_measured), method='leastsq')  # leastsq nelder
  # check results of the fit
  data_fitted = g(t, y0, result.params)
  
  plt.plot(t, data_fitted[:, 0], '-', linewidth=2, color='red', label='fitted data')
  plt.legend()
  # display fitted statistics
  report_fit(result)
  
  plt.show()
  

  回复于 2024-05-02T10:47:52+08:00
• 0

  您无法在数值积分器评估ODE函数的哪些点知道先验 . 积分器（ odeint 和其他非明确的"fixed step-size"）动态生成一个内部点列表，这些点可能比给定的采样点列表具有更小或有时更大的步长 . 输出的值从内部列表中插入 .

  如果要用函数替换ODE的一部分，则必须将样本数据转换为函数 . 这可以通过插值来完成 . 使用scipy.interpolate.interp1函数生成函数对象，然后可以像任何其他标量函数一样使用它们 .

  回复于 2024-05-02T10:47:52+08:00
• 1

  要专门做我在问题中提出的问题，即使用值代替其中一个ODES你将需要使用一个循环，你使用odesolver来解决你的系统1秒，然后将输出作为下一次循环迭代的初始条件 . 这种方法的代码如下 . 然而正如许多人所指出的那样，在大多数情况下，如Cleb和其他人所描述的那样使用插值会更好

  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import odeint
  
  Si = [345, 299, 933, 444, 265, 322] # replaced an equation with list
  
  #Parameters
  P = 0           # birth rate
  d = 0.0001  # natural death percent (per day)
  B = 0.0095  # transmission percent  (per day)
  G = 0.0001  # resurect percent (per day)
  A = 0.0001  # destroy percent  (per day)
  
  # solve the system dy/dt = f(y, t)
  def f(y, t):
      Zi = y[0]
      Ri = y[1]
      # the model equations (see Munz et al. 2009)
      f0 = B*Si*Zi + G*Ri - A*Si*Zi
      f1 = d*Si + A*Si*Zi - G*Ri
      return [f0, f1]
  
  # initial conditions
  Z0 = 0                      # initial zombie population
  R0 = 0                      # initial death population
  y0 = [Z0, R0]   # initial condition vector
  # a timestep of 1 forces the odesolve to use your inputs at the beginning and provide outputs at the end of the timestep. 
  # In this way the problem that LutzL has described is avoided.
  t  = np.linspace(0, 1, 2)       
  Si =np.array(Si).T
  
  #create a space for storing your outputdata
  dataZ =[]
  dataR =[]
  #use a for loop to use your custom inputs for Si
  for Si in Si:
      y0 = [Z0, R0]
      soln = odeint(f, y0, t)
      Z = soln[:, 0]
      R = soln[:, 1]
      #define your outputs as the initial conditions for the next iteration of the loop
      Z_0 = Z[1]
      R_0 = R[1]
      #store your outputs 
      dataZ.append(Z[1])
      dataR.append(R[1])
  
  
  print (dataZ)
  print (dataR)
  

  回复于 2024-05-02T10:47:52+08:00