构建需反复求解常数的无限递归分段函数技术问题
问题背景
需要构建可运行至t极大值（如36,000,000秒）的分段连续函数，每段结束时需重新计算常数C以保证函数连续性。目前可处理前1200秒的情况，但在多次递归计算常数C时遇到困难，最终需实现Mn55(t)与Mn56(t)的绘图。
给定参数与初始条件
# 初始条件 Mn55_0 = 0 Mn56_0 = 0 # 参数 p = 1E13 # 对应微分方程中的phi u = 1.10872E-05 lambda_val = 0.26659 # 避免与Python关键字lambda冲突 Q = 1E13 # 对应解析解中的q
控制微分方程
分段微分方程如下：
Mn55的微分方程
当满足
t%1200 == 1时：
$$\frac{dMn55}{dt} = -u \cdot p \cdot Mn55 + Q$$
否则：
$$\frac{dMn55}{dt} = Q$$
Mn56的微分方程
当满足
t%1200 == 1时：
$$\frac{dMn56}{dt} = u \cdot p \cdot Mn55(t) - \lambda \cdot Mn56(t)$$
否则：
$$\frac{dMn56}{dt} = - \lambda \cdot Mn56(t)$$
解析解与连续性实现思路
核心是分区间迭代计算，每个1200秒为一个周期，每个周期内分两个子区间处理：
- 对每个子区间，以上一个区间结束时的Mn55、Mn56值作为初始条件，求解当前区间的常数C
- 计算当前区间内所有时间点的函数值
- 迭代处理所有周期直到达到目标时间
Mn55的解析解
当
t%1200 == 1时：
$$Mn55(t) = -C_1 \cdot e^{-p u t} + \frac{Q}{p u}$$
利用区间起始点的Mn55值求解$C_1$：$C_1 = \frac{Q}{p u} - Mn55(t_{\text{start}}) \cdot e^{p u t_{\text{start}}}$否则：
$$Mn55(t) = Q \cdot t + C_2$$
利用区间起始点的Mn55值求解$C_2$：$C_2 = Mn55(t_{\text{start}}) - Q \cdot t_{\text{start}}$
Mn56的解析解
当
t%1200 == 1时：
$$Mn56(t) = C_3 \cdot e^{-\lambda t} + \frac{p u}{\lambda} \cdot Mn55(t) \cdot t - \frac{p u}{\lambda^2} \cdot Mn55(t)$$
利用区间起始点的Mn56值求解$C_3$：$C_3 = \left[ Mn56(t_{\text{start}}) - \frac{p u}{\lambda} \cdot Mn55(t_{\text{start}}) \cdot t_{\text{start}} + \frac{p u}{\lambda^2} \cdot Mn55(t_{\text{start}}) \right] \cdot e^{\lambda t_{\text{start}}}$否则：
$$Mn56(t) = C_4 \cdot e^{-\lambda t}$$
利用区间起始点的Mn56值求解$C_4$：$C_4 = Mn56(t_{\text{start}}) \cdot e^{\lambda t_{\text{start}}}$
Python实现代码
import numpy as np import matplotlib.pyplot as plt # 初始化参数 Mn55_prev = 0.0 Mn56_prev = 0.0 p = 1E13 u = 1.10872E-05 lambda_val = 0.26659 Q = 1E13 max_t = 36000000 # 目标极大时间 period = 1200 # 周期长度 step = 1 # 时间步长 # 存储结果的列表 t_list = [] mn55_list = [] mn56_list = [] current_t = 0 while current_t <= max_t: # 阶段1：当前时间处于周期内的激活段（对应t%1200 ==1 逻辑） phase1_end = ((current_t // period) + 1) * period if current_t % period !=0 else current_t +1 phase1_end = min(phase1_end, max_t) # 处理阶段1 t_phase1 = np.arange(current_t, phase1_end, step) if len(t_phase1) >0: # 计算C1 C1 = (Q/(p*u)) - Mn55_prev * np.exp(p*u*current_t) # 计算Mn55 mn55_phase1 = -C1 * np.exp(-p*u*t_phase1) + Q/(p*u) # 计算C3 term = (p*u/lambda_val)*Mn55_prev*current_t - (p*u/(lambda_val**2))*Mn55_prev C3 = (Mn56_prev - term) * np.exp(lambda_val*current_t) # 计算Mn56 mn56_phase1 = C3 * np.exp(-lambda_val*t_phase1) + (p*u/lambda_val)*mn55_phase1*t_phase1 - (p*u/(lambda_val**2))*mn55_phase1 # 添加到结果列表 t_list.extend(t_phase1) mn55_list.extend(mn55_phase1) mn56_list.extend(mn56_phase1) # 更新前值 Mn55_prev = mn55_phase1[-1] Mn56_prev = mn56_phase1[-1] current_t = phase1_end # 阶段2：当前时间处于周期内的非激活段 phase2_end = ((current_t // period) +1)*period phase2_end = min(phase2_end, max_t) if current_t >= phase2_end: break t_phase2 = np.arange(current_t, phase2_end, step) if len(t_phase2) >0: # 计算C2 C2 = Mn55_prev - Q * current_t # 计算Mn55 mn55_phase2 = Q * t_phase2 + C2 # 计算C4 C4 = Mn56_prev * np.exp(lambda_val*current_t) # 计算Mn56 mn56_phase2 = C4 * np.exp(-lambda_val*t_phase2) # 添加到结果列表 t_list.extend(t_phase2) mn55_list.extend(mn55_phase2) mn56_list.extend(mn56_phase2) # 更新前值 Mn55_prev = mn55_phase2[-1] Mn56_prev = mn56_phase2[-1] current_t = phase2_end # 绘图 plt.figure(figsize=(12,8)) plt.subplot(2,1,1) plt.plot(t_list, mn55_list, label='Mn55(t)') plt.xlabel('时间 t (秒)') plt.ylabel('Mn55 浓度') plt.title('Mn55 随时间演化') plt.legend() plt.subplot(2,1,2) plt.plot(t_list, mn56_list, label='Mn56(t)', color='orange') plt.xlabel('时间 t (秒)') plt.ylabel('Mn56 浓度') plt.title('Mn56 随时间演化') plt.legend() plt.tight_layout() plt.show()
代码说明
- 采用迭代方式处理每个1200秒周期，每个周期内分两个阶段计算
- 每个阶段开始时，用上一阶段结束的Mn55、Mn56值求解当前阶段的常数C，保证函数连续性
- 使用numpy批量计算区间内的函数值，提升计算效率，适合处理极大时间值
- 最后用matplotlib绘制两个核素的演化曲线
内容的提问来源于stack exchange，提问作者Seth Dayawansa
