3. 数值积分方法基础

非线性振子方程通常写成二阶形式 \(\ddot{x}+f(x,\dot{x})=g(t)\)，数值模拟的第一步是把它化为一阶方程组。以 Duffing 振子 \(\ddot{x}+\delta\dot{x}+\alpha x+\beta x^3=\gamma\cos\omega t\) 为例，令 \(x_1=x\)、\(x_2=\dot{x}\)，得到

dx1/dt = x2
dx2/dt = -δ*x2 - α*x1 - β*x1^3 + γ*cos(ω*t)

此后所有积分方法都作用在这个状态向量 \(\mathbf{y}=[x_1,x_2]^T\) 上。记步长为 \(h\)，第 \(n\) 步状态为 \(\mathbf{y}_n\)，则各方法的核心差别在于“用哪些时刻的斜率来推进”。

三种方法的推导与实现

显式欧拉法只用当前点的斜率：\(\mathbf{y}_{n+1}=\mathbf{y}_n+h\mathbf{f}(t_n,\mathbf{y}_n)\)。它是一阶方法，实现只需一次函数求值，但误差随步长线性下降，且对刚性或大阻尼系统容易发散。

改进欧拉法（Heun 法）先做一次预测 \(\tilde{\mathbf{y}}=\mathbf{y}_n+h\mathbf{f}(t_n,\mathbf{y}_n)\)，再用起点与预测点的斜率平均：

k1 = f(t_n, y_n)
k2 = f(t_n + h, y_n + h*k1)
y_{n+1} = y_n + h*(k1 + k2)/2

它是二阶方法，每步两次求值，精度明显优于显式欧拉，且仍保持显式结构，适合教学与快速原型。

经典四阶龙格-库塔法（RK4）在区间内取四个斜率加权平均：

k1 = f(t_n, y_n)
k2 = f(t_n + h/2, y_n + h*k1/2)
k3 = f(t_n + h/2, y_n + h*k2/2)
k4 = f(t_n + h,   y_n + h*k3)
y_{n+1} = y_n + h*(k1 + 2*k2 + 2*k3 + k4)/6

其局部截断误差为 \(O(h^5)\)，全局误差为 \(O(h^4)\)，是工程中最常用的通用积分器。下面给出可直接运行的 Python 实现（示例参数为教学假设值，非实测结果）：

import numpy as np

def f(t, y, delta=0.2, alpha=-1.0, beta=1.0, gamma=0.3, omega=1.2):
    x1, x2 = y
    return np.array([x2,
                     -delta*x2 - alpha*x1 - beta*x1**3 + gamma*np.cos(omega*t)])

def rk4(y0, t0, t1, h):
    t, y = t0, np.array(y0, float)
    ts, ys = [t], [y.copy()]
    while t < t1 - 1e-12:
        k1 = f(t, y)
        k2 = f(t + h/2, y + h*k1/2)
        k3 = f(t + h/2, y + h*k2/2)
        k4 = f(t + h,   y + h*k3)
        y = y + h*(k1 + 2*k2 + 2*k3 + k4)/6
        t += h
        ts.append(t); ys.append(y.copy())
    return np.array(ts), np.array(ys)

局部截断误差与全局误差

局部截断误差衡量“假设第 \(n\) 步精确，走一步产生的误差”，对 \(p\) 阶方法为 \(O(h^{p+1})\)；全局误差是累积到终点后的总偏差，为 \(O(h^p)\)。实践中的判断方式是步长减半检验：若步长从 \(h\) 减到 \(h/2\)，全局误差约缩小 \(2^p\) 倍。欧拉法约 2 倍、改进欧拉约 4 倍、RK4 约 16 倍，这一比例可直接用来验证实现是否正确。

稳定性区域概念

把方法作用于测试方程 \(\dot{y}=\lambda y\)（\(\lambda\) 为复数），得到放大因子 \(R(z)\)，其中 \(z=\lambda h\)。方法稳定的条件是 \(|R(z)|\le 1\)，满足该条件的 \(z\) 在复平面上的集合即绝对稳定区域。显式欧拉的 \(R=1+z\)，稳定域是以 \(-1\) 为圆心的单位圆；RK4 的稳定域在实轴上约为 \([-2.785,0]\)，并包含虚轴附近一小段，因此对含阻尼的振子方程通常要求 \(h\) 落在该范围内。对非线性振子，可先用当前状态的雅可比矩阵特征值估计 \(\lambda\)，再据此选步长。

实际实现时，建议先用改进欧拉法做快速验证，再用 RK4 做正式计算；步长选择可先按 \(h\approx T/100\)（\(T\) 为激励周期）试算，然后做步长减半检验，若两次结果在关注量（如振幅、相图闭合程度）上的差异小于允许阈值，即可认为步长足够。对长时间仿真，还需关注能量漂移：无阻尼振子在 RK4 下能量应近似守恒，若出现单调漂移，说明步长偏大或方程右端实现有误。