Scipy odeint积分在最大值处终止，无法生成完整宇宙尺度因子演化图

看起来你遇到的是闭宇宙（k=+1）积分终止时机和数值精度的问题，结合Friedmann方程的特性和SciPy积分工具的特点，我整理了几个关键修复点和完整的解决方案：

核心问题分析 #

闭宇宙的尺度因子会达到峰值后收缩至0，但scipy.integrate.odeint默认会积分到你指定的固定时间终点，不会自动在a(t)=0时停止。另外，当a趋近于0时，Friedmann方程的导数会出现奇点，若积分精度不足或未处理根号内的负数情况，会导致曲线偏离预期甚至出现数值异常。

解决方案步骤 #

1. 改用更适合的积分工具：solve_ivp #

odeint是较旧的积分接口，推荐使用scipy.integrate.solve_ivp，它支持事件触发终止，可以精准在a(t)=0时停止积分，完美适配闭宇宙的演化终点。

2. 修正Friedmann方程的实现 #

确保方程形式正确（以物质主导宇宙为例），并处理根号内可能出现的负数情况：
$$\dot{a} = \sqrt{\frac{8\pi G \rho_0}{3a} - k c^2}$$

3. 提升积分精度 #

通过设置rtol和atol参数降低数值误差，避免曲线偏离理论解。

完整修复代码 #

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 宇宙学参数设置（统一为国际单位）
H0 = 67.8  # 哈勃常数，单位：km/s/Mpc
H0_si = H0 * 1000 / 3.086e22  # 转换为s⁻¹
G = 6.67430e-11  # 引力常数，单位：m³/kg/s²
c = 3e8  # 光速，单位：m/s
Omega0 = 1.2  # 物质密度参数，闭宇宙需Ω₀>1
rho_c = 3 * H0_si**2 / (8 * np.pi * G)  # 临界密度
rho0 = Omega0 * rho_c  # 当前物质密度

# 单位转换工具：秒 ↔ 千兆年（Gyr）
sec_to_gyr = lambda t: t / (3.154e16 * 1e9)
gyr_to_sec = lambda t: t * 3.154e16 * 1e9

def friedmann(t, a, k):
    # solve_ivp要求函数格式为f(t, y)，与odeint的f(y, t)相反
    term = (8 * np.pi * G * rho0) / (3 * a) - k * c**2
    # 避免根号内出现负数，当term≤0时返回0（对应收缩到终点）
    return np.sqrt(term) if term > 0 else 0.0

def collapse_event(t, a, k):
    # 定义终止事件：当尺度因子a趋近于0时停止积分
    return a[0]
collapse_event.terminal = True  # 触发事件时终止积分
collapse_event.direction = -1  # 仅当a从正变负时触发（避免初始小a误触发）

# 初始条件：a(t=0)取极小值（避免a=0的奇点）
a0 = [1e-10]
# 积分时间范围：0到50Gyr（足够覆盖闭宇宙的寿命）
t_span = [0, gyr_to_sec(50)]

# 不同曲率的宇宙参数
curvatures = [(-1, 'red', '开放宇宙 (k=-1)'),
              (0, 'green', '平坦宇宙 (k=0)'),
              (1, 'lightblue', '闭合宇宙 (k=1)')]

plt.figure(figsize=(10, 6))
for k, color, label in curvatures:
    # 执行积分，加入事件触发和高精度设置
    sol = solve_ivp(friedmann, t_span, a0, args=(k,),
                    events=collapse_event, rtol=1e-10, atol=1e-12)
    # 绘制曲线，转换时间单位为Gyr
    plt.plot(sec_to_gyr(sol.t), sol.y[0], color=color, label=label)

plt.xlabel('时间 (Gyr)')
plt.ylabel('尺度因子 a(t)')
plt.title('不同曲率下宇宙尺度因子随时间演化')
plt.legend()
plt.grid(True)
plt.xlim(0, 50)
plt.ylim(0, None)
plt.show()

额外验证建议 #

1. 检查闭宇宙的理论寿命：用解析公式计算闭宇宙的最大寿命，确认是否与你预期的43Gyr一致：
   $$t_{max} = \frac{2\pi}{3H_0} \sqrt{\frac{\Omega_0}{\Omega_0 - 1}}$$
   代入参数计算后，若结果接近43Gyr，说明参数设置正确。
2. 调整初始条件：若初始a0过小，可能导致积分初期数值不稳定，可尝试微调为1e-8或1e-9。

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