物理引擎的数值积分有哪些好的算法？

Question

我已经在网上寻找了一段时间的物理引擎的集成方法，我正在尝试编写代码以获取乐趣（一定喜欢那里的书呆子：P）。我找到了欧拉方法、RK4 和 Verlet（以及时间校正版本）。我也一直在尝试想出一些自己的方法。我想知道您是否知道其他任何您觉得直观或有帮助的人。谢谢。

编辑：感谢您迄今为止的所有帮助。至于澄清：也许我的意思是数字积分。令人惊讶的是，在我所有的研究中，我还没有找到我想要做的事情的技术名称！也许描述我的具体问题会让我的问题更清楚。假设我想模拟一个球穿过圆形（或实现 3d 后的球形）引力场。该球将遇到力矢量，该力矢量可用于计算球在该特定刻度上所在点的相应加速度矢量。从你的物理课上，你知道速度 = 加速度 * 时间，但我的问题是，从技术上讲，球在该点上只停留了一瞬间，在微积分中用 dt 表示。显然，我不能在 C++ 中使用无穷小的数字，所以我必须使用瞬时积分方法（我在一些阅读中听到的术语，但我可能完全错误）或你认为所谓的数值积分（你可能是是的，所以我改变了标题）。

这是我（成功的）实现数值积分的欧拉方法的尝试：

    //For console output. Note: I know I could just put "using namespace std;" but I hate doing that.
    #include <iostream>
    using std::cout;
    using std::system;
    using std::endl;

    //Program entry
    int main (void)
    {
        //Variable decleration;
        double time = 0;
        double position = 0;
        double velocity = 0;
        double acceleration = 2;
        double dt = 0.000001; //Here is the "instantanious" change in time I was talking about.
        double count = 0; //I use count to make sure I am only displaying the data at whole numbers.

        //Each irritation of this loop is one tick
        while (true)
        {

            //This next bit is a simplified form of Euler's method. It is what I want to "upgrade"
            velocity += acceleration * dt;
            position += velocity * dt;

            if (count == 1/dt) //"count == 1/dt" will only return true if time is a whole number.
            {

                //Simple output to console
                cout << "Time: " << time << endl;
                cout << "Position: " << position << endl;
                cout << "------------------" << endl;
                system ("pause");

                count = 0; //To reset the counter.

            }

            //Update the counters "count" and "time"
            count++;
            time += dt;

        }
        return 1; //Program exit
    }

因为加速度是恒定的，并且这个微分实际上是可解的（为什么我用它来测试，解是位置=时间^ 2，这是相当准确的，但是如果你让它变得更复杂，例如，让加速度随时间变化，算法会很快失去准确性，再次感谢！

Answer 1

您有一个二阶微分方程 (ODE) x''=f(x,x',t)。 x 可以是向量，x' 和 x'' 是相对于时间的一阶和二阶导数。在您的情况下，x 是位置，x' 是速度，x'' 是加速度。通常，通过引入 X=x,Y=x' 将二阶 ODE 转换为一阶 ODE，然后得到

X'=Y
Y'=f(X,Y)

然后您可以使用经典方案来求解 ODE，例如 Runge-Kutta、Dormand-Prince、Adams-Bashforth...

其中许多方法都在 odeint 非常容易使用。

Answer 2

有许多不同的 ODE 数值积分算法。请参阅这篇维基百科文章了解概述。哪种算法合适很大程度上取决于您要求解的 ODE 的属性。欧拉方法很少表现良好，因为您通常需要非常小的步长来获得对解决方案的良好近似（但这非常容易实现，因此可能适合第一次尝试）。还有一些变体，例如向后欧拉方法，可以做得更好一些。

龙格-库塔方法是一类广泛的方法，其中包括欧拉方法。随着方法阶数的增加，通常需要更少的时间步长来达到相同的精度，但在每个时间步长执行计算的成本会变得越来越高 - RK4 经常使用，因为它往往能达到良好的平衡。

在实践中，通常使用自适应步长技术来控制时间步长以达到所需的精度。

ODE 求解器有许多现有的实现，人们为此投入了大量的工作 - 虽然您有兴趣了解它们的工作原理是件好事，但这些求解器可能会变得相当复杂，因此如果您对这些求解器不满意如果您从自己的尝试中获得结果，那么研究现有例程可能是一个更好的主意，例如 GNU 科学图书馆中的内容。