非线性单摆方程的闭式解Closed-form solution to the nonlinear pendulum equation

前一篇帖子探讨了非线性单摆方程，并分析了对方程进行线性化后对解的影响。

如果初始位移足够小，可以直接用 θ 替换 sin θ；若初始位移较大，则可通过求解线性化方程并调整周期来显著提高精度。

也可以求得精确解，但无法用初等函数表示，必须使用雅可比椭圆函数。这些函数与三角函数有一定类比关系，但试图明确其对应关系并无帮助。例如，雅可比函数 sn 在某些方面类似于正弦函数，但在其他方面差异很大，具体取决于自变量的取值范围。

我们从微分方程出发：

θ″(t) + c² sin( θ(t) ) = 0

其中 c² = g/L，即重力常数除以摆长，初始条件为 θ(0) = θ₀ 且 θ′(0) = 0，并假设 −π < θ₀ < π。

解为：

θ(t) = 2 arcsin( a cd(ct | m ) )

其中 a = sin(θ₀/2)，m = a²，cd 是十二种雅可比椭圆函数之一。注意，cd（与其他雅可比函数一样）具有一个自变量和一个参数。在上式中，自变量是 ct，参数是 m。

前一篇帖子的最后一个图具有误导性，它大致展示了数值求解微分方程时真实差异与误差各占一半的情况。以下是用于求解非线性方程的代码。

from scipy.special import ellipj, ellipk
from numpy import sin, cos, pi, linspace, arcsin
from scipy.integrate import solve_ivp

def exact_period(θ):
    return 2*ellipk(sin(θ/2)**2)/pi

def nonlinear_ode(t, z):
    x, y = z
    return [y, -sin(x)]    

theta0 = pi/3
b = 2*pi*exact_period(theta0)
t = linspace(0, 2*b, 2000)

sol = solve_ivp(nonlinear_ode, [0, 2*b], [theta0, 0], t_eval=t)

解包含在 sol.y[0] 中。

让我们将数值解与精确解进行比较。

def f(t, c, theta0):
    a = sin(theta0/2)
    m = a**2
    sn, cn, dn, ph = ellipj(c*t, m)
    return 2*arcsin(a*cn/dn)

关于这段代码有两点需要注意：首先，SciPy 没有实现 cd 函数，但可以表示为 cn/dn；其次，ellipj 函数一次返回四个函数，因为计算全部四个所需时间与单独计算其中一个相当。

这是求解微分方程时的误差图。

这是非线性单摆方程的精确解与线性方程拉伸解之间的差异图。