下面是一个弹簧摆的实现：
Spring Pendulum Animation
本文将介绍如何使用拉格朗日方程来描述单摆和弹簧摆的运动，并通过数值模拟展示其动态行为。
拉格朗日力学
拉格朗日力学 (Lagrangian Mechanics) 是一种描述物理系统动力学行为的理论框架。它基于拉格朗日量 (Lagrangian) \(\mathcal{L}\) 的概念，定义为系统的动能 \(T\) 与势能 \(U\) 之差：
\[\mathcal{L} = T - U\]通过拉格朗日量，可以导出描述系统运动的欧拉-拉格朗日方程 (Euler-Lagrange Equation)：
\[\frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{q}_i}\right) - \frac{\partial \mathcal{L}}{\partial q_i} = 0\]其中，\(q_i\) 是系统的广义坐标，\(\dot{q}_i\) 是其对应的广义速度。
欧拉-拉格朗日方程这样的常微分方程 (Ordinary Differential Equation, ODE) 可以通过数值方法进行求解，从而模拟系统的动态行为。
在实际计算中，我们通常需要先表达系统的动能和势能，然后计算拉格朗日量，最后求解欧拉-拉格朗日方程以获得运动方程。
Runge-Kutta 方法
我们考虑常微分方程的一般形式：
\[\frac{dy}{dt} = f\left(t, y\right), \quad y\left(t_0\right) = y_0\]对于该方程，欧拉方法 (Euler Method) 提供了一种简单的数值解法：
\[y_{n+1} = y_n + f\left(t_n, y_n\right) \Delta t\]但欧拉方法的精度较低
- 如果 \(\Delta t\) 较大，容易产生累积误差
- 如果 \(\Delta t\) 较小，计算量会显著增加
- 对于刚性方程，欧拉方法可能不稳定
为提高精度，可以采用更高级的数值积分方法，如 Runge-Kutta 方法 (Runge-Kutta Method)。
相比于欧拉方法使用区间起点 \(\left(t_n, y_n\right)\) 的斜率，Runge-Kutta 方法通过在区间内多个点计算斜率的加权平均来更新 \(y\) 的值，从而提高精度。
Runge-Kutta 方法的经典形式是四阶 Runge-Kutta 方法 (RK4)，其更新公式为：
\[\begin{aligned} k_1 &= f\left(t_n, y_n\right) \\ k_2 &= f\left(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_1\right) \\ k_3 &= f\left(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_2\right) \\ k_4 &= f\left(t_n + \Delta t, y_n + \Delta t k_3\right) \\ y_{n+1} &= y_n + \frac{\Delta t}{6} \left(k_1 + 2 k_2 + 2 k_3 + k_4\right) \end{aligned}\]其中
- \(k_1\) 是区间起点的斜率
- \(k_2\) 是使用 \(k_1\) 估计的区间中点的斜率
- \(k_3\) 是使用 \(k_2\) 估计的区间中点的斜率
- \(k_4\) 是使用 \(k_2\) 估计的区间终点的斜率
最终，根据泰勒展开，得到各个斜率所占的权重，综合计算出 \(y_{n+1}\) 的值。
这里的泰勒展开具体算法如下。
对于 RK4 的数值解 \(y_{n+1} = y_n + \Delta t \left(a_1 k_1 + a_2 k_2 + a_3 k_3 + a_4 k_4\right)\)，我们希望其与解析解的泰勒展开在 \(\Delta t\) 的各阶项上尽可能匹配。
我们有：
- \(k_1 = f\left(t_n, y_n\right)\)
- \(k_2 = f\left(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_1\right)\)
- \(k_3 = f\left(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_2\right)\)
- \(k_4 = f\left(t_n + \Delta t, y_n + \Delta t k_3\right)\)
将展开后的 \(k_1\)、\(k_2\)、\(k_3\)、\(k_4\) 代入 \(y_{n+1}\) 的表达式中，并与解析解的泰勒展开进行比较，可以得到一组关于 \(\Delta t\) 各阶系数的方程。
比较 \(\Delta t\) 的系数，我们有：
- 一阶项：\(a_1 + a_2 + a_3 + a_4 = 1\)
- 二阶项：\(\frac{1}{2} a_2 + \frac{1}{2} a_3 + a_4 = \frac{1}{2}\)
- 三阶项：\(\frac{1}{4} a_2 + \frac{1}{4} a_3 + a_4 = \frac{1}{6}\)
- 四阶项：\(\frac{1}{8} a_2 + \frac{1}{8} a_3 + a_4 = \frac{1}{24}\)
解这组方程，得到：
\[a_1 = \frac{1}{6}, \quad a_2 = \frac{1}{3}, \quad a_3 = \frac{1}{3}, \quad a_4 = \frac{1}{6}\]因此，RK4 方法的更新公式为：
\[y_{n+1} = y_n + \frac{\Delta t}{6} \left(k_1 + 2 k_2 + 2 k_3 + k_4\right)\]
RK4 方法的局部截断误差为 \(O\left(\Delta t^5\right)\)，全局截断误差为 \(O\left(\Delta t^4\right)\)，显著优于欧拉方法的 \(O\left(\Delta t^2\right)\)。
单摆
假设一个质量为 \(m\) 的物体通过一根长度为 \(L\) 的不可伸长且无质量的绳子悬挂在固定点上，形成一个单摆系统。设摆锤与竖直方向的夹角为 \(\theta\)，则摆锤的笛卡尔坐标可以表示为：
\[\begin{aligned} x_m &= L \sin{\theta} \\ y_m &= -L \cos{\theta} \end{aligned}\]求其对于时间的一阶导数：
\[\begin{aligned} \dot{x}_m &= L \dot{\theta}\cos{\theta} \\ \dot{y}_m &= L \dot{\theta}\sin{\theta} \end{aligned}\]动能为：
\[\begin{aligned} T &= \frac{1}{2} m {\nu}^2 \\ &= \frac{1}{2} m \left(\dot{x}_m^2 + \dot{y}_m^2\right) \\ &= \frac{1}{2} m \left(L^2 \dot{\theta}^2 \cos^2{\theta} + L^2 \dot{\theta}^2 \sin^2{\theta}\right) \\ &= \frac{1}{2} m L^2 \dot{\theta}^2 \end{aligned}\]势能为：
\[\begin{aligned} U &= m g h \\ &= m g \left(-L \cos{\theta}\right) \\ &= -m g L \cos{\theta} \end{aligned}\]由此计算出拉格朗日量：
\[\begin{aligned} \mathcal{L} &= T - U \\ &= \frac{1}{2} m L^2 \dot{\theta}^2 -\left(- m g L \cos{\theta}\right) \\ &= \frac{1}{2} m L^2 \dot{\theta}^2 + m g L \cos{\theta} \\ &= m L \left(\frac{1}{2} L \dot{\theta}^2 + g \cos{\theta}\right) \end{aligned}\]求解欧拉-拉格朗日微分方程1：
\[\begin{aligned} \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}}\right) - \frac{\partial \mathcal{L}}{\partial \theta} &= 0 \\ \frac{d}{dt} \left(\frac{\partial }{\partial \dot{\theta}} \left[m L \left(\frac{1}{2} L \dot{\theta}^2 + g \cos{\theta}\right)\right]\right) \\ - \frac{\partial }{\partial \theta} \left[m L \left(\frac{1}{2} L \dot{\theta}^2 + g \cos{\theta}\right)\right] &= 0 \\ \frac{d}{dt} \left[\frac{\partial }{\partial \dot{\theta}} \left(\frac{1}{2} L \dot{\theta}^2 + g \cos{\theta}\right)\right] \\ - \frac{\partial }{\partial \theta} \left(\frac{1}{2} L \dot{\theta}^2 + g \cos{\theta}\right) &= 0 \\ \frac{d}{dt} \left[\frac{\partial }{\partial \dot{\theta}} \left(\frac{1}{2} L \dot{\theta}^2\right) + \frac{\partial }{\partial \dot{\theta}} \left(g \cos{\theta}\right)\right] \\ - \left[\frac{\partial }{\partial \theta} \left(\frac{1}{2} L \dot{\theta}^2\right) + \frac{\partial }{\partial \theta} \left(g \cos{\theta}\right)\right] &= 0 \\ \frac{d}{dt} \left(L \dot{\theta} + 0\right) - \left[0 + \left(- g \sin{\theta}\right)\right] &= 0 \\ \frac{d}{dt} \left(L \dot{\theta}\right) + g \sin{\theta} &= 0 \\ L \ddot{\theta} + g \sin{\theta} &= 0 \end{aligned}\] \[\ddot{\theta} = -\frac{g \sin{\theta}}{L}\]给定初始条件 \(\theta(0) = \theta_0\)，\(\dot{\theta}(0) = \omega_0\)，可以利用上面得到的 \(\ddot{\theta}\) 方程进行数值模拟，迭代计算 \(\theta(t)\) 和 \(\dot{\theta}(t)\) 的值：
let theta = theta0; // 初始角度
let omega = omega0; // 初始角速度
let dt = 0.01; // 时间步长
function simulatePendulum() {
let alpha = - g * Math.sin(theta) / L; // 计算角加速度
omega += alpha * dt; // 更新角速度
theta += omega * dt; // 更新角度
}
while (true) {
simulatePendulum(); // 持续模拟
// 在这里可以添加代码来绘制或记录 theta 和 omega 的值
}
但这种方法在角度较大时会产生较大误差，可以使用更精确的数值积分方法如 Runge-Kutta 方法来提高模拟精度。
Runge-Kutta 方法利用四阶近似来更新状态变量：
\[\begin{aligned} k_1^{\theta} &= \omega \\ k_1^{\omega} &= -\frac{g}{L} \sin{\theta} \\ k_2^{\theta} &= \omega + 0.5 k_1^{\omega} dt \\ k_2^{\omega} &= -\frac{g}{L} \sin{\left(\theta + 0.5 k_1^{\theta} dt\right)} \\ k_3^{\theta} &= \omega + 0.5 k_2^{\omega} dt \\ k_3^{\omega} &= -\frac{g}{L} \sin{\left(\theta + 0.5 k_2^{\theta} dt\right)} \\ k_4^{\theta} &= \omega + k_3^{\omega} dt \\ k_4^{\omega} &= -\frac{g}{L} \sin{\left(\theta + k_3^{\theta} dt\right)} \\ \theta &= \theta + \frac{1}{6} \left(k_1^{\theta} + 2 k_2^{\theta} + 2 k_3^{\theta} + k_4^{\theta}\right) \\ \omega &= \omega + \frac{1}{6} \left(k_1^{\omega} + 2 k_2^{\omega} + 2 k_3^{\omega} + k_4^{\omega}\right) \end{aligned}\]let theta = theta0; // 初始角度
let omega = omega0; // 初始角速度
let dt = 0.01; // 时间步长
function simulatePendulum() {
let k1_theta = omega;
let k1_omega = - (g / L) * Math.sin(theta);
let k2_theta = omega + 0.5 * k1_omega * dt;
let k2_omega = - (g / L) * Math.sin(theta + 0.5 * k1_theta * dt);
let k3_theta = omega + 0.5 * k2_omega * dt;
let k3_omega = - (g / L) * Math.sin(theta + 0.5 * k2_theta * dt);
let k4_theta = omega + k3_omega * dt;
let k4_omega = - (g / L) * Math.sin(theta + k3_theta * dt);
theta += (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta) / 6;
omega += (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega) / 6;
}
while (true) {
simulatePendulum(); // 持续模拟
// 在这里可以添加代码来绘制或记录 theta 和 omega 的值
}
下面是完整的代码示例：
Pendulum Simulation
弹簧摆
假设一个质量为 \(m\) 的物体通过一根弹簧连接在固定点上，弹簧的自然长度为 \(L_0\)，弹簧常数为 \(k\)。设摆锤与竖直方向的夹角为 \(\theta\)，弹簧的伸长量为 \(L^\prime\)，则摆锤的笛卡尔坐标可以表示为：
\[\begin{aligned} x_m &= (L_0 + L^\prime) \sin{\theta} \\ y_m &= -(L_0 + L^\prime) \cos{\theta} \end{aligned}\]求其对于时间的一阶导数：
\[\begin{aligned} \dot{x}_m &= \dot{L}^\prime \sin{\theta} + \left(L_0 + L^\prime\right) \dot{\theta} \cos{\theta} \\ \dot{y}_m &= -\dot{L}^\prime \cos{\theta} + \left(L_0 + L^\prime\right) \dot{\theta} \sin{\theta} \end{aligned}\]动能为：
\[\begin{aligned} T =& \frac{1}{2} m {\nu}^2 \\ =& \frac{1}{2} m \left(\dot{x}_m^2 + \dot{y}_m^2\right) \\ =& \frac{1}{2} m \\ & \left({\left[\dot{L}^\prime \sin{\theta} + \left(L_0 + L^\prime\right) \dot{\theta} \cos{\theta}\right]}^2 + {\left[-\dot{L}^\prime \cos{\theta} + \left(L_0 + L^\prime\right) \dot{\theta} \sin{\theta}\right]}^2\right) \\ =& \frac{1}{2} m \left[\dot{L}^{\prime 2} \left(\sin^2{\theta} + \cos^2{\theta}\right) + \left(L_0 + L^\prime\right)^2 \dot{\theta}^2 \left(\cos^2{\theta} + \sin^2{\theta}\right)\right] \\ =& \frac{1}{2} m \left[\dot{L}^{\prime 2} + \left(L_0 + L^\prime\right)^2 \dot{\theta}^2\right] \end{aligned}\]势能为：
\[\begin{aligned} U =& U_{\text{gravity}} + U_{\text{spring}} \\ =& m g h + \frac{1}{2} k L^{\prime 2} \\ =& m g \left[-\left(L_0 + L^\prime\right) \cos{\theta}\right] + \frac{1}{2} k L^{\prime 2} \\ =& -m g \left(L_0 + L^\prime\right) \cos{\theta} + \frac{1}{2} k L^{\prime 2} \end{aligned}\]由此计算出拉格朗日量：
\[\begin{aligned} \mathcal{L} =& T - U \\ =& \frac{1}{2} m \left[\dot{L}^{\prime 2} + \left(L_0 + L^\prime\right)^2 \dot{\theta}^2\right] \\ & - \left[-m g \left(L_0 + L^\prime\right) \cos{\theta} + \frac{1}{2} k L^{\prime 2}\right] \\ =& \frac{1}{2} m \left[\dot{L}^{\prime 2} + \left(L_0 + L^\prime\right)^2 \dot{\theta}^2\right] \\ & + m g \left(L_0 + L^\prime\right) \cos{\theta} - \frac{1}{2} k L^{\prime 2} \end{aligned}\]求解欧拉-拉格朗日微分方程：
\[\begin{aligned} \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot{L}^\prime}\right) - \frac{\partial \mathcal{L}}{\partial L^\prime} &= 0 \\ \frac{d}{dt} \left(\frac{\partial }{\partial \dot{L}^\prime} \left[\frac{1}{2} m \left(\dot{L}^{\prime 2} + \left(L_0 + L^\prime\right)^2 \dot{\theta}^2\right) + m g \left(L_0 + L^\prime\right) \cos{\theta} - \frac{1}{2} k L^{\prime 2}\right]\right) \\ - \frac{\partial }{\partial L^\prime} \left[\frac{1}{2} m \left(\dot{L}^{\prime 2} + \left(L_0 + L^\prime\right)^2 \dot{\theta}^2\right) + m g \left(L_0 + L^\prime\right) \cos{\theta} - \frac{1}{2} k L^{\prime 2}\right] &= 0 \\ \frac{d}{dt} \left[m \dot{L}^\prime\right] - \left[m \left(L_0 + L^\prime\right) \dot{\theta}^2 + m g \cos{\theta} - k L^\prime\right] &= 0 \\ m \ddot{L}^\prime - m \left(L_0 + L^\prime\right) \dot{\theta}^2 - m g \cos{\theta} + k L^\prime &= 0 \end{aligned}\] \[\begin{aligned} \ddot{L}^\prime &= \left(L_0 + L^\prime\right) \dot{\theta}^2 + g \cos{\theta} - \frac{k}{m} L^\prime \\ \ddot{\theta} &= -\frac{2 \dot{L}^\prime \dot{\theta} + g \sin{\theta}}{L_0 + L^\prime} \end{aligned}\]给定初始条件 \(L^\prime(0) = L_0^\prime\)，\(\dot{L}^\prime(0) = v_0\)，\(\theta(0) = \theta_0\)，\(\dot{\theta}(0) = \omega_0\)，可以利用上面得到的 \(\ddot{L}^\prime\) 和 \(\ddot{\theta}\) 方程进行数值模拟，迭代计算 \(L^\prime(t)\)、\(\dot{L}^\prime(t)\)、\(\theta(t)\) 和 \(\dot{\theta}(t)\) 的值：
let L0 = 100; // 弹簧自然长度
let k = 0.5; // 弹簧常数
let m = 1; // 质量
let g = 9.81; // 重力加速度
let L_prime = L0_prime; // 初始弹簧伸长量
let v = v0; // 初始伸长速度
let theta = theta0; // 初始角度
let omega = omega0; // 初始角速度
let dt = 0.01; // 时间步长
function simulateSpringPendulum() {
let alpha_L = (L0 + L_prime) * omega * omega + g * Math.cos(theta) - (k / m) * L_prime; // 计算伸长加速度
let alpha_theta = - (2 * v * omega + g * Math.sin(theta)) / (L0 + L_prime); // 计算角加速度
v += alpha_L * dt; // 更新伸长速度
L_prime += v * dt; // 更新伸长量
omega += alpha_theta * dt; // 更新角速度
theta += omega * dt; // 更新角度
}
while (true) {
simulateSpringPendulum(); // 持续模拟
// 在这里可以添加代码来绘制或记录 L_prime、v、theta 和 omega 的值
}
同样的，可以使用 Runge-Kutta 方法来提高模拟精度。
\[\begin{aligned} k_1^{L} =& v \\ k_1^{v} =& \left(L_0 + L^\prime\right) \omega^2 + g \cos{\theta} - \frac{k}{m} L^\prime \\ k_1^{\theta} =& \omega \\ k_1^{\omega} =& -\frac{2 v \omega + g \sin{\theta}}{L_0 + L^\prime} \\ k_2^{L} =& v + 0.5 k_1^{v} dt \\ k_2^{v} =& \left(L_0 + L^\prime + 0.5 k_1^{L} dt\right) \left(\omega + 0.5 k_1^{\omega} dt\right)^2 \\ & + g \cos{\left(\theta + 0.5 k_1^{\theta} dt\right)} - \frac{k}{m} \left(L^\prime + 0.5 k_1^{L} dt\right) \\ k_2^{\theta} =& \omega + 0.5 k_1^{\omega} dt \\ k_2^{\omega} =& -\frac{2 \left(v + 0.5 k_1^{v} dt\right) \left(\omega + 0.5 k_1^{\omega} dt\right) + g \sin{\left(\theta + 0.5 k_1^{\theta} dt\right)}}{L_0 + L^\prime + 0.5 k_1^{L} dt} \\ k_3^{L} =& v + 0.5 k_2^{v} dt \\ k_3^{v} =& \left(L_0 + L^\prime + 0.5 k_2^{L} dt\right) \left(\omega + 0.5 k_2^{\omega} dt\right)^2 \\ & + g \cos{\left(\theta + 0.5 k_2^{\theta} dt\right)} - \frac{k}{m} \left(L^\prime + 0.5 k_2^{L} dt\right) \\ k_3^{\theta} =& \omega + 0.5 k_2^{\omega} dt \\ k_3^{\omega} =& -\frac{2 \left(v + 0.5 k_2^{v} dt\right) \left(\omega + 0.5 k_2^{\omega} dt\right) + g \sin{\left(\theta + 0.5 k_2^{\theta} dt\right)}}{L_0 + L^\prime + 0.5 k_2^{L} dt} \\ k_4^{L} =& v + k_3^{v} dt \\ k_4^{v} =& \left(L_0 + L^\prime + k_3^{L} dt\right) \left(\omega + k_3^{\omega} dt\right)^2 \\ & + g \cos{\left(\theta + k_3^{\theta} dt\right)} - \frac{k}{m} \left(L^\prime + k_3^{L} dt\right) \\ k_4^{\theta} =& \omega + k_3^{\omega} dt \\ k_4^{\omega} =& -\frac{2 \left(v + k_3^{v} dt\right) \left(\omega + k_3^{\omega} dt\right) + g \sin{\left(\theta + k_3^{\theta} dt\right)}}{L_0 + L^\prime + k_3^{L} dt} \\ L^\prime =& L^\prime + \frac{1}{6} \left(k_1^{L} + 2 k_2^{L} + 2 k_3^{L} + k_4^{L}\right) \\ v =& v + \frac{1}{6} \left(k_1^{v} + 2 k_2^{v} + 2 k_3^{v} + k_4^{v}\right) \\ \theta =& \theta + \frac{1}{6} \left(k_1^{\theta} + 2 k_2^{\theta} + 2 k_3^{\theta} + k_4^{\theta}\right) \\ \omega =& \omega + \frac{1}{6} \left(k_1^{\omega} + 2 k_2^{\omega} + 2 k_3^{\omega} + k_4^{\omega}\right) \end{aligned}\]let L0 = 100; // 弹簧自然长度
let k = 0.5; // 弹簧常数
let m = 1; // 质量
let g = 9.81; // 重力加速度
let L_prime = L0_prime; // 初始弹簧伸长量
let v = v0; // 初始伸长速度
let theta = theta0; // 初始角度
let omega = omega0; // 初始角速度
let dt = 0.01; // 时间步长
function simulateSpringPendulum() {
const k1_L = v;
const k1_v = (L0 + L_prime) * omega * omega + g * Math.cos(theta) - (k / m) * L_prime;
const k1_theta = omega;
const k1_omega = - (2 * v * omega + g * Math.sin(theta)) / (L0 + L_prime);
const k2_L = v + 0.5 * k1_v * dt;
const k2_v = (L0 + L_prime + 0.5 * k1_L * dt) * (omega + 0.5 * k1_omega * dt) ** 2
+ g * Math.cos(theta + 0.5 * k1_theta * dt)
- (k / m) * (L_prime + 0.5 * k1_L * dt);
const k2_theta = omega + 0.5 * k1_omega * dt;
const k2_omega = - (2 * (v + 0.5 * k1_v * dt) * (omega + 0.5 * k1_omega * dt)
+ g * Math.sin(theta + 0.5 * k1_theta * dt))
/ (L0 + L_prime + 0.5 * k1_L * dt);
const k3_L = v + 0.5 * k2_v * dt;
const k3_v = (L0 + L_prime + 0.5 * k2_L * dt) * (omega + 0.5 * k2_omega * dt) ** 2
+ g * Math.cos(theta + 0.5 * k2_theta * dt)
- (k / m) * (L_prime + 0.5 * k2_L * dt);
const k3_theta = omega + 0.5 * k2_omega * dt;
const k3_omega = - (2 * (v + 0.5 * k2_v * dt) * (omega + 0.5 * k2_omega * dt)
+ g * Math.sin(theta + 0.5 * k2_theta * dt))
/ (L0 + L_prime + 0.5 * k2_L * dt);
const k4_L = v + k3_v * dt;
const k4_v = (L0 + L_prime + k3_L * dt) * (omega + k3_omega * dt) ** 2
+ g * Math.cos(theta + k3_theta * dt)
- (k / m) * (L_prime + k3_L * dt);
const k4_theta = omega + k3_omega * dt;
const k4_omega = - (2 * (v + k3_v * dt) * (omega + k3_omega * dt)
+ g * Math.sin(theta + k3_theta * dt))
/ (L0 + L_prime + k3_L * dt);
L_prime += (k1_L + 2 * k2_L + 2 * k3_L + k4_L) / 6 * dt;
v += (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6 * dt;
theta += (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta) / 6 * dt;
omega += (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega) / 6 * dt;
}
while (true) {
simulateSpringPendulum(); // 持续模拟
// 在这里可以添加代码来绘制或记录 L_prime、v、theta 和 omega 的值
}
在实现前，我们还需要注意弹簧绳只能产生拉力，不能产生压力，因此当弹簧伸长量 \(L^\prime < 0\) 时，弹力应为 \(0\)，此时物体将只受重力作用。
下面是完整的代码示例：
Spring Pendulum Simulation
我们这边使用绳子的颜色来表示弹簧的紧绷程度。
你可以试着将初始角度调大一些，然后观察弹簧摆在弹簧绳有松弛时的运动情况。
这段代码似乎还有点 bug，小球如果超出画布范围就会卡住，可以尝试修复它。
讨论
在上面的内容中，我们看到似乎只能使用数值方法来模拟弹簧摆的运动轨迹，是否无法得到解析解呢？
事实上，弹簧摆的运动方程是一个耦合的非线性微分方程组，一般情况下是无法通过解析方法求解的。不过，在某些特定条件下，可以对方程进行线性化处理，从而得到近似的解析解。例如，当摆动角度较小且弹簧伸长量较小时，可以将 \(\sin{\theta} \approx \theta\) 和 \(L^\prime \approx 0\) 代入方程，得到简化后的线性微分方程组，从而求解出近似的解析解。
但这种近似解仅在小角度和小伸长量范围内有效，对于大角度或大伸长量的情况，仍然需要依赖数值方法进行模拟。
对于数值解法的选择，除了上面介绍的 Runge-Kutta 方法外，还可以考虑其他数值积分方法，如 Verlet 积分法、Leapfrog 方法等。这些方法在处理物理系统时具有较好的能量守恒性质，适合用于长时间模拟。
不过，即便选用了更高级的数值方法，我们最终的结果仍然可能出现问题，这是由于浮点数运算的精度限制所导致的。例如，对于 JavaScript 有一个著名的浮点数精度问题 \(0.1 + 0.2 = 0.30000000000000004\)：
Ch3nyang's Online Code Runner
这样的精度误差在长时间的数值模拟中会逐渐积累，最终导致结果偏离真实情况。因此，在进行数值模拟时，可能还需要配合一些数值稳定化技术，如自适应时间步长调整、能量修正等，以减小误差的积累。
Comments