在 python 中查找过零

Question

我写了下面的代码来查看 t 我的 ODE "exponential_decay" 穿过零线。这是 Brent Method.

odr, hr, dr, cr, m = np.genfromtxt('data.txt',unpack=True)
n=0
with open('RDE_nob_trans.txt', 'w') as d:
  for i in range(len(dr)):
      c = cr[i]
      initp = dr[i]
      exponential_decay = lambda t, y: -(1/(1+t)) * (2 *(1-y)* (-2 + (y/c)) + 3 - 3 * y)
      t_span = [0, 1] # Interval of integration
      y0 = [initp] # Initial state: y(t=t_span[0])=2
      desired_answer = odr[i]
      sol_ode = solve_ivp(exponential_decay, t_span, y0) # IVP solution

      f_sol_ode = interp1d(sol_ode.t, sol_ode.y) # Build interpolated function
      ans = brentq(lambda x: f_sol_ode(x) - desired_answer, t_span[0], t_span[1])
      d.write("{0} {1} {2} {3} {4}\n".format(hr[i], dr[i], cr[i], m[i], ans))

在这段代码中我们知道初始点initp = dr[i]，我们知道微分方程在过零处的值desired_answer = odr[i]，我们愿意找到其中t我们有这个答案。没关系，我们通过这段代码得到了答案。 ans 是零交叉处的 t。

我的问题是，如果我们现在的 ODE 答案 desired_answer = odr[i] 不仅仅是一个数字，而且是 t 的值，我们应该怎么做。

我的意思是使用 odr[i] 我们正在读取数据文件并随后读取数字。现在考虑我们有一些像 odr = 0.1 * t, 0.12 *t, 0.23 *t etc. 这样的东西，它不是一个数字，而是一个 t.

的函数

Answer 1

这不是 solve_ivp 界面的最有效使用方式。您可以使用 event 机制自动获取结果。

def event(t,y): return y-desired_answer;
event.terminal = True;

sol_ode = solve_ivp(exponential_decay, t_span, y0, events=event) # IVP solution
ans = sol_ode.t[-1]

即使您想使用自己的求解器（或求解器调用），您也可以使用密集输出选项获得 method-specific 和准确的分段多项式插值，

sol_ode = solve_ivp(exponential_decay, t_span, y0, dense_output=True) # IVP solution
f_sol_ode = sol_ode.sol

要获得依赖于时间的函数的根，您只需编写时间依赖性代码，例如

def event(t,y): return y-0.23*t;

或

ans = brentq(lambda t: f_sol_ode(t) - 0.23*t, t_span[0], t_span[1])

如果您想要可靠的结果，您应该明确控制公差 atol, rtol。

在 python 中查找过零

Finding zero crossing in python

python

ode

genfromtxt

odeint