生成代码 限制性三体问题

以下是 Python 中实现限制性三体问题的基本代码：

import numpy as np
import matplotlib.pyplot as plt

# 定义常数
G = 6.6743e-11  # 万有引力常数
m1 = 5.974e24  # 地球质量
m2 = 7.348e22  # 月球质量
d = 384400000  # 地球和月球之间的距离

# 定义初始条件
x1 = 0
y1 = 0
vx1 = 0
vy1 = 0
x2 = d
y2 = 0
vx2 = 0
vy2 = 1022

# 定义时间步长和总时间
dt = 1000
t = np.arange(0, 28 * 24 * 3600, dt)

# 定义数组来保存位置和速度
x1_arr = np.zeros(len(t))
y1_arr = np.zeros(len(t))
vx1_arr = np.zeros(len(t))
vy1_arr = np.zeros(len(t))
x2_arr = np.zeros(len(t))
y2_arr = np.zeros(len(t))
vx2_arr = np.zeros(len(t))
vy2_arr = np.zeros(len(t))

# 初始位置和速度
x1_arr[0] = x1
y1_arr[0] = y1
vx1_arr[0] = vx1
vy1_arr[0] = vy1
x2_arr[0] = x2
y2_arr[0] = y2
vx2_arr[0] = vx2
vy2_arr[0] = vy2

# 迭代计算位置和速度
for i in range(1, len(t)):
    # 计算地球和月球之间的距离
    r = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
    
    # 计算加速度
    ax1 = G * m2 * (x2 - x1) / r ** 3
    ay1 = G * m2 * (y2 - y1) / r ** 3
    ax2 = -G * m1 * (x2 - x1) / r ** 3
    ay2 = -G * m1 * (y2 - y1) / r ** 3
    
    # 计算速度
    vx1 += ax1 * dt
    vy1 += ay1 * dt
    vx2 += ax2 * dt
    vy2 += ay2 * dt
    
    # 计算位置
    x1 += vx1 * dt
    y1 += vy1 * dt
    x2 += vx2 * dt
    y2 += vy2 * dt
    
    # 保存位置和速度
    x1_arr[i] = x1
    y1_arr[i] = y1
    vx1_arr[i] = vx1
    vy1_arr[i] = vy1
    x2_arr[i] = x2
    y2_arr[i] = y2
    vx2_arr[i] = vx2
    vy2_arr[i] = vy2

# 绘制轨迹
plt.plot(x1_arr, y1_arr, label='Earth')
plt.plot(x2_arr, y2_arr, label='Moon')
plt.legend()
plt.show()

这个代码使用了欧拉方法来迭代计算位置和速度，并且使用了 NumPy 和 Matplotlib 库来处理数组和绘制图表。你可以根据需要进行修改，以适应你的具体应用场景。