地球与超级木星的引力模拟

简介

在本项目中，我们将使用 Python 开发一个引力模拟程序，展示地球与一颗假设的“超级木星”之间的相互作用。“超级木星”的质量是木星的 500 倍。该模拟旨在展示如此巨大的天体对地球运动的影响，考虑到其中起作用的巨大引力。本项目适合对物理、天文学和计算模拟充满热情的学生和爱好者。为此，我们将使用 NumPy 等 Python 库进行数值计算，并使用 Matplotlib 来可视化行星的动态运动。

👀 预览

🎯 任务

在本项目中，你将学习：

• 如何理解并应用牛顿万有引力定律来模拟天体之间的相互作用。
• 如何使用 Python 编程创建引力系统的计算模型。
• 如何在 Python 中使用 NumPy 库进行高效的数值计算。
• 如何在存在质量为木星 500 倍的“超级木星”的情况下模拟地球的轨道力学。
• 如何分析和解释模拟结果，以了解巨大天体对轨道动力学的影响。
• 如何使用 Matplotlib 实现模拟的可视化表示，展示行星的轨道路径和相对位置。
• 如何在宇宙背景下探索力、质量和加速度的概念。
• 如何针对不同场景微调模拟参数，如质量、距离和时间步长。
• 如何培养调试和优化用于科学计算的 Python 代码的技能。

🏆 成果

完成本项目后，你将能够：

• 在实际的计算环境中应用物理的基本原理，特别是牛顿万有引力定律。
• 使用 Python 创建并运行基于物理的模拟。
• 熟练使用 NumPy 高效处理大规模数值计算。
• 使用 Matplotlib 可视化复杂数据和模拟，增强科学结果的可解释性。
• 理解行星运动的动力学以及巨大天体引力的影响。
• 分析和解释模拟结果，得出关于天体力学的有意义结论。
• 调整和试验模拟参数，从而更深入地理解轨道力学。
• 在编程环境中，特别是在科学计算的背景下，展示增强的问题解决和调试技能。
• 展示关于引力如何塑造天体运动的基础知识，为进一步探索天体物理学和计算建模铺平道路。

设置环境

首先，我们需要设置项目文件并导入必要的库。

创建一个名为 simulation.py 的 Python 文件：

touch ~/project/simulation.py

然后导入以下库：

## 导入所需的库
from typing import Tuple
import numpy as np
import pylab as py
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.lines import Line2D
from tqdm import trange

在引力模拟项目的这一初始步骤中，导入了几个重要的 Python 库来为我们的模拟设置环境：

1. typing：具体来说，从 typing 模块中导入了 Tuple。这用于 Python 中的类型提示，有助于提高代码的可读性和调试性。它允许指定可以在元组中分组的数值类型，确保函数定义的一致性和清晰度。
2. numpy：以 np 导入的 numpy 库是 Python 中科学计算的基础包。它支持大型多维数组和矩阵，并提供大量用于对这些数组进行操作的高级数学函数。对于处理模拟中的数值计算至关重要。
3. pylab：以 py 导入的 pylab 是 Matplotlib 中的一个模块，与 Matplotlib 一起安装。它提供了类似 MATLAB 的接口，对于交互式计算和绘图特别有用。然而，一般不鼓励使用它，而是建议显式导入 matplotlib.pyplot 和 numpy。
4. matplotlib.pyplot：这个模块以 plt 导入，用于在 Python 中创建静态、动画和交互式可视化。它是本项目绘制天体轨迹的核心部分。
5. matplotlib.animation：Matplotlib 的这个子模块提供了创建动画的功能。对于在引力模拟中可视化行星的动态运动至关重要。
6. matplotlib.lines.Line2D：Matplotlib 库中的这个类用于创建可以绘制在图表上的线条对象。对于模拟中更详细或定制的绘图很有用。
7. tqdm.trange：tqdm 模块中的 trange 是 Python 的 range 函数的一个变体，带有内置的进度条。它对于显示模拟中循环的进度很有用，特别是在处理大量计算时。

通过导入这些库和模块，我们为引力模拟项目中所需的复杂数值计算和可视化奠定了基础。此设置对于后续实现实际物理和动画逻辑的步骤至关重要。

定义常量

定义将在我们的模拟中使用的常量。这包括引力常数、天文单位以及地球、木星和太阳的归一化质量。

## 常量
G = 6.673e-11                 ## 引力常数
AU = 1.496e11                 ## 以千米为单位的天文单位
YEAR = 365*24*60*60.0         ## 一年中的秒数
MM = 6e24                     ## 归一化质量
ME = 6e24/MM                  ## 地球的归一化质量
MS = 2e30/MM                  ## 太阳的归一化质量
MJ = 500*1.9e27/MM            ## 木星的归一化质量
GG = (MM*G*YEAR**2)/(AU**3)   ## 模拟用的引力常数

在这一步中，我们为引力模拟定义了几个基本常量。这些常量对于在模拟中准确建模天体力学和引力相互作用至关重要：

1. 引力常数（G）：设置为 6.673e-11（国际单位制），这个常数在牛顿万有引力定律中至关重要。它代表引力的强度，用于计算两个质量之间的引力。
2. 天文单位（AU）：定义为 1.496e11 千米，它表示地球到太阳的平均距离。这个单位用于以更易于管理的方式表示太阳系内的距离。
3. 年（YEAR）：计算为 365*24*60*60.0，它表示一年中的秒数。这用于将与时间相关的计算转换为年，这是天文模拟中更直观的时间单位。
4. 归一化质量（MM）：设置为 6e24，这个值用作参考质量，用于在模拟中归一化其他质量。它大致相当于地球的质量。
5. 地球的归一化质量（ME）：计算为 6e24/MM，它表示以归一化单位表示的地球质量。由于归一化质量（MM）基于地球质量，所以它实际上被设置为 1。
6. 太阳的归一化质量（MS）：设置为 2e30/MM，它表示以归一化单位表示的太阳质量。太阳质量是模拟行星引力作用的关键因素。
7. 木星的归一化质量（MJ）：在这个模拟中，独特的是，木星的质量被放大到 500*1.9e27/MM 以表示“超级木星”的情况。这个大得多的质量将展示一个质量大得多的行星对地球轨道的影响。
8. 模拟用的引力常数（GG）：这是根据模拟的尺度和单位调整后的引力常数，计算为 (MM*G*YEAR**2)/(AU**3)。它使现实世界的引力常数适应模拟中使用的归一化单位。

通过定义这些常量，我们具备了准确模拟地球和这个假设的“超级木星”等天体的引力和运动所需的参数。这一步对于确保模拟计算反映现实的天文动力学至关重要。同时，你可以根据需要随意修改这些常量，以在最后观察不同的结果。

创建引力函数

现在，我们将定义一个函数来计算两个天体之间的引力。这个函数对于确定行星的运动至关重要。

## 计算引力的函数
def gravitational_force(m1: float, m2: float, r: np.ndarray) -> np.ndarray:
    """
    计算两个天体之间的引力。
    """
    F_mag = GG * m1 * m2 / (np.linalg.norm(r) + 1e-20)**2
    theta = np.arctan2(np.abs(r[1]), np.abs(r[0]) + 1e-20)
    F = F_mag * np.array([np.cos(theta), np.sin(theta)])
    F *= -np.sign(r)
    return F

在这一步中，我们定义了一个函数 gravitational_force 来计算两个天体之间的引力。这个函数是模拟的关键组成部分，因为它应用牛顿万有引力定律来确定空间中任意两个质量之间施加的力。以下是该函数的详细说明：

1. 函数定义：函数 gravitational_force 接受三个参数：
   • m1 和 m2：两个天体的质量（浮点值）。
   • r：一个 NumPy 数组，表示两个天体之间的位移向量。
2. 计算力的大小（F_mag）：
   • 函数首先使用公式 F_mag = GG * m1 * m2 / (np.linalg.norm(r) + 1e-20)**2 计算引力的大小。
   • GG 是模拟用的引力常数。
   • np.linalg.norm(r) 计算两个天体之间的欧几里得距离。添加一个小数字（1e-20）以防止除以零。
3. 确定方向（theta 和 F）：
   • 使用 np.arctan2 计算角度 theta，以找到极坐标中力向量的角度。同样，在分母中添加一个小数字（1e-20），以避免当 x 方向的位移（r[0]）为零时除以零。
   • 然后使用 F_mag 和角度 theta 计算力向量 F，其分量沿 x 和 y 轴（[np.cos(theta), np.sin(theta)]）。
4. 调整力的方向：
   • 力向量 F 乘以 -np.sign(r)，以确保根据引力定律力始终是吸引的（指向另一个天体）。
5. 返回力向量：
   • 最后，函数返回引力向量 F。

通过实现这个函数，我们学习了如何在计算环境中应用物理定律。该函数对于计算模拟中地球和“超级木星”之间的引力相互作用至关重要，这将影响它们各自的轨迹。

实现四阶龙格 - 库塔求解器

实现四阶龙格 - 库塔（Runge - Kutta）求解器来求解运动的微分方程。此方法广泛用于数值求解常微分方程。

## 四阶龙格 - 库塔求解器
def RK4Solver(t: float, r: np.ndarray, v: np.ndarray, h: float, planet: str, r_other: np.ndarray, v_other: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    用于行星运动的四阶龙格 - 库塔求解器。
    """
    def dr_dt(v: np.ndarray) -> np.ndarray:
        return v

    def dv_dt(r: np.ndarray, planet: str) -> np.ndarray:
        if planet == 'earth':
            return (gravitational_force(ME, MS, r) + gravitational_force(ME, MJ, r - r_other)) / ME
        elif planet == 'jupiter':
            return (gravitational_force(MJ, MS, r) - gravitational_force(MJ, ME, r - r_other)) / MJ

    k11 = dr_dt(v)
    k21 = dv_dt(r, planet)

    k12 = dr_dt(v + 0.5 * h * k21)
    k22 = dv_dt(r + 0.5 * h * k11, planet)

    k13 = dr_dt(v + 0.5 * h * k22)
    k23 = dv_dt(r + 0.5 * h * k12, planet)

    k14 = dr_dt(v + h * k23)
    k24 = dv_dt(r + h * k13, planet)

    y0 = r + h * (k11 + 2 * k12 + 2 * k13 + k14) / 6
    y1 = v + h * (k21 + 2 * k22 + 2 * k23 + k24) / 6

    return y0, y1

在这一步中，我们实现了一个名为 RK4Solver 的函数，它使用四阶龙格 - 库塔（RK4）方法来求解微分方程。此方法对于在我们的引力模型中准确模拟行星运动至关重要。以下是 RK4Solver 函数的概述：

1. 函数定义：
   • RK4Solver 接受几个参数：
     • t：当前时间。
     • r：行星的当前位置向量。
     • v：行星的当前速度向量。
     • h：模拟的时间步长。
     • planet：一个字符串，指示正在模拟的是哪颗行星（地球或木星）。
     • r_other 和 v_other：另一颗行星的位置和速度向量。
2. RK4Solver 中的辅助函数：
   • dr_dt(v: np.ndarray) -> np.ndarray：此函数返回位置的变化率，即速度 v。
   • dv_dt(r: np.ndarray, planet: str) -> np.ndarray：此函数计算由于引力作用行星的加速度。它使用前面定义的 gravitational_force 函数来计算太阳和另一颗行星（地球或木星）施加的力，并返回加速度。
3. 龙格 - 库塔计算：
   • RK4 方法涉及在时间步长内的不同点计算四个“斜率”（k11，k21，k12，k22，k13，k23，k14，k24），然后将它们组合起来以获得下一个时间步长位置和速度的准确估计。
   • 这些斜率是使用当前位置和速度以及位置的导数（dr_dt）和速度的导数（dv_dt）来计算的。
4. 更新位置和速度：
   • 函数根据 RK4 公式使用这些斜率的加权平均值来计算行星的下一个位置（y0）和速度（y1）。
5. 返回结果：
   • 最后，函数返回行星更新后的位置和速度向量（y0，y1）。

通过实现 RK4Solver 函数，我们学习了一种用于数值求解常微分方程的高效且准确的方法。这在涉及行星运动等复杂系统的模拟中尤为重要，其中精度是获得现实结果的关键。RK4 方法在计算效率和准确性之间取得了良好的平衡，使其成为许多科学和工程应用中的热门选择。

设置动画

在运行模拟之前，我们需要设置动画。这一步涉及创建一个绘图并初始化表示行星的线条和标记。

## 设置动画
def setup_animation() -> Tuple[py.Figure, py.Axes, Line2D, Line2D, py.Text]:
    """
    设置动画绘图。
    """
    ## 创建绘图图形和坐标轴
    fig, ax = py.subplots()

    ## 设置坐标轴范围和刻度
    ax.axis('square')
    ax.set_xlim((-7.2, 7.2))
    ax.set_ylim((-7.2, 7.2))
    ax.get_xaxis().set_ticks([])
    ax.get_yaxis().set_ticks([])

    ## 绘制太阳
    ax.plot(0, 0, 'o', markersize=9, markerfacecolor="#FDB813",
            markeredgecolor="#FD7813")

    ## 初始化地球和木星的线条
    line_earth, = ax.plot([], [], 'o-', color='#d2eeff',
                          markevery=10000, markerfacecolor='#0077BE', lw=2)
    line_jupiter, = ax.plot([], [], 'o-', color='#e3dccb', markersize=8,
                            markerfacecolor='#f66338', lw=2, markevery=10000)
    ## 添加一个文本对象
    ttl = ax.text(0.24, 1.05, '', transform=ax.transAxes, va='center')

    ## 返回组件
    return fig, ax, line_earth, line_jupiter, ttl

在这一步中，我们定义了一个名为 setup_animation 的函数，该函数为引力模拟的动画准备绘图。此函数设置了模拟的视觉组件，如绘图区域、行星标记和初始定位。

这个函数对于可视化地球和木星的模拟轨道至关重要。它为模拟的动态部分奠定了基础，在该部分中，这些行星的位置将逐帧更新以创建动画。

创建动画函数

创建一个函数，用于更新动画每一帧中行星的位置。

## 动画函数
def animate(i: int) -> Tuple[Line2D, Line2D, py.Text]:
    """
    行星运动的动画函数。
    """
    earth_trail, jupiter_trail = 40, 200
    tm_yr = 'Elapsed time ='+ str(round(t[i], 1)) +'years'
    ttl.set_text(tm_yr)
    line_earth.set_data(r[i:max(1, i - earth_trail):-1, 0],
                        r[i:max(1, i - earth_trail):-1, 1])
    line_jupiter.set_data(r_jupiter[i:max(
        1, i - jupiter_trail):-1, 0], r_jupiter[i:max(1, i - jupiter_trail):-1, 1])
    return line_earth, line_jupiter, ttl

在这一步中，我们定义了 animate 函数，它是创建行星运动模拟动画的核心。这个函数会被反复调用，以更新绘图上地球和木星的位置，随着时间推移创建运动效果。以下是 animate 函数的工作原理：

1. 地球和木星的轨迹长度：
   • 设置了变量 earth_trail 和 jupiter_trail，它们决定了每个行星显示多少个先前的位置（轨迹）。这会产生视觉上的轨迹效果，展示每个行星走过的路径。
2. 更新时间显示：
   • 计算经过的时间（tm_yr），并将其设置为 ttl 文本对象的文本。这会在动画上显示以年为单位的模拟时间。
3. 更新地球的位置：
   • line_earth.set_data 更新动画中地球的位置。它使用位置数组 r 的一个切片来创建轨迹效果，展示地球最近在其轨道上的位置。
4. 更新木星的位置：
   • 类似地，line_jupiter.set_data 使用 r_jupiter 数组的一个切片来更新木星的位置。木星较长的轨迹长度反映了它较大的轨道。
5. 返回更新后的对象：
   • 函数返回更新后的 Line2D 对象（line_earth，line_jupiter）和文本对象（ttl）。这些是动画每一帧中会变化的元素。

这个 animate 函数对于可视化模拟至关重要，它展示了地球和木星在相互引力影响下如何随时间移动。通过每一帧动态更新绘图，它让模拟变得生动起来。

初始化模拟参数

在开始模拟之前，初始化地球和木星的时间参数、位置和速度。

## 初始化
ti, tf = 0, 120  ## 以年为单位的初始时间和最终时间
N = 100 * tf     ## 每年100个点
t = np.linspace(ti, tf, N)  ## 时间数组
h = t[1] - t[0]  ## 时间步长

## 位置和速度初始化
r = np.zeros([N, 2])         ## 地球的位置
v = np.zeros([N, 2])         ## 地球的速度
r_jupiter = np.zeros([N, 2])  ## 木星的位置
v_jupiter = np.zeros([N, 2])  ## 木星的速度

## 初始条件
r[0] = [1496e8 / AU, 0]
r_jupiter[0] = [5.2, 0]
v[0] = [0, np.sqrt(MS * GG / r[0, 0])]
v_jupiter[0] = [0, 13.06e3 * YEAR / AU]

在这一步中，我们设置了模拟的初始条件和参数，这对于启动地球和木星运动的动画至关重要。以下是所做工作的总结：

1. 时间初始化：
   • ti 和 tf 被设置为模拟的初始时间和最终时间，以年为单位（0到120年）。
   • N 被定义为模拟中的总点数，计算为每年100个点。
   • t 是使用 np.linspace 创建的数组，用于表示从初始时间到最终时间的时间步长。
   • h 是时间步长，计算为 t 数组前两个元素的差值。
2. 位置和速度初始化：
   • 数组 r 和 v 被初始化为存储地球在每个时间步的位置和速度。它们被初始化为零，并且具有适合每个时间步两个坐标（x和y）的形状。
   • 类似地，r_jupiter 和 v_jupiter 被初始化为存储木星的位置和速度。
3. 设置初始条件：
   • 地球的初始位置（r[0]）被设置为其与太阳的距离，通过天文单位（AU）进行归一化。
   • 木星的初始位置（r_jupiter[0]）被设置为距离太阳5.2 AU，反映了它在太阳系中的实际位置。
   • 地球的初始速度（v[0]）是根据太阳施加的引力计算得出的。
   • 木星的初始速度（v_jupiter[0]）被设置为一个特定值，反映了它的实际轨道速度。

这一步是基础，因为它为模拟建立了起点。初始位置和速度对于计算两颗行星在引力作用下的后续运动至关重要。

运行模拟

使用四阶龙格 - 库塔（RK4）求解器在定义的时间范围内执行模拟，并更新行星的位置。

## 运行模拟
for i in trange(N - 1, desc="Generating Animation"):
    r[i + 1], v[i + 1] = RK4Solver(t[i], r[i],
                                   v[i], h, 'earth', r_jupiter[i], v_jupiter[i])
    r_jupiter[i + 1], v_jupiter[i +
                                1] = RK4Solver(t[i], r_jupiter[i], v_jupiter[i], h, 'jupiter', r[i], v[i])

在这一步中，我们执行主循环来运行引力模拟。在这里，使用 RK4Solver 函数在指定的时间段内计算地球和木星的运动。以下是这个过程的简单解释：

1. 运行模拟循环：
   • 设置一个循环，运行 N - 1 次迭代，其中 N 是模拟中的总时间步数。这个循环对于在每个时间步推进模拟至关重要。
   • 使用 tqdm 模块中的 trange 函数来遍历时间步。这提供了一个进度条（描述为“生成动画”）以显示模拟的进度。
2. 更新地球和木星的位置与速度：
   • 在循环的每次迭代中，RK4Solver 函数被调用两次：
     • 一次用于更新地球在下一个时间步的位置（r[i + 1]）和速度（v[i + 1]）。
     • 一次用于更新木星在下一个时间步的位置（r_jupiter[i + 1]）和速度（v_jupiter[i + 1]）。
   • RK4Solver 函数接受行星的当前位置和速度，以及另一颗行星的位置和速度，来计算新的状态。
3. 模拟相互作用：
   • 这个循环本质上在每个小时间步模拟地球和木星之间的引力相互作用（考虑它们彼此之间的影响以及太阳的影响），从而随着时间推移真实地呈现它们的轨道。

在这个循环结束时，我们得到了数组 r、v、r_jupiter 和 v_jupiter，它们在每个时间步都填充了地球和木星的位置和速度，准备好在项目的下一阶段进行动画处理。这一步是应用核心物理和数值方法来模拟天体动力学的地方。

显示动画

最后，运行动画并展示结果。

## 设置动画
fig, ax, line_earth, line_jupiter, ttl = setup_animation()
## 添加刻度和标签
ax.plot([-6,-5],[6.5,6.5],'r-')
ax.text(-4.5,6.3,r'1 AU = $1.496 \times 10^8$ km')

ax.plot(-6,-6.2,'o', color = '#d2eeff', markerfacecolor = '#0077BE')
ax.text(-5.5,-6.4,'地球')

ax.plot(-3.3,-6.2,'o', color = '#e3dccb',markersize = 8, markerfacecolor = '#f66338')
ax.text(-2.9,-6.4,'超级木星（质量为500倍）')

ax.plot(5,-6.2,'o', markersize = 9, markerfacecolor = "#FDB813",markeredgecolor ="#FD7813")
ax.text(5.5,-6.4,'太阳')

## 创建动画
anim = animation.FuncAnimation(
    fig, animate, frames=4000, interval=1, blit=False)

## 显示动画
plt.show()

在项目的最后一步，我们展示引力模拟的动画，该动画展示了地球和假设的“超级木星”的运动。这一步涉及可视化模拟期间计算出的行星位置。

现在我们已经完成了所有步骤，可以在桌面环境中使用以下命令运行代码：

cd ~/project
python simulation.py

在这一步，我们可以直观地欣赏模拟结果，观察巨大的“超级木星”如何影响地球的轨道。

总结

在这个项目中，我们有效地使用Python开发了一个以地球和木星为特色的引力模拟。我们采用了四阶龙格 - 库塔（RK4）方法进行精确的数值积分，并使用Matplotlib进行清晰且引人入胜的可视化。这个项目不仅提供了物理和编码原理的实际应用，还鼓励学生通过修改代码来模拟不同质量的天体，从而探索一系列引力相互作用和动力学，进行实验。