使用scipy.solve_bvp求解微分方程组时遇溢出问题求助

需要求解的微分方程组：

d(phi)/dz = Rcl*i
d(i)/dz = i0*exp(phi*2.3/b)

边界条件：$z=0$时$i=il=0$；$z=1$时$i=ir=1$

原solve_bvp代码的核心问题 #

• 初始猜测过于粗糙：全零初始猜测无法引导非线性算法（指数项存在）找到正确解的趋势。从你的打靶法结果可知，$\phi$初始值约为0.31，$i$从0单调增长到1，零猜测与真实解偏差过大，导致算法收敛失败。
• 参数命名冲突：代码中将区间上限命名为b，覆盖了之前定义的参数b = 50e-3，直接导致后续指数项计算错误。
• 未校验求解状态：未判断solution.success属性就直接调用解，若求解失败会引发后续绘图错误。

改进后的solve_bvp实现 #

1. 修正参数命名冲突 #

将区间上限变量名改为z_end，避免覆盖参数b。

2. 构造合理初始猜测 #

基于打靶法的已知结果，构造接近真实解的初始猜测：

• $i$的初始猜测设为从0到1的线性函数，符合边界条件趋势
• $\phi$的初始猜测基于打靶法得到的初始值，结合$\frac{d\phi}{dz}=Rcl*i$的积分关系近似生成

3. 可选调整求解参数 #

对于非线性较强的问题，增大max_nodes允许算法使用更多节点，提升收敛概率。

修正后的完整代码 #

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

# 定义参数
b = 50e-3
Rcl = 0.4
i0 = 1e-7

il = 0
ir = 1

# 定义ODE系统
def fun(z, y):
    phi = y[0]
    i = y[1]
    didz = i0 * np.exp(phi * 2.3 / b)
    dphidz = Rcl * i
    return np.vstack([dphidz, didz])  # 确保输出为二维数组，符合solve_bvp要求

# 定义边界条件
def boundary_conditions(ya, yb):
    return np.array([ya[1] - il, yb[1] - ir])

# 数值参数
a = 0
z_end = 1  # 修正命名冲突
num_points = 100

# 构造初始猜测
x_mesh = np.linspace(a, z_end, num_points)
y_initial_guess = np.zeros((2, num_points))

# i的初始猜测：从il到ir线性增长
y_initial_guess[1] = np.linspace(il, ir, num_points)

# phi的初始猜测：基于打靶法的初始值，结合i的积分趋势
phi_start = 0.31  # 来自打靶法结果
dz = x_mesh[1] - x_mesh[0]
phi_guess = phi_start + Rcl * np.cumsum(y_initial_guess[1][:-1]) * dz
y_initial_guess[0] = np.concatenate([[phi_start], phi_guess])

# 求解边值问题，增大最大节点数提升收敛性
solution = solve_bvp(fun, boundary_conditions, x_mesh, y_initial_guess, max_nodes=500)

# 校验求解状态并绘图
if solution.success:
    print("求解成功")
    x_fine = np.linspace(a, z_end, 1000)
    y_fine = solution.sol(x_fine)
    
    fig, ax = plt.subplots()
    ax.plot(x_fine, y_fine[0], label='Phi')
    ax.legend()
    
    fig, ax = plt.subplots()
    ax.plot(x_fine, y_fine[1], label='i')
    ax.legend()
    
    plt.show()
else:
    print("求解失败，原因：", solution.message)

复杂方程组扩展建议 #

对于后续更复杂的微分方程组，建议：

• 先用简化模型或打靶法得到近似解，作为solve_bvp的初始猜测，这是提升收敛概率的核心
• 利用solve_bvp的callback函数监控求解过程，动态调整参数
• 若问题刚性极强，可保留自定义打靶法框架，改用scipy.integrate.solve_ivp替代手动欧拉法，提升初值求解精度

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