c - GSL 求解钟摆运动的 ODE

Question

我正在尝试求解钟摆运动的微分方程，给定钟摆初始角度 (x)、重力加速度 (g)、线长 (l) 和时间步长 (h)。我已经使用 Euler 方法尝试了这个，一切正常。但是现在我要使用在 GSL 中实现的 Runge-Kutta 方法。我试图通过学习 gsl 手册来实现它，但我遇到了一个问题。钟摆不想停下来。假设我以 1 rad 的初始角度开始，它的峰值倾斜始终为 1 rad，无论它已经摆动了多少次。这是我用来将其提供给 GSL 的等式和函数:

x''(t) + g/l*sin(x(t)) = 0

改造它:

x''(t) = -g/l*sin(x(t))

和分解:

y(t) = x'(t)
y'(t) = -g/l*sin(x(t))

这是代码片段，如果这还不够，我可以发布整个程序(不会太长)，但也许问题出在这里:

    int func (double t, const double x[], double dxdt[], void *params){
        double l = *(double*) params;
        double g = *(double*) (params+sizeof(double));
        dxdt[0] = x[1];
        dxdt[1] = -g/l*sin(x[0]);
        return GSL_SUCCESS;
    }

参数 g 和 l 已正确传递给函数，我已经检查过了。

Answer 1

正如 Barton Chittenden 在评论中指出的那样，钟摆应该在没有摩擦的情况下继续运转。这是预期的。

至于为什么用欧拉法会变慢甚至停止，这涉及到一个微妙而有趣的话题。一个(理想的、无摩擦的)物理摆具有系统中的能量守恒的特性。不同的集成方案在不同程度上保留了该属性。通过一些集成方案，系统中的能量将增长，钟摆将逐渐摆动得更高。对于其他人，能量会丢失，钟摆会停止。这些发生的速度部分取决于方法的顺序；更准确的方法通常会更慢地损失能量。

您可以通过绘制系统中不同集成方案的总能量(势能 + 动能)来轻松观察到这一点。

最后，积分方法还有一个令人着迷的子领域，它保留了像这样的系统的某些守恒量，称为辛方法。