基于重力原子谱的同位素量化代码报错：零除与ODEint非法输入问题

问题背景 #

正在编写代码定义"null"和"void"取值，量化同位素ley线、光子反应物性质，将时间作为同位素密度衰减参数，引入"重力原子谱（gravatomic spectrum）"的基本重力性质作为常量。代码运行时出现以下错误：

RuntimeWarning: divide by zero encountered in scalar divide sigma = F / dphi_dt
ODEintWarning: Illegal input detected (internal error). Run with full_output = 1 to get quantitative information. sol = odeint(system, y0, t, args=(phi0, lambda_val, epsilon))

原代码：

import numpy as np
from scipy.integrate import odeint

# Define the function representing the differential equations
def system(y, t, phi0, lambda_val, epsilon):
    rho, phi, r = y
    
    # Isotope density decay equation
    drho_dt = -lambda_val * rho
    
    # Change in energy of photonic reactants equation
    dphi_dt = phi - phi0
    
    # Gravitational force equation
    F = G * epsilon / r**2
    
    # Properties of the gravatomic spectrum equation
    sigma = F / dphi_dt
    
    return [drho_dt, dphi_dt, sigma]

# Define initial conditions
rho0 = 1.0  # Initial isotope density
phi0 = 0.0  # Baseline energy of photonic reactants
r0 = 1.0    # Initial distance

# Set parameters
lambda_val = 0.1  # Decay constant
G = 6.67430e-11   # Gravitational constant
epsilon = 1.0     # Efficiency factor of gravitational interaction

# Create time points for integration
t = np.linspace(0, 10, 100)

# Initial values
y0 = [rho0, phi0, r0]

# Solve the system of ODEs
sol = odeint(system, y0, t, args=(phi0, lambda_val, epsilon))

# Extract the solutions
rho_sol = sol[:, 0]
phi_sol = sol[:, 1]
sigma_sol = sol[:, 2]

错误原因分析 #

1. 除零错误：初始条件中phi0 = 0.0且初始phi取值等于phi0，导致dphi_dt = phi - phi0 = 0，计算sigma = F / dphi_dt时触发除以零的运算。
2. ODE系统定义错误：odeint要求状态向量的每个元素必须对应其时间导数。当前状态向量是[rho, phi, r]，但返回值的第三个元素是sigma（重力原子谱性质），而非r的时间导数dr_dt。这违反了ODE积分的核心要求，导致数值积分过程中出现非法输入。

解决步骤 #

1. 修正ODE系统的返回值 #

将sigma从ODE系统的返回列表中移除，改为在积分完成后基于结果计算。同时补充r的时间导数方程（需根据实际物理模型定义，示例中使用重力作用下的运动方程作为参考，可按需替换）。

2. 规避除零情况 #

可通过两种方式处理：

• 调整初始条件，让初始phi与phi0存在微小差值，比如phi0 = 0.0，初始phi设为1e-8；
• 若物理模型允许，修改dphi_dt的方程，避免初始时为零，例如加入阻尼项：dphi_dt = (phi - phi0) - k*phi（k为阻尼系数）。

3. 修正后的代码 #

import numpy as np
from scipy.integrate import odeint

# 定义微分方程组，状态向量为[rho, phi, r]，返回对应时间导数
def system(y, t, phi0, lambda_val, epsilon, G):
    rho, phi, r = y
    
    # 同位素密度衰减方程
    drho_dt = -lambda_val * rho
    
    # 光子反应物能量变化方程
    dphi_dt = phi - phi0
    # 可选：加入阻尼项防止长期发散，比如dphi_dt = (phi - phi0) - 0.1*phi
    
    # r的时间导数：示例为重力加速度对应的运动（假设质量为1）
    # 若有其他物理规则，替换此处的dr_dt计算逻辑
    F = G * epsilon / r**2
    dr_dt = F  # 假设加速度等于力（质量=1）
    
    return [drho_dt, dphi_dt, dr_dt]

# 初始条件
rho0 = 1.0
phi0 = 0.0
# 初始phi设为微小非零值，避免除零
initial_phi = 1e-8
r0 = 1.0

# 参数
lambda_val = 0.1
G = 6.67430e-11
epsilon = 1.0

# 时间点
t = np.linspace(0, 10, 100)

# 初始状态向量
y0 = [rho0, initial_phi, r0]

# 求解ODE
sol = odeint(system, y0, t, args=(phi0, lambda_val, epsilon, G))

# 提取结果
rho_sol = sol[:, 0]
phi_sol = sol[:, 1]
r_sol = sol[:, 2]

# 积分完成后计算重力原子谱性质sigma，加入判断避免除零
sigma_sol = []
for phi, r in zip(phi_sol, r_sol):
    dphi_dt = phi - phi0
    if abs(dphi_dt) < 1e-10:
        # 处理接近零的情况，设为默认值或按需调整
        sigma_sol.append(0.0)
    else:
        F = G * epsilon / r**2
        sigma_sol.append(F / dphi_dt)
sigma_sol = np.array(sigma_sol)

说明 #

• 若r的时间导数有特定物理模型，请替换示例中的dr_dt计算逻辑；
• 除零处理中的阈值（如1e-10）可根据精度需求调整；
• 若dphi_dt的物理模型需要保持原形式，必须确保初始时刻phi与phi0不相等，或在计算sigma时加入分支判断处理零值情况。

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