编辑代码

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 常量设置
G = 6.67430e-11  # 引力常数 (m³ kg⁻¹ s⁻²)
dt = 3600 * 24   # 时间步长 (1天)
total_time = 3600 * 24 * 365 * 0.1  # 总时长 (约36.5天)

# 天体初始条件 [质量 (kg), x (m), y (m), vx (m/s), vy (m/s)]
bodies = np.array([
    [1.9885e30, 0.0, 0.0, 0.0, 0.0],        # 中心恒星
        [5.972e24, 1.496e11, 0.0, 0.0, 29.78e3],  # 类地球天体
            [5.972e24, -1.496e11, 0.0, 0.0, -29.78e3] # 反向天体
            ])

            def acceleration(bodies):
                n = len(bodies)
                    acc = np.zeros((n, 2))
                        for i in range(n):
                                for j in range(n):
                                            if i != j:
                                                            dx = bodies[j,1] - bodies[i,1]
                                                                            dy = bodies[j,2] - bodies[i,2]
                                                                                            r = np.sqrt(dx**2 + dy**2)
                                                                                                            F = G * bodies[j,0] / r**3
                                                                                                                            acc[i,0] += F * dx
                                                                                                                                            acc[i,1] += F * dy
                                                                                                                                                return acc

                                                                                                                                                def rk4_step(bodies, dt):
                                                                                                                                                    k1v = acceleration(bodies)
                                                                                                                                                        k1r = bodies[:,3:5]
                                                                                                                                                            
                                                                                                                                                                k2v = acceleration(bodies + np.hstack([np.zeros((3,1)), k1r*dt/2, np.zeros((3,2))]))
                                                                                                                                                                    k2r = bodies[:,3:5] + k1v*dt/2
                                                                                                                                                                        
                                                                                                                                                                            k3v = acceleration(bodies + np.hstack([np.zeros((3,1)), k2r*dt/2, np.zeros((3,2))]))
                                                                                                                                                                                k3r = bodies[:,3:5] + k2v*dt/2
                                                                                                                                                                                    
                                                                                                                                                                                        k4v = acceleration(bodies + np.hstack([np.zeros((3,1)), k3r*dt, np.zeros((3,2))]))
                                                                                                                                                                                            k4r = bodies[:,3:5] + k3v*dt

                                                                                                                                                                                                bodies[:,1:3] += dt/6 * (k1r + 2*k2r + 2*k3r + k4r)
                                                                                                                                                                                                    bodies[:,3:5] += dt/6 * (k1v + 2*k2v + 2*k3v + k4v)
                                                                                                                                                                                                        return bodies

                                                                                                                                                                                                        # 初始化绘图
                                                                                                                                                                                                        fig, ax = plt.subplots(figsize=(8,8))
                                                                                                                                                                                                        ax.set_xlim(-3e11, 3e11)
                                                                                                                                                                                                        ax.set_ylim(-3e11, 3e11)
                                                                                                                                                                                                        colors = ['gold', 'blue', 'red']
                                                                                                                                                                                                        traces = [ax.plot([], [], '-', color=c, alpha=0.5)[0] for c in colors]
                                                                                                                                                                                                        points = [ax.plot([], [], 'o', color=c)[0] for c in colors]

                                                                                                                                                                                                        # 存储轨迹数据
                                                                                                                                                                                                        history = [[bodies[i,1], bodies[i,2]] for i in range(3)]

                                                                                                                                                                                                        def update(frame):
                                                                                                                                                                                                            global bodies
                                                                                                                                                                                                                bodies = rk4_step(bodies, dt)
                                                                                                                                                                                                                    
                                                                                                                                                                                                                        for i in range(3):
                                                                                                                                                                                                                                history[i][0].append(bodies[i,1])
                                                                                                                                                                                                                                        history[i][1].append(bodies[i,2])
                                                                                                                                                                                                                                                traces[i].set_data(history[i][0][-1000:], history[i][1][-1000:])
                                                                                                                                                                                                                                                        points[i].set_data(bodies[i,1], bodies[i,2])
                                                                                                                                                                                                                                                            
                                                                                                                                                                                                                                                                return traces + points

                                                                                                                                                                                                                                                                ani = FuncAnimation(fig, update, frames=range(500), interval=20, blit=True)
                                                                                                                                                                                                                                                                plt.show()