物理模拟在计算简单轨迹的位置时给出（非常）不准确的结果

Question

物理模拟在计算简单轨迹的位置时给出（非常）不准确的结果

4

我希望在游戏中实现物理引擎，以计算施加力后的物体运动轨迹。该引擎将根据物体的先前状态计算每个状态。当然，这意味着需要在两个时间单位之间进行大量计算，以保证足够精确。

为了更好地做到这一点，我首先想知道使用此方法获取位置与使用运动学方程之间的差异有多大。因此，我编写了此代码，将模拟和方程所给出的位置（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 = 时间）。图表如下：

黑色+虚线：根据运动方程计算的值
红色：每秒100个数据点
橙色：每秒1000个数据点
绿色：每秒10000个数据点

- T0T0R

1

请参考以下链接： https://dev59.com/G3RB5IYBdhLWcg3wj36c https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html - Jesper Juhl

2

我认为这更像是一个模拟问题，而不是浮点数问题。例如，如果您使用pos+=v_old*t+a/2*t^2公式，可以改善模拟。另一个提示是减小deltaT。如果模拟变得更好，那么肯定是一个模拟问题。 - geza

1

您可能还想更改您正在使用的数值方法。 - Bob__

如果你想要一个好的演示，可以尝试计算一个沿着圆形路径运动的物体（在所有点上作用于当前速度垂直的向心力）。看看你会得到什么... - EOF

1

这不是浮点数问题。对于前向欧拉（您当前使用的积分方案），累积误差为O(Δt)，这与您的误差图相吻合（每种情况下它都减少了一个数量级，因为您在每种情况下将步长减小了一个数量级）。这是数值积分的基本问题。 - Kyle

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Kyle · Accepted Answer

这不是浮点数问题。实际上，即使您使用精确算法，您仍会看到相同的问题。

这个错误实际上是与数值积分本身以及您正在使用的特定方法和ODE有关的基本问题。

在这种情况下，您正在使用称为前向欧拉的积分方案。这可能是解决一阶ODE最简单的方法。当然，这也意味着它具有一些不良特征。

首先，它会在每个步骤引入误差。误差的大小为 O(Δt²)。这意味着在单个时间步长内的误差大致成比例于时间步长大小的平方。因此，如果将时间步长减半，则逐步误差大约降低到1/4。

但是，由于您缩小了时间步长，因此必须进行更多的步骤来模拟相同的时间。因此，您会增加更多但更小的误差。这就是累积误差为 O(Δt) 的原因。因此，如果采取的时间步长减半，那么在整个模拟时间内，您将获得一半的累积误差。

最终，这就是您所看到的累积误差。并且您可以在错误图中看到，每次将时间步数增加10倍时，最终误差会减少约10倍：因为时间步长缩小了10倍，所以总误差最终约小了10倍。

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

x' = -k * x

其中 k 是某个常数。这个 ODE 的精确解为 x(t) = x(0) * exp( -k * t )。因此，只要 k 为正，x 在时间增加时应该趋近于 0。

然而，如果我们尝试使用前向欧拉法来近似，我们得到的结果看起来像这样：

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

这是一个简单的递归关系，我们可以解决它：

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

现在，我们知道当t变大时，精确解趋近于0。但是，如果|1-k*Δt|<1，那么向前欧拉法解才会收敛于0。请注意这个式子取决于步长和 ODE 中的k项。如果k非常大，我们需要非常小的时间步长才能避免解析解爆炸。这就是它所拥有的条件稳定性的原因：解的稳定性取决于时间步长。

除此之外还存在其他问题，但这是一个广泛的话题，我不能在单个答案中涵盖所有内容。