问题
我实现了速度 verlet 算法来计算 2 个物体在重力作用下相互作用的轨迹（仅限牛顿重力）。轨道较小的物体质量很小，位于轨道中心的物体质量很大。
理论上，Velocity Verlet 不应该改变系统的总能量（它会振荡，但随着时间的推移，平均值将保持接近初始能量）。
然而在实践中，我观察到能量随着时间的推移而增加。
结果
以下是一些说明问题的结果。所有模拟均以时间步长 dt=0.001 进行。轨道物体的质量为 1000，宇宙的万有引力常数设定为 G=1.0
在所有情况下，较小的身体初始位置是 {0, 0, 1}，它的初始速度是 {0, 32, 0}。较大物体的初始速度为{0,0,0}。
案例 1（小体重 = 0.00001）
这是较小的身体的轨迹：
这是超过 100k 步的能量。
正如我们所见，能量变化不大。小的变化可能是由于计算不准确造成的。
案例 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;
}
这是计算作用在物体上的净重力的代码：
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 算法存在一些陷阱，我已经陷入其中。
- 上面的代码中某处存在实现错误，我没有看到。
- 浮点数计算导致的误差由于大体的小运动而累积。（可能不是这种情况，请参阅下面的编辑。）
- 在尝试调试这个的过程中，我遇到
Energy drift了分子动力学模拟。也许这就是这里发生的事情？
轨道似乎没有分崩离析，但这不是我预期的结果，我想知道为什么。
有人可以帮我解开这个谜吗？
编辑：
我测试了双精度，唯一的变化是现在最小轨道质量的能量更加稳定。
现在即使对于最小的质量也可以看到增加的趋势。这暗示这不是计算精度的问题。
