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

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来说是非常重要的.

在这种情况下,您使用的是称为Forward Euler的集成方案.这可能是解决一阶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随着时间的推移而增加.

但是,如果我们尝试使用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)

现在,我们知道确切的解决方案从10到0 t变大.但是,前向欧拉解决方案只会这样做|1 - k * ?t| < 1.请注意该表达式如何取决于步长以及kODE中的术语.如果k真的非常大,我们需要一个非常小的时间步骤来防止解决方案爆炸.这就是为什么它拥有所谓的条件稳定性:解决方案的稳定性取决于时间步长.

还有其他一些问题,但这是一个广泛的主题,我无法在一个答案中涵盖所有内容.