C++模拟三体运动

C++模拟三体运动

受力计算

对于三个物体 b 1 b_1 b1, b 2 b_2 b2, b 3 b_3 b3, 其质量为 m 1 m_1 m1, m 2 m_2 m2, m 3 m_3 m3, t t t时坐标为 p 1 ⃗ \vec{p_1} p1 , p 2 ⃗ \vec{p_2} p2 , p 3 ⃗ \vec{p_3} p3 , 加在 b 1 b_1 b1上的力为 F ⃗ = G m 1 m 2 ∣ p 1 ⃗ − p 2 ⃗ ∣ 2 ⋅ ( p 2 ⃗ − p 1 ⃗ ) ^ + G m 1 m 3 ∣ p 1 ⃗ − p 3 ⃗ ∣ 2 ⋅ ( p 3 ⃗ − p 1 ⃗ ) ^ \vec{F} = \frac{G m_1 m_2}{|\vec{p_1} - \vec{p_2}|^2} \cdot \widehat{(\vec{p_2} - \vec{p_1})} + \frac{G m_1 m_3}{|\vec{p_1} - \vec{p_3}|^2} \cdot \widehat{(\vec{p_3} - \vec{p_1})} F =p1 p2 2Gm1m2(p2 p1 ) +p1 p3 2Gm1m3(p3 p1 ) .

x x x 分量为 F x = G m 1 m 2 ∣ p 1 ⃗ − p 2 ⃗ ∣ 2 ⋅ p 1 x − p 2 x ∣ p 1 ⃗ − p 2 ⃗ ∣ F_x = \frac{G m_1 m_2}{|\vec{p_1} - \vec{p_2}|^2} \cdot \frac{p_{1x} - p_{2x}}{|\vec{p_1} - \vec{p_2}|} Fx=p1 p2 2Gm1m2p1 p2 p1xp2x. 对于其 y y y, z z z 分量同理.

因此, 可用以下函数计算 b 1 b_1 b1受到的重力:

Vector g_force(float m1, Vector p1, float m2, Vector p2, float m3, Vector p3)
{
    float r = sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) + (p1.z - p2.z) * (p1.z - p2.z));
    Vector force;
    float f = G * (m1 * m2) / (r * r);
    force.x = f / r * (p2.x - p1.x);
    force.y = f / r * (p2.y - p1.y);
    force.z = f / r * (p2.z - p1.z);

    r = sqrt((p1.x - p3.x) * (p1.x - p3.x) + (p1.y - p3.y) * (p1.y - p3.y) + (p1.z - p3.z) * (p1.z - p3.z));
    f = G * (m1 * m3) / (r * r);
    force.x += f / r * (p3.x - p1.x);
    force.y += f / r * (p3.y - p1.y);
    force.z += f / r * (p3.z - p1.z);
    return force;
}

运动计算

根据受力可以算出加速度, 乘上 Δ t \Delta t Δt得到速度增量;
速度乘上 Δ t \Delta t Δt得到坐标增量.

void Body::move(Vector force, float time)
{
    Vector accel;
    accel.x = force.x / _mass;
    accel.y = force.y / _mass;
    accel.z = force.z / _mass;

    _speed.x += accel.x * time;
    _speed.y += accel.y * time;
    _speed.z += accel.z * time;

    _coordinate.x += _speed.x * time;
    _coordinate.y += _speed.y * time;
    _coordinate.z += _speed.z * time;
}

数值模拟

对每一个body执行move, 用很小的时间步长, 和实时的重力, 并且输出结果.

    for (int i = 0; i < 100; i++)
    {
        b1.move(g_force(m1, b1.get_coord(), m2, b2.get_coord(), m3, b3.get_coord()), 0.1);
        b2.move(g_force(m2, b2.get_coord(), m1, b1.get_coord(), m3, b3.get_coord()), 0.1);
        b3.move(g_force(m3, b3.get_coord(), m2, b2.get_coord(), m1, b1.get_coord()), 0.1);
        printf("Time: %f\n", (float)i / 10.0);
        output_vector(b1.get_coord(), "b1 coord");
        output_vector(b2.get_coord(), "b2 coord");
        output_vector(b3.get_coord(), "b3 coord");
        printf("\n");
    }

全部代码

body.h
#include <cstdio>
#include <cstdlib>

struct Vector
{
    float x;
    float y;
    float z;
};


class Body
{
private:
    Vector _coordinate;
    Vector _speed;
    float _mass;
public:
    Body(float mass, Vector coord);
    void move(Vector force, float time);
    Vector get_speed();
    Vector get_coord();
};

body.cpp
#include "body.h"

Body::Body(float mass, Vector coord)
{
    _coordinate = coord;
    _mass = mass;
    _speed.x = 0;
    _speed.y = 0;
    _speed.z = 0;
}

void Body::move(Vector force, float time)
{
    Vector accel;
    accel.x = force.x / _mass;
    accel.y = force.y / _mass;
    accel.z = force.z / _mass;

    _speed.x += accel.x * time;
    _speed.y += accel.y * time;
    _speed.z += accel.z * time;

    _coordinate.x += _speed.x * time;
    _coordinate.y += _speed.y * time;
    _coordinate.z += _speed.z * time;
}

Vector Body::get_speed()
{
    return _speed;
}

Vector Body::get_coord()
{
    return _coordinate;
}
main.cpp
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "body.h"

const float G = 6.67e-11;

void output_vector(Vector vec, char *comment)
{
    printf("%s: ", comment);
    printf("x: %f, y: %f, z: %f\n", vec.x, vec.y, vec.z);
}

Vector g_force(float m1, Vector p1, float m2, Vector p2, float m3, Vector p3)
{
    float r = sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) + (p1.z - p2.z) * (p1.z - p2.z));
    Vector force;
    float f = G * (m1 * m2) / (r * r);
    force.x = f / r * (p2.x - p1.x);
    force.y = f / r * (p2.y - p1.y);
    force.z = f / r * (p2.z - p1.z);

    r = sqrt((p1.x - p3.x) * (p1.x - p3.x) + (p1.y - p3.y) * (p1.y - p3.y) + (p1.z - p3.z) * (p1.z - p3.z));
    f = G * (m1 * m3) / (r * r);
    force.x += f / r * (p3.x - p1.x);
    force.y += f / r * (p3.y - p1.y);
    force.z += f / r * (p3.z - p1.z);
    return force;
}

int main()
{

    float m1 = 1e10;
    float m2 = 1e10;
    float m3 = 1e10;

    Vector p1;
    p1.x = 0;
    p1.y = 1;
    p1.z = 0;
    Body b1 = Body(m1, p1);

    Vector p2;
    p2.x = 0.866;
    p2.y = -0.5;
    p2.z = 0;
    Body b2 = Body(m2, p2);

    Vector p3;
    p3.x = -0.866;
    p3.y = -0.5;
    p3.z = 0;
    Body b3 = Body(m3, p3);


    for (int i = 0; i < 100; i++)
    {
        b1.move(g_force(m1, b1.get_coord(), m2, b2.get_coord(), m3, b3.get_coord()), 0.1);
        b2.move(g_force(m2, b2.get_coord(), m1, b1.get_coord(), m3, b3.get_coord()), 0.1);
        b3.move(g_force(m3, b3.get_coord(), m2, b2.get_coord(), m1, b1.get_coord()), 0.1);
        printf("Time: %f\n", (float)i / 10.0);
        output_vector(b1.get_coord(), "b1 coord");
        output_vector(b2.get_coord(), "b2 coord");
        output_vector(b3.get_coord(), "b3 coord");
        printf("\n");
    }
    return 0;
}

模拟结果

Time: 0.000000
b1 coord: x: 0.000000, y: 0.996149, z: 0.000000
b2 coord: x: 0.862658, y: -0.498068, z: 0.000000
b3 coord: x: -0.862650, y: -0.498066, z: 0.000000

Time: 0.100000
b1 coord: x: -0.000000, y: 0.988417, z: 0.000000
b2 coord: x: 0.855943, y: -0.494184, z: 0.000000
b3 coord: x: -0.855908, y: -0.494173, z: 0.000000

Time: 0.200000
b1 coord: x: -0.000000, y: 0.976743, z: 0.000000
b2 coord: x: 0.845793, y: -0.488308, z: 0.000000
b3 coord: x: -0.845705, y: -0.488283, z: 0.000000

Time: 0.300000
b1 coord: x: -0.000000, y: 0.961032, z: 0.000000
b2 coord: x: 0.832117, y: -0.480387, z: 0.000000
b3 coord: x: -0.831938, y: -0.480335, z: 0.000000

Time: 0.400000
b1 coord: x: -0.000000, y: 0.941150, z: 0.000000
b2 coord: x: 0.814792, y: -0.470345, z: 0.000000
b3 coord: x: -0.814470, y: -0.470250, z: 0.000000

Time: 0.500000
b1 coord: x: -0.000000, y: 0.916917, z: 0.000000
b2 coord: x: 0.793648, y: -0.458080, z: 0.000000
b3 coord: x: -0.793117, y: -0.457922, z: 0.000000

Time: 0.600000
b1 coord: x: -0.000000, y: 0.888099, z: 0.000000
b2 coord: x: 0.768468, y: -0.443462, z: 0.000000
b3 coord: x: -0.767643, y: -0.443215, z: 0.000000

Time: 0.700000
b1 coord: x: -0.000001, y: 0.854390, z: 0.000000
b2 coord: x: 0.738968, y: -0.426320, z: 0.000000
b3 coord: x: -0.737737, y: -0.425950, z: 0.000000

Time: 0.800000
b1 coord: x: -0.000002, y: 0.815392, z: 0.000000
b2 coord: x: 0.704775, y: -0.406433, z: 0.000000
b3 coord: x: -0.702991, y: -0.405891, z: 0.000000

Time: 0.900000
b1 coord: x: -0.000003, y: 0.770579, z: 0.000000
b2 coord: x: 0.665396, y: -0.383504, z: 0.000000
b3 coord: x: -0.662864, y: -0.382726, z: 0.000000

Time: 1.000000
b1 coord: x: -0.000006, y: 0.719242, z: 0.000000
b2 coord: x: 0.620159, y: -0.357130, z: 0.000000
b3 coord: x: -0.616608, y: -0.356023, z: 0.000000

Time: 1.100000
b1 coord: x: -0.000010, y: 0.660393, z: 0.000000
b2 coord: x: 0.568116, y: -0.326741, z: 0.000000
b3 coord: x: -0.563157, y: -0.325169, z: 0.000000

Time: 1.200000
b1 coord: x: -0.000018, y: 0.592594, z: 0.000000
b2 coord: x: 0.507860, y: -0.291487, z: 0.000000
b3 coord: x: -0.500890, y: -0.289228, z: 0.000000

Time: 1.300000
b1 coord: x: -0.000032, y: 0.513591, z: 0.000000
b2 coord: x: 0.437121, y: -0.249992, z: 0.000000
b3 coord: x: -0.427120, y: -0.246651, z: 0.000000

Time: 1.400000
b1 coord: x: -0.000062, y: 0.419464, z: 0.000000
b2 coord: x: 0.351748, y: -0.199727, z: 0.000000
b3 coord: x: -0.336664, y: -0.194446, z: 0.000000

Time: 1.500000
b1 coord: x: -0.000138, y: 0.301971, z: 0.000000
b2 coord: x: 0.242102, y: -0.134813, z: 0.000000
b3 coord: x: -0.216203, y: -0.124872, z: 0.000000

Time: 1.600000
b1 coord: x: -0.000406, y: 0.135095, z: 0.000000
b2 coord: x: 0.066865, y: -0.031526, z: 0.000000
b3 coord: x: 0.012879, y: 0.013173, z: 0.000000

Time: 1.700000
b1 coord: x: 0.124700, y: -0.664174, z: 0.000000
b2 coord: x: -1.152653, y: 0.921198, z: 0.000000
b3 coord: x: 0.241856, y: 0.139134, z: 0.000000

Time: 1.800000
b1 coord: x: 0.250257, y: -1.452175, z: 0.000000
b2 coord: x: -2.369449, y: 1.871891, z: 0.000000
b3 coord: x: 0.470281, y: 0.262836, z: 0.000000

Time: 1.900000
b1 coord: x: 0.375867, y: -2.237670, z: 0.000000
b2 coord: x: -3.585548, y: 2.822047, z: 0.000000
b3 coord: x: 0.698421, y: 0.385628, z: 0.000000

Time: 2.000000
b1 coord: x: 0.501495, y: -3.022091, z: 0.000000
b2 coord: x: -4.801333, y: 3.771961, z: 0.000000
b3 coord: x: 0.926391, y: 0.507933, z: 0.000000

Time: 2.100000
b1 coord: x: 0.627130, y: -3.805917, z: 0.000000
b2 coord: x: -6.016941, y: 4.721736, z: 0.000000
b3 coord: x: 1.154251, y: 0.629934, z: 0.000000

Time: 2.200000
b1 coord: x: 0.752770, y: -4.589365, z: 0.000000
b2 coord: x: -7.232433, y: 5.671422, z: 0.000000
b3 coord: x: 1.382032, y: 0.751729, z: 0.000000

Time: 2.300000
b1 coord: x: 0.878412, y: -5.372554, z: 0.000000
b2 coord: x: -8.447846, y: 6.621046, z: 0.000000
b3 coord: x: 1.609755, y: 0.873373, z: 0.000000

Time: 2.400000
b1 coord: x: 1.004056, y: -6.155553, z: 0.000000
b2 coord: x: -9.663200, y: 7.570623, z: 0.000000
b3 coord: x: 1.837433, y: 0.994904, z: 0.000000

Time: 2.500000
b1 coord: x: 1.129702, y: -6.938406, z: 0.000000
b2 coord: x: -10.878510, y: 8.520166, z: 0.000000
b3 coord: x: 2.065075, y: 1.116346, z: 0.000000

Time: 2.600000
b1 coord: x: 1.255348, y: -7.721145, z: 0.000000
b2 coord: x: -12.093782, y: 9.469681, z: 0.000000
b3 coord: x: 2.292689, y: 1.237715, z: 0.000000

Time: 2.700000
b1 coord: x: 1.380996, y: -8.503791, z: 0.000000
b2 coord: x: -13.309027, y: 10.419173, z: 0.000000
b3 coord: x: 2.520278, y: 1.359026, z: 0.000000

Time: 2.800000
b1 coord: x: 1.506644, y: -9.286361, z: 0.000000
b2 coord: x: -14.524246, y: 11.368647, z: 0.000000
b3 coord: x: 2.747847, y: 1.480288, z: 0.000000

Time: 2.900000
b1 coord: x: 1.632293, y: -10.068867, z: 0.000000
b2 coord: x: -15.739446, y: 12.318105, z: 0.000000
b3 coord: x: 2.975399, y: 1.601507, z: 0.000000

Time: 3.000000
b1 coord: x: 1.757941, y: -10.851317, z: 0.000000
b2 coord: x: -16.954628, y: 13.267550, z: 0.000000
b3 coord: x: 3.202936, y: 1.722690, z: 0.000000

Time: 3.100000
b1 coord: x: 1.883591, y: -11.633721, z: 0.000000
b2 coord: x: -18.169796, y: 14.216982, z: 0.000000
b3 coord: x: 3.430460, y: 1.843842, z: 0.000000

Time: 3.200000
b1 coord: x: 2.009241, y: -12.416084, z: 0.000000
b2 coord: x: -19.384951, y: 15.166405, z: 0.000000
b3 coord: x: 3.657972, y: 1.964967, z: 0.000000

Time: 3.300000
b1 coord: x: 2.134891, y: -13.198411, z: 0.000000
b2 coord: x: -20.600094, y: 16.115818, z: 0.000000
b3 coord: x: 3.885474, y: 2.086067, z: 0.000000

Time: 3.400000
b1 coord: x: 2.260541, y: -13.980706, z: 0.000000
b2 coord: x: -21.815228, y: 17.065224, z: 0.000000
b3 coord: x: 4.112967, y: 2.207145, z: 0.000000

Time: 3.500000
b1 coord: x: 2.386191, y: -14.762973, z: 0.000000
b2 coord: x: -23.030352, y: 18.014622, z: 0.000000
b3 coord: x: 4.340452, y: 2.328204, z: 0.000000

Time: 3.600000
b1 coord: x: 2.511842, y: -15.545215, z: 0.000000
b2 coord: x: -24.245468, y: 18.964014, z: 0.000000
b3 coord: x: 4.567930, y: 2.449246, z: 0.000000

Time: 3.700000
b1 coord: x: 2.637492, y: -16.327433, z: 0.000000
b2 coord: x: -25.460577, y: 19.913401, z: 0.000000
b3 coord: x: 4.795401, y: 2.570271, z: 0.000000

Time: 3.800000
b1 coord: x: 2.763143, y: -17.109631, z: 0.000000
b2 coord: x: -26.675680, y: 20.862782, z: 0.000000
b3 coord: x: 5.022865, y: 2.691283, z: 0.000000

Time: 3.900000
b1 coord: x: 2.888794, y: -17.891809, z: 0.000000
b2 coord: x: -27.890776, y: 21.812159, z: 0.000000
b3 coord: x: 5.250324, y: 2.812281, z: 0.000000

Time: 4.000000
b1 coord: x: 3.014446, y: -18.673971, z: 0.000000
b2 coord: x: -29.105865, y: 22.761532, z: 0.000000
b3 coord: x: 5.477778, y: 2.933267, z: 0.000000

Time: 4.100000
b1 coord: x: 3.140097, y: -19.456116, z: 0.000000
b2 coord: x: -30.320951, y: 23.710901, z: 0.000000
b3 coord: x: 5.705228, y: 3.054242, z: 0.000000

Time: 4.200000
b1 coord: x: 3.265748, y: -20.238247, z: 0.000000
b2 coord: x: -31.536032, y: 24.660267, z: 0.000000
b3 coord: x: 5.932673, y: 3.175206, z: 0.000000

Time: 4.300000
b1 coord: x: 3.391400, y: -21.020365, z: 0.000000
b2 coord: x: -32.751110, y: 25.609629, z: 0.000000
b3 coord: x: 6.160114, y: 3.296161, z: 0.000000

Time: 4.400000
b1 coord: x: 3.517051, y: -21.802469, z: 0.000000
b2 coord: x: -33.966183, y: 26.558987, z: 0.000000
b3 coord: x: 6.387551, y: 3.417107, z: 0.000000

Time: 4.500000
b1 coord: x: 3.642703, y: -22.584562, z: 0.000000
b2 coord: x: -35.181252, y: 27.508343, z: 0.000000
b3 coord: x: 6.614985, y: 3.538045, z: 0.000000

Time: 4.600000
b1 coord: x: 3.768354, y: -23.366644, z: 0.000000
b2 coord: x: -36.396317, y: 28.457695, z: 0.000000
b3 coord: x: 6.842415, y: 3.658976, z: 0.000000

Time: 4.700000
b1 coord: x: 3.894006, y: -24.148716, z: 0.000000
b2 coord: x: -37.611378, y: 29.407045, z: 0.000000
b3 coord: x: 7.069842, y: 3.779899, z: 0.000000

Time: 4.800000
b1 coord: x: 4.019658, y: -24.930779, z: 0.000000
b2 coord: x: -38.826439, y: 30.356394, z: 0.000000
b3 coord: x: 7.297266, y: 3.900815, z: 0.000000

Time: 4.900000
b1 coord: x: 4.145310, y: -25.712831, z: 0.000000
b2 coord: x: -40.041496, y: 31.305740, z: 0.000000
b3 coord: x: 7.524688, y: 4.021725, z: 0.000000

Time: 5.000000
b1 coord: x: 4.270962, y: -26.494877, z: 0.000000
b2 coord: x: -41.256550, y: 32.255085, z: 0.000000
b3 coord: x: 7.752107, y: 4.142629, z: 0.000000

Time: 5.100000
b1 coord: x: 4.396614, y: -27.276915, z: 0.000000
b2 coord: x: -42.471600, y: 33.204426, z: 0.000000
b3 coord: x: 7.979523, y: 4.263527, z: 0.000000

Time: 5.200000
b1 coord: x: 4.522265, y: -28.058945, z: 0.000000
b2 coord: x: -43.686649, y: 34.153767, z: 0.000000
b3 coord: x: 8.206938, y: 4.384420, z: 0.000000

Time: 5.300000
b1 coord: x: 4.647918, y: -28.840967, z: 0.000000
b2 coord: x: -44.901695, y: 35.103104, z: 0.000000
b3 coord: x: 8.434350, y: 4.505308, z: 0.000000

Time: 5.400000
b1 coord: x: 4.773570, y: -29.622982, z: 0.000000
b2 coord: x: -46.116741, y: 36.052441, z: 0.000000
b3 coord: x: 8.661760, y: 4.626191, z: 0.000000

Time: 5.500000
b1 coord: x: 4.899222, y: -30.404991, z: 0.000000
b2 coord: x: -47.331783, y: 37.001774, z: 0.000000
b3 coord: x: 8.889169, y: 4.747069, z: 0.000000

Time: 5.600000
b1 coord: x: 5.024875, y: -31.186995, z: 0.000000
b2 coord: x: -48.546825, y: 37.951107, z: 0.000000
b3 coord: x: 9.116575, y: 4.867944, z: 0.000000

Time: 5.700000
b1 coord: x: 5.150527, y: -31.968992, z: 0.000000
b2 coord: x: -49.761864, y: 38.900440, z: 0.000000
b3 coord: x: 9.343980, y: 4.988814, z: 0.000000

Time: 5.800000
b1 coord: x: 5.276179, y: -32.750984, z: 0.000000
b2 coord: x: -50.976902, y: 39.849770, z: 0.000000
b3 coord: x: 9.571383, y: 5.109680, z: 0.000000

Time: 5.900000
b1 coord: x: 5.401832, y: -33.532970, z: 0.000000
b2 coord: x: -52.191936, y: 40.799099, z: 0.000000
b3 coord: x: 9.798783, y: 5.230543, z: 0.000000

Time: 6.000000
b1 coord: x: 5.527484, y: -34.314953, z: 0.000000
b2 coord: x: -53.406971, y: 41.748428, z: 0.000000
b3 coord: x: 10.026183, y: 5.351401, z: 0.000000

Time: 6.100000
b1 coord: x: 5.653136, y: -35.096931, z: 0.000000
b2 coord: x: -54.622002, y: 42.697754, z: 0.000000
b3 coord: x: 10.253581, y: 5.472257, z: 0.000000

Time: 6.200000
b1 coord: x: 5.778789, y: -35.878902, z: 0.000000
b2 coord: x: -55.837032, y: 43.647079, z: 0.000000
b3 coord: x: 10.480978, y: 5.593109, z: 0.000000

Time: 6.300000
b1 coord: x: 5.904441, y: -36.660870, z: 0.000000
b2 coord: x: -57.052063, y: 44.596405, z: 0.000000
b3 coord: x: 10.708373, y: 5.713958, z: 0.000000

Time: 6.400000
b1 coord: x: 6.030093, y: -37.442833, z: 0.000000
b2 coord: x: -58.267090, y: 45.545727, z: 0.000000
b3 coord: x: 10.935767, y: 5.834804, z: 0.000000

Time: 6.500000
b1 coord: x: 6.155746, y: -38.224792, z: 0.000000
b2 coord: x: -59.482117, y: 46.495049, z: 0.000000
b3 coord: x: 11.163160, y: 5.955647, z: 0.000000

Time: 6.600000
b1 coord: x: 6.281398, y: -39.006748, z: 0.000000
b2 coord: x: -60.697144, y: 47.444370, z: 0.000000
b3 coord: x: 11.390552, y: 6.076488, z: 0.000000

Time: 6.700000
b1 coord: x: 6.407050, y: -39.788700, z: 0.000000
b2 coord: x: -61.912167, y: 48.393692, z: 0.000000
b3 coord: x: 11.617942, y: 6.197326, z: 0.000000

Time: 6.800000
b1 coord: x: 6.532702, y: -40.570648, z: 0.000000
b2 coord: x: -63.127190, y: 49.343010, z: 0.000000
b3 coord: x: 11.845331, y: 6.318161, z: 0.000000

Time: 6.900000
b1 coord: x: 6.658355, y: -41.352592, z: 0.000000
b2 coord: x: -64.342209, y: 50.292328, z: 0.000000
b3 coord: x: 12.072720, y: 6.438994, z: 0.000000

Time: 7.000000
b1 coord: x: 6.784008, y: -42.134537, z: 0.000000
b2 coord: x: -65.557228, y: 51.241646, z: 0.000000
b3 coord: x: 12.300107, y: 6.559824, z: 0.000000

Time: 7.100000
b1 coord: x: 6.909661, y: -42.916477, z: 0.000000
b2 coord: x: -66.772247, y: 52.190964, z: 0.000000
b3 coord: x: 12.527493, y: 6.680653, z: 0.000000

Time: 7.200000
b1 coord: x: 7.035314, y: -43.698414, z: 0.000000
b2 coord: x: -67.987267, y: 53.140282, z: 0.000000
b3 coord: x: 12.754879, y: 6.801478, z: 0.000000

Time: 7.300000
b1 coord: x: 7.160966, y: -44.480347, z: 0.000000
b2 coord: x: -69.202286, y: 54.089596, z: 0.000000
b3 coord: x: 12.982264, y: 6.922302, z: 0.000000

Time: 7.400000
b1 coord: x: 7.286619, y: -45.262276, z: 0.000000
b2 coord: x: -70.417305, y: 55.038910, z: 0.000000
b3 coord: x: 13.209647, y: 7.043124, z: 0.000000

Time: 7.500000
b1 coord: x: 7.412272, y: -46.044205, z: 0.000000
b2 coord: x: -71.632324, y: 55.988224, z: 0.000000
b3 coord: x: 13.437030, y: 7.163944, z: 0.000000

Time: 7.600000
b1 coord: x: 7.537925, y: -46.826130, z: 0.000000
b2 coord: x: -72.847343, y: 56.937538, z: 0.000000
b3 coord: x: 13.664412, y: 7.284762, z: 0.000000

Time: 7.700000
b1 coord: x: 7.663578, y: -47.608051, z: 0.000000
b2 coord: x: -74.062355, y: 57.886852, z: 0.000000
b3 coord: x: 13.891792, y: 7.405579, z: 0.000000

Time: 7.800000
b1 coord: x: 7.789230, y: -48.389973, z: 0.000000
b2 coord: x: -75.277367, y: 58.836166, z: 0.000000
b3 coord: x: 14.119173, y: 7.526393, z: 0.000000

Time: 7.900000
b1 coord: x: 7.914883, y: -49.171890, z: 0.000000
b2 coord: x: -76.492378, y: 59.785477, z: 0.000000
b3 coord: x: 14.346553, y: 7.647205, z: 0.000000

Time: 8.000000
b1 coord: x: 8.040536, y: -49.953808, z: 0.000000
b2 coord: x: -77.707390, y: 60.734787, z: 0.000000
b3 coord: x: 14.573932, y: 7.768016, z: 0.000000

Time: 8.100000
b1 coord: x: 8.166188, y: -50.735722, z: 0.000000
b2 coord: x: -78.922401, y: 61.684097, z: 0.000000
b3 coord: x: 14.801310, y: 7.888825, z: 0.000000

Time: 8.200000
b1 coord: x: 8.291841, y: -51.517632, z: 0.000000
b2 coord: x: -80.137413, y: 62.633408, z: 0.000000
b3 coord: x: 15.028687, y: 8.009633, z: 0.000000

Time: 8.300000
b1 coord: x: 8.417494, y: -52.299541, z: 0.000000
b2 coord: x: -81.352425, y: 63.582718, z: 0.000000
b3 coord: x: 15.256064, y: 8.130439, z: 0.000000

Time: 8.400000
b1 coord: x: 8.543147, y: -53.081448, z: 0.000000
b2 coord: x: -82.567436, y: 64.532028, z: 0.000000
b3 coord: x: 15.483440, y: 8.251244, z: 0.000000

Time: 8.500000
b1 coord: x: 8.668800, y: -53.863354, z: 0.000000
b2 coord: x: -83.782448, y: 65.481339, z: 0.000000
b3 coord: x: 15.710816, y: 8.372046, z: 0.000000

Time: 8.600000
b1 coord: x: 8.794454, y: -54.645256, z: 0.000000
b2 coord: x: -84.997459, y: 66.430649, z: 0.000000
b3 coord: x: 15.938191, y: 8.492848, z: 0.000000

Time: 8.700000
b1 coord: x: 8.920107, y: -55.427158, z: 0.000000
b2 coord: x: -86.212471, y: 67.379959, z: 0.000000
b3 coord: x: 16.165565, y: 8.613648, z: 0.000000

Time: 8.800000
b1 coord: x: 9.045760, y: -56.209057, z: 0.000000
b2 coord: x: -87.427475, y: 68.329269, z: 0.000000
b3 coord: x: 16.392939, y: 8.734447, z: 0.000000

Time: 8.900000
b1 coord: x: 9.171413, y: -56.990955, z: 0.000000
b2 coord: x: -88.642479, y: 69.278580, z: 0.000000
b3 coord: x: 16.620312, y: 8.855246, z: 0.000000

Time: 9.000000
b1 coord: x: 9.297067, y: -57.772850, z: 0.000000
b2 coord: x: -89.857483, y: 70.227882, z: 0.000000
b3 coord: x: 16.847685, y: 8.976042, z: 0.000000

Time: 9.100000
b1 coord: x: 9.422720, y: -58.554745, z: 0.000000
b2 coord: x: -91.072487, y: 71.177185, z: 0.000000
b3 coord: x: 17.075056, y: 9.096837, z: 0.000000

Time: 9.200000
b1 coord: x: 9.548373, y: -59.336639, z: 0.000000
b2 coord: x: -92.287491, y: 72.126488, z: 0.000000
b3 coord: x: 17.302427, y: 9.217631, z: 0.000000

Time: 9.300000
b1 coord: x: 9.674026, y: -60.118530, z: 0.000000
b2 coord: x: -93.502495, y: 73.075790, z: 0.000000
b3 coord: x: 17.529799, y: 9.338425, z: 0.000000

Time: 9.400000
b1 coord: x: 9.799680, y: -60.900421, z: 0.000000
b2 coord: x: -94.717499, y: 74.025093, z: 0.000000
b3 coord: x: 17.757170, y: 9.459216, z: 0.000000

Time: 9.500000
b1 coord: x: 9.925333, y: -61.682308, z: 0.000000
b2 coord: x: -95.932503, y: 74.974396, z: 0.000000
b3 coord: x: 17.984539, y: 9.580007, z: 0.000000

Time: 9.600000
b1 coord: x: 10.050986, y: -62.464195, z: 0.000000
b2 coord: x: -97.147507, y: 75.923698, z: 0.000000
b3 coord: x: 18.211908, y: 9.700796, z: 0.000000

Time: 9.700000
b1 coord: x: 10.176640, y: -63.246082, z: 0.000000
b2 coord: x: -98.362511, y: 76.873001, z: 0.000000
b3 coord: x: 18.439278, y: 9.821585, z: 0.000000

Time: 9.800000
b1 coord: x: 10.302293, y: -64.027969, z: 0.000000
b2 coord: x: -99.577515, y: 77.822304, z: 0.000000
b3 coord: x: 18.666647, y: 9.942372, z: 0.000000

Time: 9.900000
b1 coord: x: 10.427946, y: -64.809853, z: 0.000000
b2 coord: x: -100.792519, y: 78.771606, z: 0.000000
b3 coord: x: 18.894014, y: 10.063159, z: 0.000000

可视化: 略.

结论

通过数值模拟, 可以得出三体运动是混沌的.

  • 2
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
三体运动是指三个物体之间的相互作用和运动,可以使用Python进行模拟。其中,可以使用scipy库中的odeint函数解微分方程组,来模拟三体运动。 以下是一个简单的三体运动模拟代码示例: ```python import numpy as np from scipy.integrate import odeint # 定义微分方程组 def three_body_equations(w, t, G, m1, m2, m3): x1, y1, vx1, vy1, x2, y2, vx2, vy2, x3, y3, vx3, vy3 = w r12 = np.sqrt((x2 - x1)**2 + (y2 - y1)**2) r13 = np.sqrt((x3 - x1)**2 + (y3 - y1)**2) r23 = np.sqrt((x3 - x2)**2 + (y3 - y2)**2) dx1dt = vx1 dy1dt = vy1 dvx1dt = G * m2 * (x2 - x1) / r12**3 + G * m3 * (x3 - x1) / r13**3 dvy1dt = G * m2 * (y2 - y1) / r12**3 + G * m3 * (y3 - y1) / r13**3 dx2dt = vx2 dy2dt = vy2 dvx2dt = G * m1 * (x1 - x2) / r12**3 + G * m3 * (x3 - x2) / r23**3 dvy2dt = G * m1 * (y1 - y2) / r12**3 + G * m3 * (y3 - y2) / r23**3 dx3dt = vx3 dy3dt = vy3 dvx3dt = G * m1 * (x1 - x3) / r13**3 + G * m2 * (x2 - x3) / r23**3 dvy3dt = G * m1 * (y1 - y3) / r13**3 + G * m2 * (y2 - y3) / r23**3 return dx1dt, dy1dt, dvx1dt, dvy1dt, dx2dt, dy2dt, dvx2dt, dvy2dt, dx3dt, dy3dt, dvx3dt, dvy3dt # 定义初始状态和参数 w0 = [1, 0, 0, 6, -1, 0, 0, -6, 0, 0, 0, 0] t = np.linspace(0, 10, 1000) G = 1 m1 = 1 m2 = 1 m3 = 1 # 解微分方程组 wsol = odeint(three_body_equations, w0, t, args=(G, m1, m2, m3)) # 绘制轨迹图 import matplotlib.pyplot as plt plt.plot(wsol[:,0], wsol[:,1], label='Body 1') plt.plot(wsol[:,4], wsol[:,5], label='Body 2') plt.plot(wsol[:,8], wsol[:,9], label='Body 3') plt.legend() plt.show() ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值