我正在尝试使用四阶 Runge-Kutta 方法求解非线性摆，用于 a=0.0 到 b=110 秒之间的限制，并模拟结果以观察摆的运动。但是当我在 RK4 方法中增加步数 N 时，钟摆并没有衰减，它一直在振荡，但是对于较小的 N 值，我可以观察到衰减效果。

from math import sin, pi
import numpy as np
import matplotlib.pyplot as plt


# constants
G = 9.51        # acceleration due to gravity
L = 0.1         # length of the pendulum
T = 0.5         # time period of the pendulum
ANGLE = 90.0    # initial angle of the pendulum


def f(r, t):
    """Vectorized simultaneous ODEs functions."""
    theta = r[0]
    omega = r[1]
    ftheta = omega
    fomega = - (G / L) * sin(theta)

    return np.array([ftheta, fomega], dtype=float)


a = 0.0
b = 110.0
N = 2000
h = (b - a) / N

tpoints = np.arange(a, b, h)
thetapoints = []
omegapoints = []

theta0 = (pi / 180) * ANGLE
omega0 = (2 * pi) / T
r = np.array([theta0, omega0], dtype=float)
for t in tpoints:
    thetapoints.append(r[0])
    omegapoints.append(r[1])
    k1 = h * f(r, t)
    k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
    k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
    k4 = h * f(r + k3, t + h)
    r += (k1 + 2 * k2 + 2 * k3 + k4) / 6

ypoints = -L * np.cos(thetapoints)
xpoints = L * np.sin(thetapoints)

plt.plot(tpoints, omegapoints, 'k', lw=0.8)
plt.xlabel('t')
plt.ylabel('omega')
plt.savefig('2000_steps.jpeg', dpi=300, format='jpeg')
plt.show()

欧米茄与时间的输出：

N = 1000

N = 2000

我坚信这是由于近似

对于较小的值$\theta$在常微分方程中。