【发布时间】：2015-09-24 15:15:02
【问题描述】：
我有一个一般性问题，我将在更具体的情况下提出这个问题。
如果想找到双摆的动力学，可以用数学方法推导出运动方程，将 ODE 重写为对数值计算有用的特殊形式，然后使用 C++ 中的 odeint 求解 ODE（参见堆栈溢出的示例https://stackoverflow.com/a/30582741)。
现在假设我们想对 n 个耦合摆（n 在运行时已知）做同样的事情。这需要我们写出所谓的拉格朗日（动能-势能），这个函数的不同导数将是我们需要求解的 ODE。此外，这些 ODE 必须以适合 odeint 的形式重写。这对于一般 n 来说很难手工完成。
在像 Mathematica 和 Maple 这样的程序中，这实际上很容易。可以从拉格朗日函数中显式推导出所需的常微分方程，而常微分方程求解器不需要我们将方程置于任何特殊形式（请参阅此处的数学示例https://mathematica.stackexchange.com/a/84279）。
有没有可能在c++中做这样的事情而不会遇到太多麻烦？
可能的方法：
一种可能的方法是使用 c++ 包ginac。这可以帮助我们分析推导 ODE。但我不知道如何将来自 ginac 的表达式重写为适合 odeint 中数值计算的形式。有什么想法吗？
【问题讨论】：
-
对于仅在运行时知道的 n，恐怕您面临着几乎不可能的困难。如果您在编译时知道 n，我认为您可以基于自动微分和表达式模板来实现这样的计算。从我对 ginac 的（表面）检查来看，我认为这个库也是基于表达式模板的，因此在编译时需要 n。
-
@mariomulansky 非常感谢您的评论。现在，如果我能在编译时做到这一点，我会很高兴。如何使用表达式模板将 ginac 输出转换为 odeint 可以理解的形式？是否存在一些简单的标准方法？
-
我没有使用 GiNaC 的经验，所以帮不上什么忙。然而，通过他们的教程阅读更多内容，我想到你可以在运行时构造一个表达式，区分它并评估结果。所以我上面的第一条评论可能是不正确的。尝试将 Langragian 写成 GiNaC 表达式、区分它并在 rhs 中使用生成的表达式来表示 odeint 声音听起来确实是一种可行的方法。如果你真的让它运行，我会对代码非常感兴趣。或许你可以在 odeint 的 github 上分享：github.com/headmyshoulder/odeint-v2
-
您也可以将 Lagrange 系统视为 index-3 DAE 系统，并使用通常用于约束动力学的 BDF 或 RADAU 求解器。或者手动对约束方程进行两次微分，并通过 ODE 函数中的牛顿对二阶导数的所得非线性系统进行数值求解。
-
@mariomulansky 是的，似乎可以直接评估表达式，这应该足以将其放入 odeint。当我有时间时，我会尝试实现它，我很乐意分享代码。对于这种更具体的情况（n-pendula），可能有一种更简单的方法。如果我使用普通笛卡尔坐标并使用拉格朗日乘数实现摆约束，则哈密顿量的形式为 H=T(p) + V(q)（q 之间的乘数）。也许从那里可以直接使用 odeint 的辛积分器？
标签： c++ simulation physics symbolic-math odeint