Алгоритм верлета скорости - увеличение энергии для задачи n тел

Я выяснил, в чем дело.

Проблема оказалась в обновлении положения тел по одному. Вычисление ускорения предполагает, что никакое тело не было перемещено между временными шагами, однако обновление одного за другим привело к тому, что некоторые тела имели положение от t, а некоторые от t + dt. Это различие в этой конкретной системе привело к тому, что векторная разность скорости движущегося по орбите тела не была идеально направлена ​​к центру масс. Фактически генерировалась небольшая тангенциальная составляющая, и в систему добавлялась энергия. Ошибка была небольшая, но со временем она накапливалась и становилась заметной.

Я исправил проблему, выполнив каждый этап алгоритма верлета сразу для всех тел. Вот исправленный код интегратора:

    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';
   }
}