如何在Python代码中缩小轨道半径实现双星碰撞？

问题根源分析 #

你的代码存在几个核心问题导致轨道无法衰减甚至脱轨：

1. 阻力方向错误：当前代码中阻力直接叠加在引力项上，未与速度方向绑定，相当于施加了固定方向的力，而非阻碍运动的阻力。
2. 引力模型不完整：仅考虑中心天体mS对两个天体的引力，未加入双星之间的相互引力，不符合双星系统的物理规律。
3. 位置与引力计算时序错误：计算引力时使用旧的位置变量（如x1），而非当前时刻的最新位置（xlist1[-1]），导致引力计算偏差。
4. 手动修改坐标违背物理逻辑：直接给坐标减固定值会破坏轨道动力学平衡，引发运动混乱。

具体修正步骤 #

1. 修正阻力公式 #

阻力加速度应与速度方向相反，定义阻力系数后，阻力项表达式为drag_coeff * v（v为当前速度分量），方向与速度相反，需在速度更新时减去该阻力带来的加速度（除以质量后）。

2. 添加双星间的相互引力 #

计算两个天体之间的距离，分别给双方施加对方的引力作用，这是双星系统的核心交互逻辑。

3. 修正引力计算的位置变量 #

使用每个天体当前时刻的位置（列表最后一个元素）计算引力，确保力的计算与位置匹配。

4. 调整初始参数匹配物理规律 #

确保初始速度满足圆周运动的向心力条件，施加阻力后天体才会逐渐向内螺旋运动。

修改后的代码 #

import matplotlib.pyplot as plt

G = 6.67e-11  # 引力常数
dt = 1e-3     # 时间步长
mS = 1e6      # 中心天体质量（可按需移除）
mp1 = 1e5     # 天体1质量
mp2 = 1e5     # 天体2质量
drag_coeff = 1e-4  # 阻力系数，可根据衰减速度调整

# 初始位置：双星对称分布在原点两侧
x1 = 1
y1 = 0
x2 = -1
y2 = 0

# 初始速度：满足绕原点圆周运动的横向速度，适配双星互绕
vx1 = 0
vy1 = 1.5
vx2 = 0
vy2 = -1.5

# 绘图样式设置
for param in ['figure.facecolor', 'axes.facecolor', 'savefig.facecolor']:
    plt.rcParams[param] = 'black'
plt.style.use("cyberpunk")

timestep = 500  # 视频帧步数
n = 10          # 每帧内部迭代次数
tail_length = 100

# 初始化位置和速度列表
xlist1 = [x1]
ylist1 = [y1]
vxlist1 = [vx1]
vylist1 = [vy1]

xlist2 = [x2]
ylist2 = [y2]
vxlist2 = [vx2]
vylist2 = [vy2]

for t in range(1, timestep+1):
    for _ in range(n):
        # 获取当前时刻的位置和速度
        x1_current = xlist1[-1]
        y1_current = ylist1[-1]
        vx1_current = vxlist1[-1]
        vy1_current = vylist1[-1]
        
        x2_current = xlist2[-1]
        y2_current = ylist2[-1]
        vx2_current = vxlist2[-1]
        vy2_current = vylist2[-1]
        
        # 计算天体1受到的合力：中心天体引力 + 天体2引力 + 阻力
        # 中心天体引力
        r1_to_S = (x1_current**2 + y1_current**2)**0.5
        acc_x_S1 = -G * mS * x1_current / r1_to_S**3
        acc_y_S1 = -G * mS * y1_current / r1_to_S**3
        
        # 天体2对天体1的引力
        dx = x2_current - x1_current
        dy = y2_current - y1_current
        r12 = (dx**2 + dy**2)**0.5
        acc_x_21 = G * mp2 * dx / r12**3
        acc_y_21 = G * mp2 * dy / r12**3
        
        # 阻力：与速度方向相反
        acc_x_drag1 = -drag_coeff * vx1_current / mp1
        acc_y_drag1 = -drag_coeff * vy1_current / mp1
        
        # 更新天体1的速度和位置
        new_vx1 = vx1_current + (acc_x_S1 + acc_x_21 + acc_x_drag1) * dt
        new_vy1 = vy1_current + (acc_y_S1 + acc_y_21 + acc_y_drag1) * dt
        new_x1 = x1_current + new_vx1 * dt
        new_y1 = y1_current + new_vy1 * dt
        
        xlist1.append(new_x1)
        ylist1.append(new_y1)
        vxlist1.append(new_vx1)
        vylist1.append(new_vy1)
        
        # 计算天体2受到的合力：中心天体引力 + 天体1引力 + 阻力
        # 中心天体引力
        r2_to_S = (x2_current**2 + y2_current**2)**0.5
        acc_x_S2 = -G * mS * x2_current / r2_to_S**3
        acc_y_S2 = -G * mS * y2_current / r2_to_S**3
        
        # 天体1对天体2的引力（与21大小相等方向相反）
        acc_x_12 = -acc_x_21 * mp1 / mp2
        acc_y_12 = -acc_y_21 * mp1 / mp2
        
        # 阻力
        acc_x_drag2 = -drag_coeff * vx2_current / mp2
        acc_y_drag2 = -drag_coeff * vy2_current / mp2
        
        # 更新天体2的速度和位置
        new_vx2 = vx2_current + (acc_x_S2 + acc_x_12 + acc_x_drag2) * dt
        new_vy2 = vy2_current + (acc_y_S2 + acc_y_12 + acc_y_drag2) * dt
        new_x2 = x2_current + new_vx2 * dt
        new_y2 = y2_current + new_vy2 * dt
        
        xlist2.append(new_x2)
        ylist2.append(new_y2)
        vxlist2.append(new_vx2)
        vylist2.append(new_vy2)

# 绘图示例（可结合原视频生成逻辑补充）
plt.plot(xlist1, ylist1, color='cyan', linewidth=1)
plt.plot(xlist2, ylist2, color='magenta', linewidth=1)
plt.axis('equal')
plt.show()

补充说明 #

• 阻力系数drag_coeff可根据轨道衰减快慢调整，值越大衰减速度越快。
• 若不需要中心天体mS，可直接移除其引力相关代码，变为纯粹的双星互绕系统。
• 初始速度需根据天体质量和初始位置微调，确保初始轨道稳定，否则加阻力后仍可能出现异常运动。

内容的提问来源于stack exchange，提问作者Woobzieer