含悬殊时间尺度的衰变链模拟方法技术问询

针对半衰期跨度极大的衰变链，以下几种方法可以有效解决统一时间步长效率低下的问题：

1. 事件驱动蒙特卡洛模拟 #

这种方法不依赖固定时间步，而是直接模拟每个衰变事件的发生时刻，完全适配不同时间尺度的核素。核心逻辑是：

• 对每个核素，基于其衰变常数λ（λ=ln2/半衰期），计算单个核的下一次衰变时间（服从指数分布：t_next = -ln(random.random())/λ）
• 用优先队列（比如Python的heapq）维护所有待发生的衰变事件，每次取出最早发生的事件
• 执行衰变：减少母核数量，增加子核数量（如果子核是放射性的，还要为新产生的子核计算下一次衰变时间并加入队列）
• 重复上述步骤，直到模拟时间达到目标时长

示例代码片段：

import heapq
import random
import math

# 定义核素参数：(名称, 半衰期s, 子核名称)
nuclides = {
    "Po-214": (0.0001643, "Pb-210"),
    "Ra-226": (307106512477175.9, "Rn-222"),
    "Pb-210": (223996800, "Bi-210"),
    # 补充其他核素参数...
}

# 初始化：核素数量 + 事件队列（(衰变时间, 核素名称)）
counts = {"Ra-226": 10000, "Po-214": 5000}
event_queue = []

# 初始化事件队列
for nuclide, cnt in counts.items():
    half_life, _ = nuclides[nuclide]
    lam = math.log(2) / half_life
    for _ in range(cnt):
        t_decay = -math.log(random.random()) / lam
        heapq.heappush(event_queue, (t_decay, nuclide))

total_time = 1e6  # 模拟总时长，单位秒
current_time = 0

while event_queue and current_time < total_time:
    t_decay, nuclide = heapq.heappop(event_queue)
    if t_decay > total_time:
        break
    current_time = t_decay
    
    # 母核衰变，数量减1
    counts[nuclide] -= 1
    if counts[nuclide] == 0:
        del counts[nuclide]
    
    # 生成子核
    half_life, daughter = nuclides[nuclide]
    counts[daughter] = counts.get(daughter, 0) + 1
    
    # 为子核添加衰变事件（如果子核有半衰期）
    if daughter in nuclides:
        d_half_life, _ = nuclides[daughter]
        d_lam = math.log(2) / d_half_life
        d_t_decay = current_time - math.log(random.random()) / d_lam
        heapq.heappush(event_queue, (d_t_decay, daughter))

print("模拟结束时核素数量：", counts)

2. 分阶段多时间步长模拟 #

将模拟过程分为两个阶段，适配不同时间尺度：

• 快速阶段：用极小时间步模拟短半衰期核素，直到它们的数量变化率低于设定阈值（比如<1e-6），此时短寿命核素已接近稳态
• 稳态阶段：切换为极大时间步，只模拟长半衰期核素的衰变，短寿命核素的数量用稳态公式（ secular equilibrium：N_daughter = (λ_parent / λ_daughter) * N_parent）计算，无需逐时间步更新

这种方法兼顾效率和精度，适合短寿命核素能快速达到平衡的场景。

3. 解析近似与数值模拟结合 #

对于短半衰期的核素，直接用解析解计算其稳态浓度，避免逐事件或逐时间步模拟：

• 当衰变链中存在“快衰变-慢衰变”的层级时，快衰变核素会很快进入 secular equilibrium，此时子核的衰变率等于母核的衰变率
• 模拟时只需跟踪长半衰期核素的数量变化，短寿命核素的数量直接通过平衡公式推导得出，大幅减少计算量

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