Python数值求解非线性耦合微分方程组遇阻求助

核心问题分析 #

• 平凡解陷阱：你设置的初始条件[0,0,0,0]是方程组的平凡解（所有变量为0时方程成立），因此无论用什么数值方法，解都会保持为0。
• 边值问题误用初值求解器：你的边界条件是无穷远的（$f(\infty)=1, g(\infty)=0$），这属于边值问题（BVP），而solve_ivp是用于解**初值问题（IVP）**的工具，直接用初值方法无法满足远端边界条件。
• 一阶方程组推导错误：原方程移项后的符号和你代码中的dvdr、dudr存在多处不一致，导致方程建模本身错误。

修正步骤 #

1. 修正一阶微分方程组的推导 #

对比原方程，正确的一阶方程组应为：
$$
\begin{align*}
\frac{dg}{dr} &= v \
\frac{dv}{dr} &= -\frac{2}{r}v + \frac{2}{r^2}g + \frac{3}{r}e g^2 + e^2 g^3 + \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^2 f + \lambda f^3
\end{align*}
$$
（注：原方程中$\nu^2$对应代码中的mu）

2. 使用边值问题求解器solve_bvp #

Scipy提供了solve_bvp专门处理边值问题，需要定义：

• 微分方程组的右端函数
• 边界条件的残差函数（即边界处期望的差值）
• 非零的初始猜测解（避免平凡解）

3. 处理无穷远边界条件 #

将无穷远截断到足够大的r_max（比如100），假设在r_max处：

• $f(r_{max}) = 1$，$f'(r_{max}) = 0$（趋近于常数时导数为0）
• $g(r_{max}) = 0$，$g'(r_{max}) = 0$（趋近于0时导数为0）

近端（$r=r_0=0.1$）可根据正则性要求设置边界条件（比如导数为0，避免发散），具体可根据问题物理背景调整。

修正后的代码 #

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

# Constants
e = 1
mu = 1
lambd = 1

# 定义微分方程组：dy/dr = f(r, y)，y = [g, g', f, f']
def dydr(r, y):
    g, dgdr, f, dfdr = y
    # 计算g''
    d2gdr2 = (-2/r)*dgdr + (2/r**2)*g + (3*e/r)*g**2 + e**2*g**3 + (e/(2*r))*f**2 + (e**2/4)*f**2*g
    # 计算f''
    d2fdr2 = (-2/r)*dfdr + (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, d2gdr2, dfdr, d2fdr2])

# 定义边界条件残差：边界处的期望差值为0
def bc(ya, yb):
    # ya是r=r0处的解，yb是r=rmax处的解
    # 远端边界条件：f(rmax)=1, f'(rmax)=0; g(rmax)=0, g'(rmax)=0
    # 近端边界条件：假设g'(r0)=0、f'(r0)=0（球对称中心正则性，可根据问题调整）
    return np.array([
        ya[1],          # g'(r0) = 0
        ya[3],          # f'(r0) = 0
        yb[0] - 0,      # g(rmax) = 0
        yb[1],          # g'(rmax) = 0
        yb[2] - 1,      # f(rmax) = 1
        yb[3]           # f'(rmax) = 0
    ])

# 设置r的区间
r0 = 0.1
rmax = 100
r = np.linspace(r0, rmax, 200)

# 初始猜测解：避免全0，构造接近预期的形状
guess_g = np.exp(-(r-5)**2/10)  # 高斯峰值模拟g的小峰值
guess_dgdr = - (r-5)/5 * np.exp(-(r-5)**2/10)
guess_f = 1 - np.exp(-r/5)      # 指数上升趋近于1的f
guess_dfdr = (1/5)*np.exp(-r/5)
y_guess = np.vstack([guess_g, guess_dgdr, guess_f, guess_dfdr])

# 求解边值问题
sol = solve_bvp(dydr, bc, r, y_guess, tol=1e-6, max_nodes=10000)

# 检查求解结果
if sol.success:
    print("求解成功")
    print(f"残差最大值: {np.max(sol.residual)}")
else:
    print(f"求解失败: {sol.message}")

# 绘图
plt.figure(figsize=(12, 6))
plt.subplot(1,2,1)
plt.plot(sol.x, sol.y[0], label='g(r)')
plt.xlabel('r')
plt.ylabel('g(r)')
plt.title('g(r)的解')
plt.legend()
plt.grid()

plt.subplot(1,2,2)
plt.plot(sol.x, sol.y[2], label='f(r)')
plt.xlabel('r')
plt.ylabel('f(r)')
plt.title('f(r)的解')
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

关键说明 #

• 初始猜测的重要性：非零的初始猜测是避免平凡解的核心，需根据预期解的形状构造（比如高斯峰、指数趋近函数）。
• 边界条件的调整：如果近端的正则性条件不符合你的问题，可以修改bc函数中的近端条件（比如换成$g(r0)=0$），需结合物理背景调整。
• 数值精度优化：可调整solve_bvp的tol、max_nodes参数来提升求解精度和收敛性。

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