c++ - 物理模拟为简单轨迹演算提供了(非常)不准确的位置

我想在游戏中实现一个物理引擎，以便计算施加于它们的力的物体的轨迹。
该引擎将根据对象的先前状态计算对象的每个状态。当然，这意味着要在两个时间单位之间进行大量计算才能足够精确。

为了正确地做到这一点，我首先想知道这种获取位置的方法与运动学方程之间的差异有多大。
所以我制作了这段代码，它存储了由模拟和文件中的方程给出的位置 (x, y, z)。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "header.h"


Body nouveauCorps(Body body, Vector3 force, double deltaT){
    double m = body.mass;
    double t = deltaT;

    //Newton's second law:
    double ax = force.x/m;
    double ay = force.y/m;
    double az = force.z/m;

    body.speedx += ax*t;
    body.speedy += ay*t;
    body.speedz += az*t;

    body.x +=t*body.speedx;
    body.y +=t*body.speedy;
    body.z +=t*body.speedz;

    return body;
}

int main()
{
    //Initial conditions:
    double posX = 1.4568899;
    double posY = 5.6584225;
    double posZ = -8.8944444;
    double speedX = 0.232323;
    double speedY = -1.6565656;
    double speedZ = -8.6565656;
    double mass = 558.74;

    //Force applied:
    Vector3 force = {5.8745554, -97887.568, 543.5875};

    Body body = {posX, posY, posZ, speedX, speedY, speedZ, mass};

    double duration = 10.0;
    double pointsPS = 100.0; //Points Per Second
    double pointsTot = duration * pointsPS;

    char name[20];
    sprintf(name, "BN_%fs-%fpts.txt", duration, pointsPS);

    remove(name);
    FILE* fichier = NULL;
    fichier = fopen(name, "w");

    for(int i=1; i<=pointsTot; i++){
        body = nouveauCorps(body, force, duration/pointsTot);
        double t = i/pointsPS;

        //Make a table: TIME | POS_X, Y, Z by simulation | POS_X, Y, Z by modele (reference)
        fprintf(fichier, "%e \t %e \t %e \t %e \t %e \t %e \t %e\n", t, body.x, body.y, body.z, force.x*(t*t)/2.0/mass + speedX*t + posX, force.y*(t*t)/2.0/mass + speedY*t + posY, force.z*(t*t)/2.0/mass + speedZ*t + posZ);
    }
    return 0;
}

问题是，对于简单的数字(比如在 -9.81 重力场中的简单坠落)，我得到了不错的位置，但是对于更大(且相当随机)的数字，我得到的位置不准确。

那是浮点问题吗？

这是结果，有相对误差。 (注意:标签轴是法语，Temps = Time)。

图表

Black+dashed : 来自运动方程的值

红色:每秒 100 点

橙色:每秒 1000 点

绿色:每秒 10000 点

最佳答案

这不是浮点问题。事实上，即使您使用精确算术，您也会看到同样的问题。

这个错误对于数值积分本身以及您正在使用的特定方法和您正在解决的 ODE 来说非常重要。

在这种情况下，您将使用称为 Forward Euler 的积分方案。这可能是求解一阶 ODE 的最简单方法。当然，这给它留下了一些不受欢迎的特征。

一方面，它在每一步都会引入错误。错误的大小是 O(Δt²) 。这意味着单个时间步长的误差大致与时间步长的平方成正比。因此，如果您将时间步长的大小减半，大致可以将增量误差降低到值的 1/4。

但是由于您减少了时间步长，因此您必须执行更多步骤来模拟相同的时间量。所以你增加了更多但更小的错误。这就是累积误差为 O(Δt) 的原因。因此，在整个模拟时间内，如果您采用一半大的时间步长，则会得到一半的累积误差。

最终，这个累积错误就是你所看到的。您可以在误差图中看到，每次将时间步数增加 10 倍时，最终误差最终会减少约 10 倍:因为时间步长小 10 倍，因此总误差结束小了大约 10 倍。

另一个问题是 Forward Euler 表现出所谓的条件稳定性。这意味着在某些情况下累积误差可能会无限增长。要了解原因，让我们看一个简单的 ODE:

x' = -k * x

其中 k 是某个常数。这个 ODE 的精确解是 x(t) = x(0) * exp( -k * t ) 。所以只要 k 是正数，x 应该会随着时间的增加而趋向于 0。

然而，如果我们尝试使用 Forward Euler 来近似它，我们会得到如下所示的结果:

x(t + Δt) = x(t) + Δt * ( -k * x[n] )
          = ( 1 - k * Δt ) * x(t)

这是一个我们可以解决的简单递推关系:

x(t) = ( 1 - k * Δt )^(t / Δt) * x(0)

现在，随着 t 变大，我们知道精确解为 10 到 0。但是 Forward Euler 解决方案仅在 |1 - k * Δt| < 1 时才这样做。请注意该表达式如何取决于步长以及来自 ODE 的 k 项。如果 k 真的很大，我们需要一个非常非常小的时间步来防止解决方案爆炸。这就是为什么它具有所谓的条件稳定性:解的稳定性取决于时间步长。

还有许多其他问题，但这是一个广泛的主题，我无法在一个答案中涵盖所有内容。

关于c++ - 物理模拟为简单轨迹演算提供了(非常)不准确的位置，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/45200210/