使用scipy的solve_ivp求解非线性摆运动

Question

我仍在尝试了解 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()

我错过了什么？

Answer 1

这是一道数值题。 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()

剧情：

使用scipy的solve_ivp求解非线性摆运动

Using scipy's solve_ivp to solve non linear pendulum motion

scipy

ode

python-3.x