我写了下面的代码来看看t
我的ODE“指数衰变”在哪条线上穿过了零线。这是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]
,我们知道微分方程在过零点{ans
是零交叉处的t
。在
我的问题是,如果我们的ODE的答案现在是desired_answer = odr[i]
不仅仅是一个数字,而且是t
的值,我们该怎么办。在
我的意思是使用odr[i]
读取数据文件,然后读取数字。现在我们考虑一下像odr = 0.1 * t, 0.12 *t, 0.23 *t etc.
这样的东西,它不是一个数字,是t
的函数。在
这不是
solve_ivp
接口的最有效使用。您可以使用event
机制自动获得结果。在即使您想使用自己的解算器(或解算器调用),也可以使用密集输出选项获得特定于方法的精确分段多项式插值
^{pr2}$要获得依赖于时间的函数的根,只需编写时间依赖性的代码,例如
或者
如果您想要可靠的结果,您应该显式地控制公差
atol, rtol
。在相关问题 更多 >
编程相关推荐