速度 verlet 算法 - n 体问题的能量增加

我发现了问题所在。

结果发现一个问题是逐个更新尸体的位置。加速度的计算假设在时间步长之间没有物体移动 然而，一个接一个地更新导致一些物体的位置来自 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, 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% 以内，这对于我的应用来说是可以接受的。

在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& pos, 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), Vector3D(0,0,0));
    universe.add_body(100, Vector3D(10,0,0), Vector3D(0,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, std::ios::out);
    if(!f.is_open()) throw std::runtime_error("Unable to open file " + fname);
    f << "time,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';
   }
}