c++ - 物理模拟给出(非常)不准确的简单轨迹微积分位置

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 重力场中简单下落)，我得到了不错的位置，但是对于更大(且相当随机)的数字，我得到的位置不准确。

那是 float 问题吗？

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

图表

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

Answer 1

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

这个错误对于数值积分本身以及您正在使用的特定方法和您正在求解的 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)

现在，我们知道十到 0 的精确解为 t变大。但是正向欧拉解只有在 |1 - k * Δt| < 1 时才会这样做。 .请注意该表达式如何取决于步长以及 k来 self 们的 ODE 的术语。如果k真的很大，我们需要一个非常非常小的时间步来防止解决方案爆炸。这就是它具有所谓的条件稳定性的原因:解的稳定性取决于时间步长。

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