c++ - Velocity verlet 算法 - n 体问题的能量增加

Question

问题

我实现了 Verlet 速度算法来计算 2 个物体在重力作用下相互作用的轨迹(仅限牛顿引力)。绕轨道运行的小天体质量很小，位于轨道中心的天体质量很大。

理论上，Velocity Verlet 不应改变系统的总能量(它会振荡，但随着时间的推移，平均值将保持接近初始能量)。

但在实践中，我观察到能量会随着时间的推移而增加。

结果

下面是一些说明问题的结果。所有模拟均以时间步长 dt = 0.001 执行。轨道天体质量为1000，宇宙万有引力常数G=1.0

在所有情况下，较小的物体初始位置为 {0, 0, 1}，初始速度为 {0, 32, 0}。较大物体的初始速度为 {0,0,0}。

案例 1(小体重 = 0.00001)

这是较小的 body 的轨迹:

这是超过 10 万步的能量。

正如我们所见，能量变化不大。微小的变化可能是由于计算不准确造成的。

案例 1(小体重 = 0.001)

这是轨道天体的轨迹:

这是总能量:

正如我们现在看到的，系统正在获得能量。

案例 3(小体重 = 1)

这是轨道天体的轨迹:

这是总能量:

现在我们获得了很多能量。

代码

这是执行所有计算的源代码:

高级积分器代码:

void Universe::simulation_step()
{
    for(std::size_t i=0; i<get_size(); i++)
    {
        // Verlet step 1: Compute v(t + dt/2) = v(t) + 0.5*dt*a(t)
        const Vector3D<Real> vel_half_step = {
            velocity(i, 0) + static_cast<Real>(0.5)*sim_config.timestep*acceleration(i, 0),
            velocity(i, 1) + static_cast<Real>(0.5)*sim_config.timestep*acceleration(i, 1),
            velocity(i, 2) + static_cast<Real>(0.5)*sim_config.timestep*acceleration(i, 2)
        };

        // Verlet step 2: Compute x(t + dt) = x(t) + v(t + dt/2)*dt
        position(i, 0) += vel_half_step.x*sim_config.timestep;
        position(i, 1) += vel_half_step.y*sim_config.timestep;
        position(i, 2) += vel_half_step.z*sim_config.timestep;

        // Verlet step 3: update forces and update acceleration.
        const Vector3D<Real> forces = compute_net_grawitational_force(i);
        acceleration(i, 0) = forces.x/mass(i);
        acceleration(i, 1) = forces.y/mass(i);
        acceleration(i, 2) = forces.z/mass(i);

        // Verlet step 4: update velocity to the full timestep.
        velocity(i, 0) = vel_half_step.x + static_cast<Real>(0.5)*sim_config.timestep*acceleration(i, 0);
        velocity(i, 1) = vel_half_step.y + static_cast<Real>(0.5)*sim_config.timestep*acceleration(i, 1);
        velocity(i, 2) = vel_half_step.z + static_cast<Real>(0.5)*sim_config.timestep*acceleration(i, 2);
    }

    sim_time += sim_config.timestep;
}

下面是计算作用在 body 上的净重力的代码:

Vector3D<Real> Universe::compute_net_grawitational_force(std::size_t i)
{
    Vector3D<Real> accumulated_force = {0,0,0};
    const Vector3D<Real> r2 = {
            position(i, 0),
            position(i, 1),
            position(i, 2)
    };

    const Real m1 = mass(i);

    for(std::size_t k=0; k<get_size(); k++)
    {
        if(k == i)
            continue;
        
        const Vector3D<Real> distace_vec = {
            r2.x - position(k, 0),
            r2.y - position(k, 1),
            r2.z - position(k, 2),
        };
        
        const Real distance = distace_vec.norm2();

        // Compute term that will be multipled by distance vector.
        const Real a = (-1*sim_config.G*m1*mass(k))/
        (distance*distance*distance);

        // Compute and add new force acting on the body.
        accumulated_force.x += distace_vec.x*a;
        accumulated_force.y += distace_vec.y*a;
        accumulated_force.z += distace_vec.z*a;
    }

    return accumulated_force;
}

这是实现 norm2() 的代码:

template<typename T>
struct Vector3D
{
    T x;
    T y;
    T z;
    
    T norm2() const
    {
        return sqrt(x*x + y*y + z*z);
    }
};

最后是计算先前绘制的结果的代码:

    std::vector<Real> x, y, z, energy;
    x.resize(NSTEPS);
    y.resize(NSTEPS);
    z.resize(NSTEPS);
    energy.resize(NSTEPS);

    for(std::size_t i=0; i<NSTEPS; i++)
    {
        universe.simulation_step();
        const Vector3D<Real> pos1 = 
        {
            universe.get_positions()(0, 0),
            universe.get_positions()(0, 1),
            universe.get_positions()(0, 2)
        };

        const Vector3D<Real> pos2 = 
        {
            universe.get_positions()(1, 0),
            universe.get_positions()(1, 1),
            universe.get_positions()(1, 2)
        };

        x[i] = pos2.x;
        y[i] = pos2.y;
        z[i] = pos2.z;

        // Compute total kinetic energy of the system.
        const Vector3D<Real> vel1 = 
        {
            universe.get_velocities()(0, 0),
            universe.get_velocities()(0, 1),
            universe.get_velocities()(0, 2),
        };
        const Vector3D<Real> vel2 = 
        {
            universe.get_velocities()(1, 0),
            universe.get_velocities()(1, 1),
            universe.get_velocities()(1, 2),
        };

        const Real mass1 = universe.get_masses()(0);
        const Real mass2 = universe.get_masses()(1);
        const Real spd1 = vel1.norm2();
        const Real spd2 = vel2.norm2();

        energy[i] = (spd1*spd1)*mass1*static_cast<float>(0.5);
        energy[i] += (spd2*spd2)*mass2*static_cast<float>(0.5);

        // Compute total potential energy
        const Vector3D<Real> distance_vec = 
        {
            pos1.x - pos2.x,
            pos1.y - pos2.y,
            pos1.z - pos2.z
        };

        const Real G = universe.get_sim_config().G;

        energy[i] += -G*((mass1*mass2)/distance_vec.norm2());
    }

类型 Real 是 float。

我的理论

我是数值积分的初学者(这就是我在这里发布这个问题的原因)。然而，这里有一些关于可能出错的理论:

• 当涉及到 n>=2 时，Velocity Verlet 算法存在一些缺陷，我已经陷入其中。
• 上述代码中某处存在实现错误，但我没有看到。
• float 计算的误差由于小的移动而累积庞大的身躯。 (可能不是这种情况，请参阅下面的编辑。)
• 在调试过程中，我遇到了分子动力学模拟中的能量漂移。也许这就是这里发生的事情？

轨道似乎没有解体，但不是我预期的结果我想知道为什么。

谁能帮我解开这个谜团？

编辑:

我测试了 double ，唯一的变化是现在最小轨道质量的能量更加稳定。

现在即使是最小的质量也可以看到增加的趋势。这暗示这不是计算精度的问题。

Answer 1

我发现了问题。

结果是一个问题是一个一个地更新 body 的位置。加速度的计算假设没有物体在时间步之间移动然而，一个一个地更新导致一些物体的位置来自 t 而一些物体来自 t + dt。这个特定系统中的差异导致轨道体速度的 vector 差异无法理想地指向质心。实际上，产生了一个小的切向分量，并且能量被添加到系统中。错误很小，但随着时间的推移它会累积并变得可见。

我通过一次对所有物体执行 verlet 算法的每个阶段来解决这个问题。这是集成器的修订代码:

    for(std::size_t i=0; i<get_size(); i++)
    {
        position(i, 0) += velocity(i, 0)*sim_config.timestep + acceleration(i, 0)*sim_config.timestep*sim_config.timestep*static_cast<Real>(0.5);
        position(i, 1) += velocity(i, 1)*sim_config.timestep + acceleration(i, 1)*sim_config.timestep*sim_config.timestep*static_cast<Real>(0.5);
        position(i, 2) += velocity(i, 2)*sim_config.timestep + acceleration(i, 2)*sim_config.timestep*sim_config.timestep*static_cast<Real>(0.5);
    }
    
    for(std::size_t i=0; i<get_size(); i++)
    {
        velocity(i, 0) += acceleration(i, 0)*sim_config.timestep*static_cast<Real>(0.5);
        velocity(i, 1) += acceleration(i, 1)*sim_config.timestep*static_cast<Real>(0.5);
        velocity(i, 2) += acceleration(i, 2)*sim_config.timestep*static_cast<Real>(0.5);
    }

    for(std::size_t i=0; i<get_size(); i++)
    {
        const Vector3D<Real> forces = compute_net_grawitational(i);
        acceleration(i, 0) = forces.x/mass(i);
        acceleration(i, 1) = forces.y/mass(i);
        acceleration(i, 2) = forces.z/mass(i);
    }

    for(std::size_t i=0; i<get_size(); i++)
    {
        velocity(i, 0) += acceleration(i, 0)*sim_config.timestep*static_cast<Real>(0.5);
        velocity(i, 1) += acceleration(i, 1)*sim_config.timestep*static_cast<Real>(0.5);
        velocity(i, 2) += acceleration(i, 2)*sim_config.timestep*static_cast<Real>(0.5);

现在即使对于最重的轨道天体，能量也不会增加:

能量仍在漂移，但随着时间的推移，差异似乎趋于平衡，并且相对于总值而言变化很小。绘图是自动调整范围的，所以变化看起来很大，但它们在总能量的 +- 1% 以内，这对我的应用来说是可以接受的。