问题描述
问题
我实现了veLocity verlet 算法来计算两个物体在重力作用下相互作用的轨迹（仅限牛顿重力）。绕轨道运行的小天体质量非常小，位于轨道中心的天体质量大。
理论上 VeLocity Verlet 不应改变系统的总能量（它会振荡，但随着时间的推移，平均值将保持接近初始能量）。
结果
以下是一些说明问题的结果。 所有模拟均以时间步长 dt=0.001 进行。 轨道天体质量为1000，宇宙万有引力常数G=1.0
在所有情况下，较小的身体初始位置为 {0,1}，其初始速度为 {0,32,0}。 较大物体的初始速度为 {0,0}。
案例 1（小体重 = 0.00001）
这是较小身体的轨迹：
这是超过 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),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};
const Vector3D<Real> r2 = {
position(i,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,r2.y - position(k,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;
}
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,universe.get_positions()(0,2)
};
const Vector3D<Real> pos2 =
{
universe.get_positions()(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,universe.get_veLocities()(0,};
const Vector3D<Real> vel2 =
{
universe.get_veLocities()(1,universe.get_veLocities()(1,};
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。也许这就是这里发生的事情？
轨道似乎没有分崩离析，但这不是我预期的结果 我想知道为什么。
有人能帮我解开这个谜吗？
编辑：
我已经测试过双精度，唯一的变化是现在最小轨道质量的能量更加稳定。
现在即使是最小的质量也可以看到增加的趋势。 这说明计算精度没有问题。
解决方法
我发现出了什么问题。
结果是一个问题是一个一个地更新身体的位置。加速度的计算假设没有物体在时间步之间移动 然而，一一更新导致一些主体的位置来自 t，一些主体来自 t + dt。 这个特定系统中的差异导致轨道体速度的矢量差异不理想地指向质心。 实际上，产生了一个小的切向分量，并向系统添加了能量。错误很小，但随着时间的推移它会累积并可见。
我通过一次在所有物体上执行 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,2)*sim_config.timestep*static_cast<Real>(0.5);
现在即使是最重的轨道天体，能量也没有增加：
能量仍在漂移，但随着时间的推移，差异似乎趋于平均，并且相对于总值的变化很小。绘图是自动调整范围的，因此变化看起来很大，但它们在总能量的 +- 1% 以内，这是我的应用可以接受的。
,在 OP 找到答案之前，出于好奇，我正在努力重现该问题；我正在开始调试阶段（快速浏览第一个数据会发现至少存在一个主要错误），但我想我会放弃这个任务（假期快结束了）。在从我的硬盘中删除代码之前，我想无论如何都要把它放在这里，为了后代。
#include <cmath>
#include <limits>
#include <iostream>
#include <fstream>
#include <string>
#include <exception>
#include <vector>
//----------------------------------------------------------------------
// Some floating point facilities
//----------------------------------------------------------------------
// Some floating point facilities
constexpr inline bool is_zero(double x) noexcept
{
return std::fabs(x) < std::numeric_limits<double>::epsilon();
}
constexpr inline double ratio(const double n,const double d) noexcept
{
return is_zero(d) ? n/std::numeric_limits<double>::epsilon() : n/d;
}
////////////////////////////////////////////////////////////////////////
struct Vector3D ////////////////////////////////////////////////////////
{
double x,y,z;
Vector3D() noexcept : x(0.0),y(0.0),z(0.0) {}
Vector3D(const double X,const double Y,const double Z) noexcept
: x(X),y(Y),z(Z) {}
Vector3D operator+=(const Vector3D& other) noexcept
{
x+=other.x; y+=other.y; z+=other.z;
return *this;
}
Vector3D operator-() const noexcept
{
return Vector3D{-x,-y,-z};
}
friend Vector3D operator+(const Vector3D& v1,const Vector3D& v2) noexcept
{
return Vector3D{v1.x+v2.x,v1.y+v2.y,v1.z+v2.z};
}
friend Vector3D operator-(const Vector3D& v1,const Vector3D& v2) noexcept
{
return Vector3D{v1.x-v2.x,v1.y-v2.y,v1.z-v2.z};
}
friend Vector3D operator*(double k,const Vector3D& v) noexcept
{
return Vector3D{k*v.x,k*v.y,k*v.z};
}
friend Vector3D operator/(const Vector3D& v,double k) noexcept
{
return Vector3D{v.x/k,v.y/k,v.z/k};
}
friend std::ostream& operator<<(std::ostream& os,const Vector3D& v)
{
os << v.x << ',' << v.y << ',' << v.z;
return os;
}
double norm2() const noexcept { return x*x + y*y + z*z; }
double norm() const noexcept { return std::sqrt(norm2()); }
};
////////////////////////////////////////////////////////////////////////
class Body /////////////////////////////////////////////////////////////
{
public:
Body(const double m,const Vector3D& pos,const Vector3D& spd) noexcept
: i_mass(m),i_pos(pos),i_spd(spd) {}
double mass() const noexcept { return i_mass; }
const Vector3D& position() const noexcept { return i_pos; }
const Vector3D& speed() const noexcept { return i_spd; }
const Vector3D& acceleration() const noexcept { return i_acc; }
Vector3D distance_from(const Body& other) const noexcept
{
return position() - other.position();
}
double kinetic_energy() const noexcept
{// ½ m·V²
return 0.5 * i_mass * i_spd.norm2();
}
Vector3D gravitational_force_on(const Body& other,const double G) const noexcept
{// G · M·m / d²
Vector3D disp = distance_from(other);
double d = disp.norm();
return ratio(G * mass() * other.mass(),d*d*d) * disp;
}
double gravitational_energy_with(const Body& other,const double G) const noexcept
{// U = -G · mi·mj / d
double d = distance_from(other).norm();
return ratio(G * mass() * other.mass(),d);
}
void apply_force(const Vector3D& f)
{// Newton's law: F=ma
i_acc = f / i_mass;
}
void evolve_speed(const double dt) noexcept
{
i_spd += dt * i_acc;
}
void evolve_position(const double dt) noexcept
{
i_pos += dt * i_spd;
}
private:
double i_mass;
Vector3D i_pos,// Position [<space>]
i_spd,// Speed [<space>/<time>]
i_acc; // Acceleration [<space>/<time>²]
};
////////////////////////////////////////////////////////////////////////
class Universe /////////////////////////////////////////////////////////
{
public:
Universe(const double g) noexcept : G(g) {}
void evolve(const double dt) noexcept
{
for(Body& body : i_bodies)
{
body.evolve_speed(dt/2);
body.evolve_position(dt);
}
for( auto ibody=i_bodies.begin(); ibody!=i_bodies.end(); ++ibody )
{
Vector3D f = gravitational_force_on_body(ibody);
ibody->apply_force(f);
ibody->evolve_speed(dt/2);
}
}
double kinetic_energy() const noexcept
{// K = ∑ ½ m·V²
double Ek = 0.0;
for( const Body& body : i_bodies ) Ek += body.kinetic_energy();
return Ek;
}
double gravitational_energy() const noexcept
{// U = ½ ∑ -G · mi·mj / d
double Eu = 0.0;
for( auto ibody=i_bodies.begin(); ibody!=i_bodies.end(); ++ibody )
{// Iterate over all the other bodies
for( auto ibody2=i_bodies.begin(); ibody2!=ibody; ++ibody2 )
Eu += ibody->gravitational_energy_with(*ibody2,G);
for( auto ibody2=ibody+1; ibody2!=i_bodies.end(); ++ibody2 )
Eu += ibody->gravitational_energy_with(*ibody2,G);
}
return Eu/2;
}
double total_energy() const noexcept { return kinetic_energy() + gravitational_energy(); }
Vector3D center_of_mass() const noexcept
{// U = ∑ m·vpos / M
Vector3D c;
double total_mass = 0.0;
for( const Body& body : i_bodies )
{
c += body.mass() * body.position();
total_mass += body.mass();
}
return c/total_mass;
}
Vector3D gravitational_force_on_body( std::vector<Body>::const_iterator ibody ) const noexcept
{// F = ∑ G · m·mi / di²
Vector3D f;
// Iterate over all the other bodies
for( auto ibody2=i_bodies.begin(); ibody2!=ibody; ++ibody2 )
f += ibody2->gravitational_force_on(*ibody,G);
for( auto ibody2=ibody+1; ibody2!=i_bodies.end(); ++ibody2 )
f += ibody2->gravitational_force_on(*ibody,G);
return f;
}
void add_body(const double m,const Vector3D& spd)
{
i_bodies.emplace_back(m,pos,spd);
}
const std::vector<Body>& bodies() const noexcept { return i_bodies; }
const double G; // Gravitational constant
private:
std::vector<Body> i_bodies;
};
////////////////////////////////////////////////////////////////////////
class Simulation ///////////////////////////////////////////////////////
{
public:
class Data /////////////////////////////////////////////////////////
{
public:
struct Sample ///////////////////////////////////////////////////
{
double time;
std::vector<Body> bodies; // Bodies status
Vector3D Cm; // Center of mass
double Ek,// Kinetic energy
Eu; // Potential energy
Sample(const double t,const std::vector<Body>& bd,const Vector3D& cm,const double ek,const double eu) : time(t),bodies(bd),Cm(cm),Ek(ek),Eu(eu) {}
};
void init(const std::vector<std::string>::size_type N) noexcept
{
i_samples.clear();
i_samples.reserve(N);
}
void add(Sample&& s)
{
i_samples.push_back(s);
}
void save_to_file(std::string fpath) const
{
std::ofstream f (fpath,std::ios::out);
if(!f.is_open()) throw std::runtime_error("Unable to open file " + fpath);
//f << "time,total-energy\n";
//for(const Sample& s : i_samples)
// f << s.time << ',' << (s.Ek+s.Eu) << '\n';
f << "time,bodies-xyz-pos\n";
for(const Sample& s : i_samples)
{
f << s.time;
for(const Body& body : s.bodies)
f << ',' << body.position();
f << '\n';
}
}
const std::vector<Sample>& samples() const noexcept { return i_samples; }
private:
std::vector<Sample> i_samples;
};
// Total time Time increment
void execute(Universe& universe,const double T,const double dt)
{
auto N = static_cast<std::size_t>(T/dt + 1);
i_data.init(N);
double t = 0.0;
do {
universe.evolve(dt);
i_data.add( Data::Sample(t,universe.bodies(),universe.center_of_mass(),universe.kinetic_energy(),universe.gravitational_energy() ) );
t += dt;
}
while(t<T);
}
const Data& data() const noexcept { return i_data; }
private:
Data i_data;
};
//----------------------------------------------------------------------
int main()
{
// Creating a universe with a certain gravitational constant
Universe universe(1); // Our universe: 6.67408E-11 m³/kg s²
// Adding bodies (mass,initial position and speed)
universe.add_body(1000,Vector3D(0,0),0));
universe.add_body(100,Vector3D(10,10,0));
// Simulation settings
Simulation sim;
const double T = 100; // Total time
const double dt = 0.001; // Iteration time
std::cout << "Simulating T=" << T << " dt=" << dt << "...";
sim.execute(universe,T,dt);
std::cout << "...Done";
// Now do something with the simulation data...
// ...Edit as desired
//sim.data().save_to_file("data.txt");
{// Energies
std::string fname = "energies.txt";
std::ofstream f (fname,std::ios::out);
if(!f.is_open()) throw std::runtime_error("Unable to open file " + fname);
f << "time,kinetic,potential,total\n";
for(const Simulation::Data::Sample& s : sim.data().samples())
f << s.time << ',' << s.Ek << ',' << s.Eu << ',' << (s.Ek+s.Eu) << '\n';
}
{// Positions of...
std::size_t idx = 1; // ...Second body
std::string fname = "body" + std::to_string(idx) + ".txt";
std::ofstream f (fname,body" << idx << ".x,body" << idx << ".y,body" << idx << ".z\n";
for(const Simulation::Data::Sample& s : sim.data().samples())
f << s.time << ',' << (s.bodies.begin()+idx)->position() << '\n';
}
}