耦合常微分方程下打靶法实现：石英光纤激光功率仿真求解

你目前的实现有三个核心错误导致长距离求解失败：

1. 把两点边值问题当成了初值问题处理：双向传输的光纤功率耦合方程属于典型的边值问题，两端都有反射约束，你直接把z=0处所有反向功率设为0完全不符合物理边界，90km长度下增益、损耗、散射的累积效应会直接导致数值溢出，求解器自然报错。
2. 用了不适配刚性方程的求解器：你的耦合方程里增益项、损耗项、瑞利散射项的特征尺度差了好几个数量级，属于刚性方程，solve_ivp默认的显式RK45算法处理这类问题时步长会被压到极小，要么跑不动要么直接发散。
3. 完全缺失打靶法的核心求根环节：打靶法的本质是迭代调整未知的初始边界值，让积分到另一端的结果满足边界条件，你现在的代码就是直接给了个初始值硬跑，根本不是打靶法。

1. 先明确两端边界约束 #

根据你代码里的反射率参数，先把两端的边界条件理清楚：

• z=0（输入端）：正向泵浦功率固定为输入功率Pin，所有反向波的功率满足P_b(0) = R_L * P_f(0)，你注释里的输入端FBG反射参数可以直接用在这里。
• z=L（输出端）：所有反向波的功率满足P_b(L) = R_R * P_f(L)，你代码里定义的R=0.6可以作为输出端反射率使用。

2. 打靶法实现逻辑 #

按以下步骤写迭代逻辑即可：

• 定义待优化的未知量：z=0处的3个斯托克斯正向功率、4个反向功率，共7个变量（正向泵浦是固定的Pin，不需要优化）。
• 定义残差函数：输入未知的初始值，拼接成完整的8维初始功率向量，用刚性求解器积分到z=L，计算z=L处的边界残差（即P_b(L) - R*P_f(L)的误差和）。
• 调用scipy.optimize.root或者fsolve迭代优化未知初始值，直到残差收敛到阈值（比如1e-6）以下。

3. 长距离求解优化技巧 #

• 不要设置points = fiberlength*100这么高的采样点，90km的场景用1000个采样点完全足够，过高的采样点只会大幅拖慢求解速度。
• 积分时必须使用刚性求解器，比如solve_ivp的Radau或者BDF方法，示例：

sol = solve_ivp(odes2, (0, fiberlength), P0, method='Radau', rtol=1e-6, atol=1e-9)

• 初始猜测值不要全设为0，可以先用短距离（比如500m）的求解结果作为长距离的初始猜测，或者按功率量级设为1e-6~1e-3W的小值，避免迭代发散。
• 如果用solve_bvp求解，需要单独写边界条件函数，示例模板如下：

def boundary_conditions(Pa, Pb):
    # Pa为z=0处的功率向量，Pb为z=L处的功率向量
    return [
        Pa[0] - Pin,  # 输入端正向泵浦固定为Pin
        Pa[1] - R_L * Pa[0],  # 输入端泵浦反向满足反射约束
        Pa[3] - R_L * Pa[2],  # 输入端1阶斯托克斯反向满足反射约束
        Pa[5] - R_L * Pa[4],  # 输入端2阶斯托克斯反向满足反射约束
        Pa[7] - R_L * Pa[6],  # 输入端3阶斯托克斯反向满足反射约束
        Pb[1] - R_R * Pb[0],  # 输出端泵浦反向满足反射约束
        Pb[3] - R_R * Pb[2],  # 输出端1阶斯托克斯反向满足反射约束
        Pb[5] - R_R * Pb[4],  # 输出端2阶斯托克斯反向满足反射约束
        Pb[7] - R_R * Pb[6],  # 输出端3阶斯托克斯反向满足反射约束
    ]

# 初始猜测网格，先设稀疏网格降低计算量
x = np.linspace(0, fiberlength, 100)
y_guess = np.zeros((8, x.size))
y_guess[0] = np.linspace(Pin, Pin*np.exp(-alpha*fiberlength), x.size) # 泵浦正向初始猜测为衰减曲线
sol_bvp = solve_bvp(odes2, boundary_conditions, x, y_guess, tol=1e-5, max_nodes=10000)

你可以再核对一遍ODE函数里反向波的导数符号：正向波沿z正方向传输，损耗项为负、增益项为正；反向波沿z负方向传输，导数符号和正向波相反，这里写错也很容易导致数值发散。

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