使用scipy的solve_ivp解决非线性摆运动 [英] Using scipy's solve_ivp to solve non linear pendulum motion
问题描述
我仍在尝试了解resolve_ivp如何与odeint结合使用,但就在我掌握了它的本质的同时,发生了一些事情.
I am still trying to understand how solve_ivp works against odeint, but just as I was getting the hang of it something happened.
我正在尝试解决非线性摆的运动.有了odeint,一切都会像个魅力一样,在solve_ivp上发生奇怪的事情:
I am trying to solve for the motion of a non linear pendulum. With odeint, everything works like a charm, on solve_ivp hoever something weird happens:
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.
It is a numerical problem. The default relative and absolute tolerances of solve_ivp
are 1e-3 and 1e-6, respectively. For many problems, these values are too low. The default relative tolerance for odeint
is 1.49e-8.
如果将参数rtol=1e-8
添加到solve_ivp
调用中,则表示一致:
If you add the argument rtol=1e-8
to the solve_ivp
call, the plots agree:
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()
情节:
这篇关于使用scipy的solve_ivp解决非线性摆运动的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!