使用来自 solve_ivp 和 odeint 的解决方案的曲线拟合差异
Discrepancy in curve fit using solutions from solve_ivp and odeint
这是我试图拟合示例数据的基本示例方程。
目标是为我的数据找到 k
的最佳拟合,假设数据遵循上述等式。一个明显的方法是对方程进行数值积分,然后使用曲线拟合方法最小化最小二乘并得到 k
.
使用 odeint
and ivp_solve
and using them on curve_fit
进行数值积分产生了相当大的差异。与较新的 solve_ivp
相比,较旧的 odeint
更适合。 k
的最佳拟合值也非常不同。
Best fit k using ivp = [2.74421966]
Best fit k using odeint = [161.82220545]
这背后的原因可能是什么?
## SciPy, NumPy etc imports here
N_SAMPLES = 20
rnd = np.random.default_rng(12345)
times = np.sort(rnd.choice(rnd.uniform(0, 1, 100), N_SAMPLES))
vals = np.logspace(10, 0.1, N_SAMPLES) + rnd.uniform(0, 1, N_SAMPLES)
def using_ivp(t, k):
eqn = lambda t, x, k: -k * x
y0 = vals[0]
sol = solve_ivp(eqn, t, y0=[y0], args=(k, ),
dense_output=True)
return sol.sol(t)[0]
def using_odeint(t, k):
eqn = lambda x, t: -k * x
y0 = vals[0]
sol = odeint(eqn, y0, t)
return sol[:,0]
tfit = np.linspace(min(times), max(times), 100)
#Fitting using ivp
k_fit1, kcov1 = curve_fit(using_ivp, times, vals, p0=1.3)
fit1 = using_ivp(tfit, k_fit1)
#Fitting using odeint
k_fit2, kcov2 = curve_fit(using_odeint, times, vals, p0=1.3)
fit2 = using_odeint(tfit, k_fit2)
plt.figure(figsize=(5, 5))
plt.plot(times, vals, 'ro', label='data')
plt.plot(tfit, fit1, 'r-', label='using ivp')
plt.plot(tfit, fit2, 'b-', label='using odeint')
plt.legend(loc='best');
print('Best fit k using ivp = {}\n'\
'Best fit k using odeint = {}\n'.\
format(k_fit1, k_fit2))
再次检查 solve_ivp
的输入参数是什么。积分区间由 t_span
参数中的前两个数字给出,因此在您的应用程序中,sol.sol(t)
中的大多数值都是通过外推获得的。
通过将间隔指定为 [min(t),max(t)]
.
来更正
要获得更兼容的计算,请明确设置误差容限,因为默认值不必相等。例如 atol=1e-22, rtol=1e-9
这样只有相对公差有影响。
很好奇你是如何使用args
机制的。它最近才被引入 solve_ivp
,以便与 odeint
更加兼容。在这两种情况下我都不会使用它,因为参数的定义及其使用包含在一个 3 行代码块中。它有它的用途,适当的封装和与其他代码的隔离是一个真正的问题。
这是我试图拟合示例数据的基本示例方程。
目标是为我的数据找到 k
的最佳拟合,假设数据遵循上述等式。一个明显的方法是对方程进行数值积分,然后使用曲线拟合方法最小化最小二乘并得到 k
.
使用 odeint
and ivp_solve
and using them on curve_fit
进行数值积分产生了相当大的差异。与较新的 solve_ivp
相比,较旧的 odeint
更适合。 k
的最佳拟合值也非常不同。
Best fit k using ivp = [2.74421966]
Best fit k using odeint = [161.82220545]
这背后的原因可能是什么?
## SciPy, NumPy etc imports here
N_SAMPLES = 20
rnd = np.random.default_rng(12345)
times = np.sort(rnd.choice(rnd.uniform(0, 1, 100), N_SAMPLES))
vals = np.logspace(10, 0.1, N_SAMPLES) + rnd.uniform(0, 1, N_SAMPLES)
def using_ivp(t, k):
eqn = lambda t, x, k: -k * x
y0 = vals[0]
sol = solve_ivp(eqn, t, y0=[y0], args=(k, ),
dense_output=True)
return sol.sol(t)[0]
def using_odeint(t, k):
eqn = lambda x, t: -k * x
y0 = vals[0]
sol = odeint(eqn, y0, t)
return sol[:,0]
tfit = np.linspace(min(times), max(times), 100)
#Fitting using ivp
k_fit1, kcov1 = curve_fit(using_ivp, times, vals, p0=1.3)
fit1 = using_ivp(tfit, k_fit1)
#Fitting using odeint
k_fit2, kcov2 = curve_fit(using_odeint, times, vals, p0=1.3)
fit2 = using_odeint(tfit, k_fit2)
plt.figure(figsize=(5, 5))
plt.plot(times, vals, 'ro', label='data')
plt.plot(tfit, fit1, 'r-', label='using ivp')
plt.plot(tfit, fit2, 'b-', label='using odeint')
plt.legend(loc='best');
print('Best fit k using ivp = {}\n'\
'Best fit k using odeint = {}\n'.\
format(k_fit1, k_fit2))
再次检查 solve_ivp
的输入参数是什么。积分区间由 t_span
参数中的前两个数字给出,因此在您的应用程序中,sol.sol(t)
中的大多数值都是通过外推获得的。
通过将间隔指定为 [min(t),max(t)]
.
要获得更兼容的计算,请明确设置误差容限,因为默认值不必相等。例如 atol=1e-22, rtol=1e-9
这样只有相对公差有影响。
很好奇你是如何使用args
机制的。它最近才被引入 solve_ivp
,以便与 odeint
更加兼容。在这两种情况下我都不会使用它,因为参数的定义及其使用包含在一个 3 行代码块中。它有它的用途,适当的封装和与其他代码的隔离是一个真正的问题。