代码之家 › 专栏 › 技术社区 › Ronny

如何防止程序返回“inf”?

numeric floating-point c++

Ronny · 技术社区 · 2 年前

我正在研究一个著名的经典物理问题:三体问题。我制作了以下程序:-

#include <math.h>
#include <iostream>
using namespace std;

class System_variables
{
    private:

        // These are the system variables:- position, velocity and accleration
        long double body_1_pos[3];
        long double body_2_pos[3];
        long double body_3_pos[3];

        long double body_1_v[3];
        long double body_2_v[3];
        long double body_3_v[3];

        long double acc_body_1[3];
        long double acc_body_2[3];
        long double acc_body_3[3];

        // Needed constants
        const long double body_mass = 6 * pow(10, 24);
        double time_step;
        const long double G = 6.67 * pow(10, -11);

    public:


        System_variables(long double body_1_pos_i[3], long double body_2_pos_i[3], long double body_3_pos_i[3], long double body_1_v_i[3], long double body_2_v_i[3], long double body_3_v_i[3])
        {
            // Initiates the system variables, position and velocity.
            for (int m = 0; m<=2; m++)
            {
                body_1_pos[m] = body_1_pos_i[m];
                body_2_pos[m] = body_2_pos_i[m];
                body_3_pos[m] = body_3_pos_i[m];

                body_1_v[m] = body_1_v_i[m];
                body_2_v[m] = body_2_v_i[m];
                body_3_v[m] = body_3_v_i[m];
            }
        }
        void init_time_step(double step)
        {
            // Initializes time step.
            time_step = step;
        }
        void calculate_acc()
        {   
            // Calculates the accleration on the bodies using newtons equations.
            for (int m=0; m<=2; m++)
            {
                acc_body_1[m] = G * body_mass/pow((body_2_pos[m] - body_1_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_1_pos[m]), 2);
                acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2);
                acc_body_3[m] = G * body_mass/pow((body_1_pos[m] - body_3_pos[m]), 2) + G * body_mass/pow((body_2_pos[m] - body_3_pos[m]), 2);
            };
            return;
        }

        void apply_acc_to_v()
        {   
            // Applys the accleration to the velocities of the bodies.
            for (int m=0; m<=2; m++)
            {
                body_1_v[m] += acc_body_1[m] * time_step;
                body_2_v[m] += acc_body_2[m] * time_step;
                body_3_v[m] += acc_body_3[m] * time_step;
            }
            return;
        }

        void apply_v_to_pos()
        {   
            // Applys the velocities to the positions of the bodies.
            for (int m=0; m<=2; m++)
            {
                body_1_pos[m] += body_1_v[m] * time_step;
                body_2_pos[m] += body_2_v[m] * time_step;
                body_3_pos[m] += body_3_v[m] * time_step;
            }
            return;
        }

        long double get_body_1_pos(int index)
        {
            // Returns the position of body 1 based on index.
            return body_1_pos[index];
        }

        long double get_body_2_pos(int index)
        {   
            // Returns the position of body 2 based on index.
            return body_2_pos[index];
        }

        long double get_body_3_pos(int index)
        {
            // Returns the position of body 3 based on index.
            return body_3_pos[index];
        }

        long double get_body_1_v(int index)
        {
            // Returns the velocity of body 1 based on index.
            return body_1_v[index];
        }

        long double get_body_2_v(int index)
        {
            // Returns the velocity of body 2 based on index.
            return body_2_v[index];
        }

        long double get_body_3_v(int index)
        {
            // Returns the velocity of body 3 based on index.
            return body_3_v[index];
        }


};

int main()
{
    long double array_1[3] = {100,100,100}; 
    long double array_2[3] ={10000000000, 10000000000, 10000000000};
    long double array_3[3] = {10000000000, 10000000000, -10000000000};
    long double array_0[3] =  {0,0,0};

    System_variables system(array_1, array_2, array_3, array_0, array_0, array_0);
    system.init_time_step(100);

    for (int time = 0; time <= 10; time++)
    {
        cout << "x position of body 1: " << system.get_body_1_pos(0) << endl 
        << "y position of body 1: " << system.get_body_1_pos(1) << endl 
        << "z position of body 1: " << system.get_body_1_pos(2) << endl << endl;

        cout << "x position of body 2: " << system.get_body_2_pos(0) << endl 
        << "y position of body 2: " << system.get_body_2_pos(1) << endl 
        << "z position of body 2: " << system.get_body_2_pos(2) << endl << endl;

        cout << "x position of body 3: " << system.get_body_3_pos(0) << endl 
        << "y position of body 3: " << system.get_body_3_pos(1) << endl 
        << "z position of body 3: " << system.get_body_3_pos(2) << endl << endl;

        system.calculate_acc();
        system.apply_acc_to_v();
        system.apply_v_to_pos();
    }
}

需要注意的一点是,与其使用经典 double 数据类型,我使用过 long double ,因为当我使用 双重的 数据类型,过一段时间 inf 正在输出,所以我使用 长双 认为它可以大精度地存储小数。但是,它不起作用。我也试过了 float (尽管精度较低),但输出始终相同:-

x position of body 1: 100
y position of body 1: 100
z position of body 1: 100

x position of body 2: 1e+10
y position of body 2: 1e+10
z position of body 2: 1e+10

x position of body 3: 1e+10
y position of body 3: 1e+10
z position of body 3: -1e+10

x position of body 1: 100.08
y position of body 1: 100.08
z position of body 1: 100.08

x position of body 2: inf
y position of body 2: inf
z position of body 2: 1e+10

x position of body 3: inf
y position of body 3: inf
z position of body 3: -1e+10

x position of body 1: 100.16
y position of body 1: 100.16
z position of body 1: 100.24

x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: 1e+10

x position of body 3: -nan(ind)
y position of body 3: -nan(ind)
z position of body 3: -1e+10

x position of body 1: -nan(ind)
y position of body 1: -nan(ind)
z position of body 1: 100.48

x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: 1e+10

x position of body 3: -nan(ind)
y position of body 3: -nan(ind)
z position of body 3: -1e+10

我正在考虑使用指针,但我怀疑这是否是一个正确的想法。你怎么认为?

1 回复 | 直到 2 年前

ComicSansMS 2 年前

你的计算在数值上不稳定。

在你的程序中,事情首先开始破裂的确切位置是加速度计算。例如:

acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) + 
                G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2);

Body_3和Body_2包含大小相似的非常大的数字。它们的差异将是一个非常小的数字。平方这个小数字会使它变得更小。然后你把质量除以这个非常小的数字,得到一个非常非常大的数字。这不好。

在您的情况下,body_3和body_2的数字是偶数 完全相同的 ,所以你在这里真的被0除了,这会直接把你送到 inf 。但是,即使数字只是有点相等,你也可能会在这里遇到问题。

请注意,这主要不是所用浮点数据类型精度有限的问题(尽管可以通过增加更多精度在一定程度上缓解这个问题),而是计算结构的问题。有一个完整的计算机科学学科叫做数值计算,专门处理这样的计算可能出现的微妙问题。

我建议你从对数值积分进行一些研究开始,这正是你在这里试图做的,并看看一些常见的陷阱,比如 catastrophic cancellation 问题

Vincent Saulue-Laborde 2 年前

你的 calculate_acc 不仅在数值上不稳定:而且完全不正确。

// X,Y,Z coordinates of the acceleration on "body1" caused by "body2"'s gravity according to OP
accelX = G * body_mass/pow((body_2_pos[0] - body_1_pos[0]), 2);
accelY = G * body_mass/pow((body_2_pos[1] - body_1_pos[1]), 2);
accelZ = G * body_mass/pow((body_2_pos[2] - body_1_pos[2]), 2);

上的加速度矢量 body1 由引起 body2 根据牛顿引力定律,的引力为:

// PSEUDOCODE
// p1 -> body_1_pos: vector position of body1
// p2 -> body_2_pos: vector position of body2
// m2: mass of body2
// norm(someVector): euclidian norm of someVector

accelerationNorm = G * m2 / (norm(p2 - p1))Â²; // 1d scalar, dimension: acceleration
unitVector = (p2 - p1) / norm(p2 - p1); // 3d vector, dimension: one
forceVector = accelerationNorm * unitVector; // 3d vector, dimension: acceleration

请注意 norm(p2 - p1) ,这在某种程度上从未在您的程序中计算过。你不能只除以 (p2[0]-p1[0])Â² 计算X坐标,然后 (p2[1]-p1[1])Â² 以计算Y坐标。您必须首先计算一个适当的范数: sqrt(xÂ² + yÂ² + zÂ²) ,然后使用该值来计算加速度向量的每个坐标。