GSL 求解钟摆运动的 ODE

Question

我正在尝试求解一个钟摆运动的微分方程，给定钟摆初始角度 (x)、重力加速度 (g)、线长 (l) 和时间步长 (h)。我用欧拉法试过这个，一切都很好。但现在我要使用 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))

这是代码 sn-p，如果这还不够，我可以发布整个程序（不会太长），但问题可能出在某处：

    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 在评论中指出的那样，钟摆应该在没有摩擦的情况下继续前进。这是意料之中的。

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

您可以通过绘制系统中不同积分方案的总能量（势能 + 动能）来轻松观察到这一点。

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