用C++进行N体模拟

Question

我正在尝试实现一个二维的n体模拟的OpenMP版本。但是有一个问题：我假设每个粒子的初始速度和加速度都为零。当这些粒子第一次聚集在一起时，它们会以高速分散开来，而不会再次聚集在一起。这似乎与牛顿定律不一致，对吗？有人能解释一下为什么会发生这种情况并告诉我如何修复错误吗？以下是我的代码的一部分：

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {

    # pragma omp for schedule(static)
        for ( i = 0; i < g_N; ++i ) {
            g_pv[i].a_x = 0.0;
            g_pv[i].a_y = 0.0;
            for ( j = 0; j < g_N; ++j ) {
                if (i == j)
                    continue;
                double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
            } 
            g_pv[i].v_x += period * g_pv[i].a_x;
            g_pv[i].v_y += period * g_pv[i].a_y;
        }
    # pragma omp for schedule(static)
        for ( int i = 0; i < g_N; ++i ) {
            g_pv[i].pos_x += g_pv[i].v_x * period;
            g_pv[i].pos_y += g_pv[i].v_y * period;
        }
    }
}

不用担心OpenMP，把它当作一个顺序版本即可。OpenMP不会对结果产生太大影响。

编辑：为了澄清，这是整个代码（这部分可能有一些错误，但我描述的问题应该发生在上面的代码部分）

# include <iostream>
# include <fstream>
# include <iomanip>
# include <cmath>
# include <vector>
# include <cstdlib>
# include <omp.h>

# include <GL/glew.h>
# include <GL/freeglut.h>
# include <GL/gl.h>

using namespace std;

/* the size of the opengl window */
# define WIDTH 2000
# define HEIGHT 2000

/* define the global constants */
const double G = 6.67 * pow(10, -11);
// const double G = 6.67;
const double e = 0.00001;
const double period = 1;

/* define the structure of particle */
struct particle
{
    double m;
    double pos_x;
    double pos_y;
    double v_x;
    double v_y;
    double a_x;
    double a_y;

    particle(double m = 0, double pos_x = 0, double pos_y = 0, 
            double v_x = 0, double v_y = 0, double a_x = 0, double a_y = 0)
    {
        this->m         = m;
        this->pos_x     = pos_x;
        this->pos_y     = pos_y;
        this->v_x       = v_x;
        this->v_y       = v_y;
        this->a_x       = a_x;
        this->a_y       = a_y;
    }
};

/* define the global data */
int g_N;                        // number of particles
vector<particle> g_pv;          // particle vector

void setUp();
void update();
void display();

int main(int argc, char ** argv) {

    /* set up the window */
    glutInit(&argc, argv);
    glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
    glutInitWindowSize (WIDTH, HEIGHT);
    glutInitWindowPosition (100, 100);
    glutCreateWindow ("openmp");

    /* initialize */
    setUp();

    glutDisplayFunc(display);
    glutMainLoop();

    return 0;
}

/* read the input data */
void setUp() {
    glMatrixMode (GL_PROJECTION);
    glLoadIdentity ();
    /* Sets a 2d projection matrix
     * (0,0) is the lower left corner (WIDTH, HEIGHT) is the upper right */
    glOrtho (0, WIDTH, 0, HEIGHT, 0, 1);
    glDisable(GL_DEPTH_TEST);

    ifstream inFile;
    inFile.open("input_25.txt");

    inFile >> g_N;
    g_pv.resize(g_N);
    for ( int i = 0; i < g_N; ++i )
    {
        inFile >> g_pv[i].m >> g_pv[i].pos_x >> g_pv[i].pos_y
               >> g_pv[i].v_x >> g_pv[i].v_y >> g_pv[i].a_x >> g_pv[i].a_y; 
    }
    inFile.close();
}

/* display in openGL */
void display(void) {
    glClear(GL_COLOR_BUFFER_BIT);
    for(int i = 0; i < g_pv.size(); ++i) {
        /* Get the ith particle */
        particle p = g_pv[i];

        /* Draw the particle as a little square. */
        glBegin(GL_QUADS);
        glColor3f (1.0, 1.0, 1.0);
        glVertex2f(p.pos_x + 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y - 2);
        glVertex2f(p.pos_x + 2, p.pos_y - 2);
        glEnd();
    }   
    update();
    glutPostRedisplay();
    glFlush ();
}

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {
        /* compute the force */
        # pragma omp for schedule(static)
            for ( i = 0; i < g_N; ++i ) {
                g_pv[i].a_x = 0.0;
                g_pv[i].a_y = 0.0;
                for ( j = 0; j < g_N; ++j ) {
                    if (i == j)
                        continue;
                    double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                    g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                    g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
                } 
                g_pv[i].v_x += period * g_pv[i].a_x;
                g_pv[i].v_y += period * g_pv[i].a_y;
            }

        /* compute the velocity */
        # pragma omp for schedule(static)
            for ( int i = 0; i < g_N; ++i ) {
                g_pv[i].pos_x += g_pv[i].v_x * period;
                g_pv[i].pos_y += g_pv[i].v_y * period;
            }
    }
}

- user3558391

---

你的 j 循环似乎正在访问 i 循环尚未初始化的数组元素。也许在使用 i 和 j 的双重循环之前先运行整个 i 循环以进行初始化？ - Galik

简化代码可以降低出错的概率。 - Daniel Daranas

我怀疑这里的问题至少部分原因是您正在使用的积分方法不准确（一阶欧拉方法）。您考虑过改用Verlet积分法吗？当粒子彼此碰撞并合并时，位势会达到奇点，可能会导致大量数值误差发生。 - Rufflewind

我知道，那是我的观点。在担心OpenMP之前，请确保顺序版本能够正确运行。 - Z boson

我制作了以下输入文件：2 100000000 500 500 0 0 0 0 1 530 500 0 0.005 0 0它的行为与我预期的完全一致，其中一个行星绕着另一个行星运转。现在，如果将 0.005 替换为 0，较轻的行星将与较重的行星碰撞。也就是说，力最终会变得几乎无限，导致爆炸（参见例如 http://www.gromacs.org/Documentation/Terminology/Blowing_Up）。在短距离范围内添加排斥力（例如 Lennard-Jones 类型），情况可能会好些。 - alarge

Answer 1

我正在将我的评论扩展为答案（如Z boson所建议的），并提出了几个解决问题的建议。

这个问题确实属于计算科学.SE，因为我认为代码本身没有问题，而是算法似乎有问题：当粒子靠近时，你会得到一个力 G / pow(e,1.5) ~ G * 1e7。这很大，非常非常大（与您的时间步长相比）。为什么？假设你有两个行星，一个庞大的坐在（0,0），另一个小行星在（0,1），后者向前者获得非常大的加速度。在下一步中，小行星将位于（0，-100）或其他位置，并且其上的力为零，这意味着它永远不会返回，并且现在也具有相当大的速度。您的模拟已经崩溃。

正如你所说，这与牛顿定律不符，因此表明你的数值方案失败了。不要担心，有几种方法可以解决这个问题。你已经预料到了这一点，因为你添加了e。只需将其增大，比如设为10，你就应该没问题了。或者，设置一个非常小的时间步长period。你也可以在粒子过于接近时不计算力，或者有一些启发式方法来处理行星碰撞（也许它们会爆炸和消失，或者抛出异常）。或者使用r^-2势的斥力： g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) * (1 / pow(r_2 + e,1.5) - 1 / pow(r_2 + e,2.5)); 这类似于现象学Lennard-Jones interaction中包含的由泡利排斥原理引起的排斥作用。请注意，您可以增加排斥的锐度（将2.5改为12.5），这意味着远离时排斥作用较小（好），但需要更准确地解决，导致period更小（不好）。理想情况下，你应该从不会导致碰撞的初始配置开始，但这几乎不可能预测。最终，你可能需要使用上面列出的方法的组合。

使用OpenMP可能会导致轻微加速，但你应该真正使用长程力的算法，例如Barnes-Hut simulation。更多信息，请参见2013 Summer School on Fast Methods for Long Range Interactions in Complex Particle Systems以及他们发布的关于最新发展的手册（可免费获取）。你也可能不想显示每个时间步骤：科学模拟通常保存每5000个步骤左右。如果你想要漂亮的电影效果，可以使用一些移动平均插值来平滑噪声（你的模拟没有温度或类似的东西，所以你很可能不需要平均）。此外，我不确定你的数据结构是否为该工作进行了优化，或者你是否会遇到缓存未命中的问题。我并不是一个专家，所以也许其他人可以在这方面提供帮助。最后，考虑不使用pow，而是使用快速反平方根或类似的方法。

- alarge

---

我同意 OP 对所使用的方法实现没有问题，但是 OP 期望从该方法中得到什么，这是值得商榷的。不过我认为 OP 的问题属于物理标签的范畴（而 OpenMP 是一个主要的误导因素）。 - Z boson

我建议使用微小的切向速度是为了避免使用 e 调整因子。使用切向速度可以通过使轨迹变成二维而不是一维来避免奇点。但现在我意识到，微小的切向速度可能意味着轨迹是双曲线而不是轨道。 - Z boson

我不反对你回答中的任何内容，我只是想让你更加精确地表达“这与牛顿定律不相符”。可能有几种方法可以证明这一点。 - Z boson

@Zboson 问题在于，特别是当有多个行星时，你不能假设轨道的稳定性等等。是的，对于两个行星的情况，如果你确保有足够的切向速度，使得行星不会过于接近，那么你可以相对安全地放弃“调整因子”。现在，如果有25个行星（比如说），就没有任何保证了，你必须想办法解决这个问题（其中一种方法是使用 e）。 - alarge

@Zboson - 开始的代码使用了辛积分器。基本欧拉方法使用公式x(t+Δt)=x(t)+v(t)Δt, v(t+Δt)=x(t)+a(x(t))Δt。辛欧拉方法使用公式v(t+Δt)=x(t)+a(x(t))Δt, x(t+Δt)=x(t)+v(t+Δt)Δt。OP使用后者。问题不在于缺乏辛性，而是技术的低阶性质。 - David Hammen

Answer 2

The first problem is that you are using a simple "Taylor integration scheme". The point is that when approximating the infinitely small dt to a finite time difference Δt, which is your period, you are expanding the well-known motion equations in a Taylor expansion with infinite terms, which you are truncating.

In general: x_t+Δt = x_t + x′_t Δt + 1/2 x″_t Δt² + 1/6 x′″_t Δt³ + O(Δt⁴)

在这里，x′_t代表时间t的一阶导数，即速度v；x″_t代表二阶导数，即加速度a；O(Δt⁴)是误差的阶数 - 在这个例子中，我们将展开截断到第3阶，因此有一个本地误差为第4阶。

在您的情况下，您正在使用（欧拉方法）：

x_t+Δt = x_t + v_t Δt + O(Δt²)

和

v_t+Δt = v_t + a_t Δt + O(Δt²)

在这种情况下，由于您停止了对第一项的展开，因此这是一阶近似。您将获得一个Δt²的O级局部误差，对应于Δt的O级全局误差。这就是截断误差的行为方式 - 有关详细信息，请参见reference。

我非常了解两种不同的积分方案，可以改进欧拉方法：

- Verlet Integration（通常用于分子动力学代码） - Leapfrog Integration 我还知道龙格-库塔方法（从引用的两篇文章中找到参考资料），它们甚至具有更高的阶数，但我个人从未使用过它们。

我在这里提到的顺序是近似截断的顺序，严格相关于通过截断本身执行的误差的顺序。

Leapfrog方法是二阶的。

Verlet方法是三阶的，局部误差为O(Δt⁴)。即使您在推导中没有看到任何三阶导数，这也是一种三阶方法，因为它们在推导中被抵消了，如wikipedia reference所示。

更高阶的积分方案可以让您获得更精确的轨迹，而不需要缩小时间步长Δt。

然而，这些积分方法最重要的特性，简单的Taylor前向积分法无论截断阶数如何都会缺失，包括：

- 可逆性 - 它们是Symplectic的，即它们保持能量守恒

从参考文献中您可以找到很多学习材料，但是一本好的N体书籍可以让您以更有结构的方式学习。

关于这两种方法之间非常重要的区别，只需最后注意一点：Verlet方法不能自动给出速度，需要作为随后的步骤进行评估（直接评估）。另一方面，Leapfrog方法确实同时评估位置和速度，但是在不同的时间 - 使用此方法无法同时评估这两个量。

一旦您拥有了良好的集成方案，您需要考虑能够容忍的最大误差是多少，此时需要进行一些误差分析，然后您将需要实施多时间尺度方法，即使在两个物体过于接近甚至对于O(Δt⁴)或更大的情况下（考虑到重力摆），也能进行准确的集成。

它们可以是全局的（即在整个系统中减少“周期”）或局部的（仅在某些粒子过于接近的系统分区中执行）...但在这一点上，您应该能够自行查找更多参考资料。

- Sigi

---

我很高兴你提到了OP的方法（欧拉方法）不能保持能量不变。 - Z boson

值得一提的是，Runge-Kutta方法不是辛波立克结构的，因此通常不适用于这种类型的模拟。例如，http://arxiv.org/abs/1412.5187的第18页有关于每种方法表现的精美插图（大多数介绍性书籍也有）。 - alarge

有辛普勒克 Runge-Kutta 方法（例如辛普勒克分块 Runge-Kutta）。此外，是否重视辛普勒性取决于感兴趣的时间跨度。如果时间跨度小于约一百万年（并且系统没有碰撞），高阶、高精度技术可能比任何辛普勒技术都更好。如果系统存在碰撞，仅靠辛普勒技术是不起作用的，需要进行大量的细心调试。将初始速度设为零可以确保系统没有碰撞。 - David Hammen

我给了@alarge奖励。我认为你和他的答案都同样好（投票者似乎也这么认为）。他回答得最早。如果我能分割奖励，我会这么做，但我不能。OP的问题的一部分是“这似乎与牛顿定律不一致，对吧？”我希望看到一个数值例子，说明他选择的方法为什么违反了牛顿物理学（而不仅仅是说它违反了），但至少你提到了能量守恒被破坏。 - Z boson

我突然意识到我之前的评论可能不是很清楚。我的意思是，我完全同意你的观点，即应该使用Verlet / velocity Verlet / leapfrog（或可能是更奇特的方法），这确实是我在回答中应该提到的。然而，这些方法并没有直接解决爆炸的根本问题，仅仅说这些积分器应该被使用并不能解释为什么会发生这种情况。无论如何，我对你所说的一切都不持反对态度。 - alarge