Python数值求解耦合二阶ODE的边界条件适配问题

问题背景 #

需要数值求解以下耦合二阶常微分方程组：
$$
\begin{cases}
g'' + \frac{2}{r}g' - \frac{2}{r^2}g - \frac{3e}{r}g^2 - e^2g3 - \frac{e}{2r}f^2 - \frac{e^2}{4}f2g = 0 \
f'' + \frac{2}{r}f' - \frac{2}{r^2}f - \frac{2e}{r}fg - \frac{e^2}{2}fg2 + \mu^2f - \lambda f^3 = 0
\end{cases}
$$

通过变量替换 $g'=v$、$f'=u$，转化为4个耦合一阶方程组：
$$
\begin{cases}
\frac{dg}{dr} = v \
\frac{dv}{dr} = -\frac{2}{r}v + \frac{2}{r^2}g + \frac{3e}{r}g^2 + e^2g3 + \frac{e}{2r}f^2 + \frac{e^2}{4}f2g \
\frac{df}{dr} = u \
\frac{du}{dr} = -\frac{2}{r}u + \frac{2}{r^2}f + \frac{2e}{r}fg + \frac{e^2}{2}fg2 - \mu^2f + \lambda f^3
\end{cases}
$$

物理边界条件为：

• $g(0)=0,\ f(0)=0$
• $g(\infty)=0,\ f(\infty)=\frac{\mu}{\sqrt{\lambda}}$

当前使用solve_ivp（初值求解器）猜测初始导数的方式无法满足边界条件，结果与文献预期不符。

解决方案 #

核心思路 #

这是典型的边值问题（BVP），需要使用专门的边值求解器scipy.integrate.solve_bvp，它会迭代调整初始参数，让解同时满足两端的边界条件。

实现代码 #

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

# 常数设置
e = 1
mu = 1
lambd = 1
f_inf = mu / np.sqrt(lambd)  # f在无穷远处的边界值

# 定义一阶微分方程组
def dSdr(r, S):
    g, v, f, u = S
    dgdr = v
    dvdr = -2 / r * v + 2 / r**2 * g + 3 * e / r * g**2 + \
           e**2 * g**3 + e / (2 * r) * f**2 + e**2 / 4 * f**2 * g
    dfdr = u
    dudr = -2 / r * u + 2 / r**2 * f + 2 * e / r * f * g + \
           e**2 / 2 * f * g**2 - mu**2 * f + lambd * f**3
    return np.vstack([dgdr, dvdr, dfdr, dudr])

# 定义边界条件残差：残差为0时满足边界条件
def bc(ya, yb):
    # ya是r=start处的解，yb是r=end处的解
    return np.array([
        ya[0],          # g(0) = 0
        ya[2],          # f(0) = 0
        ya[1],          # g'(0) = 0（光滑性要求）
        ya[3],          # f'(0) = 0（光滑性要求）
        yb[0],          # g(∞) = 0
        yb[2] - f_inf   # f(∞) = mu/sqrt(lambda)
    ])

# 设置r的区间：避免r=0，取极小起始值，终点取足够大以趋近稳态
r = np.linspace(1e-3, 50, 200)

# 构造初始猜测解，贴合预期趋势
g_guess = np.zeros_like(r)
v_guess = np.zeros_like(r)
f_guess = f_inf * (1 - np.exp(-r/2))  # 从0趋近f_inf的曲线
u_guess = f_inf/2 * np.exp(-r/2)      # f_guess的导数
S_guess = np.vstack([g_guess, v_guess, f_guess, u_guess])

# 求解边值问题
sol = solve_bvp(dSdr, bc, r, S_guess)

# 检查求解状态
if sol.success:
    print("求解成功")
else:
    print("求解失败：", sol.message)

# 绘图
plt.figure(figsize=(10,6))
plt.plot(sol.x, sol.y[0], label='g(r)')
plt.plot(sol.x, sol.y[2], label='f(r)')
plt.xlabel('r')
plt.ylabel('函数值')
plt.legend()
plt.title('耦合方程组的数值解')
plt.grid(True)
plt.show()

关键说明 #

1. 边界条件处理：通过bc函数定义残差，当残差全为0时，解满足所有边界条件。其中r=0处的导数为0是由物理光滑性推导的（原点处无奇异，函数对称）。
2. 初始猜测：初始猜测的曲线越贴合预期解，求解器越容易收敛。这里构造的f(r)是指数趋近稳态的曲线，g(r)初始设为0，求解器会自动调整出正确的先增后减曲线。
3. 区间选择：r起始值取1e-3避免分母为0的奇异，终点取50足够让f(r)趋近稳态值。

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