物理模拟给出（非常）不准确的简单轨迹微积分位置

Question

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

为了正确地做到这一点，我首先想知道这种获取位置的方法与运动学方程之间的差异有多大。所以我编写了这段代码，用于存储模拟和文件中的方程式给出的位置 (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）。

图表

黑色+虚线：来自运动方程的值
红色：每秒 100 点
橙色：每秒 1000 点
绿色：每秒 10000 点

Answer 1

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

这个错误对于数值积分本身以及您正在使用的特定方法和您正在求解的 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 时才这样做。请注意该表达式如何取决于步长以及 ODE 中的 k 项。如果 k 真的很大，我们需要一个非常非常小的时间步来防止解决方案爆炸。这就是它具有所谓的条件稳定性的原因：解的稳定性取决于时间步长。

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

物理模拟给出（非常）不准确的简单轨迹微积分位置

Physics simulation gives (very) inaccurate positions for simple trajectories calculus

c++

simulation

physics