我仍在尝试了解solve_ivp如何对抗odeint,但正当我掌握它的窍门时,发生了一些事情。
我正在尝试求解非线性摆的运动。有了 odeint,一切都像魅力一样,在solve_ivp 上却发生了一些奇怪的事情:
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp, odeint
g = 9.81
l = 0.1
def f(t, r):
omega = r[0]
theta = r[1]
return np.array([-g / l * np.sin(theta), omega])
time = np.linspace(0, 10, 1000)
init_r = [0, np.radians(179)]
results = solve_ivp(f, (0, 10), init_r, method="RK45", t_eval=time) #??????
cenas = odeint(f, init_r, time, tfirst=True)
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.plot(results.t, results.y[1])
ax1.plot(time, cenas[:, 1])
plt.show()
我错过了什么?
最佳答案
这是一个数值问题。 solve_ivp
的默认相对和绝对容差分别为 1e-3 和 1e-6。对于许多问题,这些值太大,应该给出更严格的误差容限。 odeint
的默认相对容差是 1.49e-8。
如果将参数 rtol=1e-8
添加到 solve_ivp
调用中,则绘图一致:
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp, odeint
g = 9.81
l = 0.1
def f(t, r):
omega = r[0]
theta = r[1]
return np.array([-g / l * np.sin(theta), omega])
time = np.linspace(0, 10, 1000)
init_r = [0, np.radians(179)]
results = solve_ivp(f, (0, 10), init_r, method='RK45', t_eval=time, rtol=1e-8)
cenas = odeint(f, init_r, time, tfirst=True)
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.plot(results.t, results.y[1])
ax1.plot(time, cenas[:, 1])
plt.show()
剧情:
关于python-3.x - 使用scipy的solve_ivp求解非线性摆运动,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/56153628/