时间戳

首页 > 解决方案 > 尝试使用 Scipy 解决 ODE 时出错

python - 尝试使用 Scipy 解决 ODE 时出错

问题描述

我正在尝试使用 SciPy 库解决一些微分方程。到目前为止，我有一个示例，您可以在下面看到。代码运行正常，但是执行时间比较长，而且在nsteps < 10^6的情况下也会报错。

错误是：

0.1 [55.4662258]
0.2 [61.50134206]
0.30000000000000004 [68.16443669]
0.4 [75.52104176]
0.5 [83.64306141]
0.6 [92.61029686]
0.7 [102.51090657]
0.7999999999999999 [113.44178205]
0.8999999999999999 [125.51003611]
0.9999999999999999 [138.83414699]
1.0999999999999999 [153.54491757]
1.2 [169.78658954]
1.3 [187.71843357]
1.4000000000000001 [207.51637757]
1.5000000000000002 [229.37467955]
1.6000000000000003 [253.507774]
1.7000000000000004 [280.15236668]
1.8000000000000005 [309.56989072]
1.9000000000000006 [342.04876583]
2.0000000000000004 [377.90754689]
2.1000000000000005 [417.49804931]
2.2000000000000006 [273.12387997]
2.3000000000000007 [165.68845619]
2.400000000000001 [100.59106997]
2.500000000000001 [61.14730695]
2.600000000000001 [37.24734963]
2.700000000000001 [22.76575439]
2.800000000000001 [13.99098983]
2.9000000000000012 [8.6741386]
3.0000000000000013 [5.45252409]
3.1000000000000014 [3.50046623]
3.2000000000000015 [2.31766504]
3.3000000000000016 [1.60097602]
3.4000000000000017 [1.16671624]
3.5000000000000018 [1.00000093]

Warning (from warnings module):
  File "C:\Users\X\AppData\Local\Programs\Python\Python37\lib\site-packages\scipy\integrate\_ode.py", line 1009
    self.messages.get(istate, unexpected_istate_msg)))
UserWarning: vode: Excess work done on this call. (Perhaps wrong MF.)
3.600000000000002 [0.99999893]

如果我将 nsteps 增加到 10^6 或更多，它工作正常，但程序完成它的工作需要更长的时间。

可以做些什么来加快程序的执行并防止此类错误的发生？例如，最大步数不必设置为高于 10^3 可以做什么？

也许通过某种方式降低准确性？结果只能精确到小数点后 1 或 2 位。

from scipy.integrate import ode

r = 1.3
K = 1000000000
Gamma = 0.5
A_m = 6
phi_PS = 0.25
theta_PS = 0.01
delta_PS = 1
beta_BS = 0.1
N_BS = 25
L = 25


def F1(time, N_PS, par):
  if N_PS >= 1.0:
      return ((r * N_PS * (1 - (N_PS) / K) - theta_PS * N_PS - phi_PS * N_PS - par[0] * delta_PS * A_m * N_PS - Gamma * N_PS * 0.1) + beta_BS * N_BS)
  elif par[1] >= 1.0:
      return beta_BS * N_BS
  else:
      return 0


t0 = 0
y0 = 50
solver = ode(F1)
solver.set_integrator('vode', nsteps=10**5, method='bdf')
solver.set_initial_value(y0, t0)

dt = 0.1
while solver.successful() and solver.t < 20:
    time = solver.t
    next_t = time + dt
    x = 1 if (time >= 2 and time <= 4) else 0
    par = (x, L)
    solver.set_f_params(par)
    next_y = solver.integrate(next_t)
    if next_y < 1:
        L = L - 1
    print(next_t, next_y)

标签： pythonpython-3.xscipyruntime-errorode

解决方案

当 时，你有一个不连续性N_PS=1。这将导致严重的积分器问题。如果这确实是需要集成的功能，我建议使用三个独立的功能，即

def F1(time, N_PS, par):
    return ((r * N_PS * (1 - (N_PS) / K) - theta_PS * N_PS - phi_PS * N_PS - par[0] * delta_PS * A_m * N_PS - Gamma * N_PS * 0.1) + beta_BS * N_BS)

def F2(time, N_PS, par):
    return beta_BS * N_BS

def F3(time, N_PS, par):
    return 0

然后，您可以使用scipy.integrate.solve_ivpwith 函数events来检测何时在不同状态之间切换。这不应该遇到与您在此处看到的相同的不连续性问题。

推荐阅读