三体运动模拟程序优化版（小心CPU爆炸）

北大西洋公约 · NATO

2025-08-20 19:22:30

发布于：广东

26阅读

0回复

0点赞

以下是一个计算/模拟三体运动的C++程序，运行可能会使CPU升温，但温度仍在CPU可承受范围内，不会对CPU造成损害。

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <ctime>
using namespace std;

// 简化为double类型，减少计算压力
typedef double real;

struct Vec3 {
    real x, y, z;
    Vec3(real x=0, real y=0, real z=0) : x(x), y(y), z(z) {}
    Vec3 operator+(const Vec3& o) const { return Vec3(x+o.x, y+o.y, z+o.z); }
    Vec3 operator-(const Vec3& o) const { return Vec3(x-o.x, y-o.y, z-o.z); }
    Vec3 operator*(real s) const { return Vec3(x*s, y*s, z*s); }
    real dot(const Vec3& o) const { return x*o.x + y*o.y + z*o.z; }
    real norm() const { return sqrt(dot(*this)); }
};

struct Body {
    Vec3 pos, vel, acc;
    real mass;
    Body(Vec3 p, Vec3 v, real m) : pos(p), vel(v), acc(Vec3()), mass(m) {}
};

struct ThreeBodySystem {
    Body bodies[3];
    real time = 0;
    const real G = 6.67430e-11;

    ThreeBodySystem() {
        // 固定初始化（避免随机数问题）
        bodies[0] = Body(Vec3(0, 0, 0), Vec3(0, 0, 0), 1e30);
        bodies[1] = Body(Vec3(100, 0, 0), Vec3(0, 10, 0), 1e30);
        bodies[2] = Body(Vec3(-50, 86.6, 0), Vec3(0, -5, 0), 1e30);
    }

    void compute_acc() {
        for (int i=0; i<3; i++) {
            bodies[i].acc = Vec3();
            for (int j=0; j<3; j++) {
                if (i != j) {
                    Vec3 r = bodies[j].pos - bodies[i].pos;
                    real r_norm = r.norm();
                    if (r_norm > 1e-9) {
                        bodies[i].acc += r * (G * bodies[j].mass / (r_norm * r_norm * r_norm));
                    }
                }
            }
        }
    }

    void rk4_step(real dt) {
        Body temp[3] = {bodies[0], bodies[1], bodies[2]};
        Vec3 k1v[3], k1r[3], k2v[3], k2r[3], k3v[3], k3r[3], k4v[3], k4r[3];

        compute_acc();
        for (int i=0; i<3; i++) { k1v[i] = bodies[i].acc; k1r[i] = bodies[i].vel; }

        for (int i=0; i<3; i++) { bodies[i].pos = temp[i].pos + k1r[i]*(dt/2); bodies[i].vel = temp[i].vel + k1v[i]*(dt/2); }
        compute_acc();
        for (int i=0; i<3; i++) { k2v[i] = bodies[i].acc; k2r[i] = bodies[i].vel; }

        for (int i=0; i<3; i++) { bodies[i].pos = temp[i].pos + k2r[i]*(dt/2); bodies[i].vel = temp[i].vel + k2v[i]*(dt/2); }
        compute_acc();
        for (int i=0; i<3; i++) { k3v[i] = bodies[i].acc; k3r[i] = bodies[i].vel; }

        for (int i=0; i<3; i++) { bodies[i].pos = temp[i].pos + k3r[i]*dt; bodies[i].vel = temp[i].vel + k3v[i]*dt; }
        compute_acc();
        for (int i=0; i<3; i++) { k4v[i] = bodies[i].acc; k4r[i] = bodies[i].vel; }

        for (int i=0; i<3; i++) {
            bodies[i].pos = temp[i].pos + (k1r[i] + k2r[i]*2 + k3r[i]*2 + k4r[i])*(dt/6);
            bodies[i].vel = temp[i].vel + (k1v[i] + k2v[i]*2 + k3v[i]*2 + k4v[i])*(dt/6);
        }
        time += dt;
    }

    void run(real total_t, real dt) {
        int steps = total_t / dt;
        for (int s=0; s<steps; s++) {
            rk4_step(dt);
            if (s % 100 == 0) {
                cout << "Time: " << fixed << setprecision(2) << time 
                     << " | Pos1: (" << bodies[0].pos.x << "," << bodies[0].pos.y << ")\r" << flush;
            }
        }
    }
};

int main() {
    cout << "基础版三体模拟（单线程）" << endl;
    ThreeBodySystem system;
    system.run(1000, 0.1); // 总时间1000，步长0.1（计算量小）
    cout << endl << "模拟结束" << endl;
    return 0;
}