使用solve_ivp求解电子基态能量时遇到数值异常问题

你遇到的核心问题是初始条件错误，同时数值积分的稳定性也需要优化，以下是具体分析和修正方案：

问题重现 #

尝试用solve_ivp求解电子在谐振子势能中的基态能量，薛定谔方程为：

ψ''(x) = -2 * m_e*(E - V(x)) * ψ(x)
初始条件为ψ(-10a)=0，ψ'(-10a)=1，势能V(x)=V₀*(x/a)²。采用自然单位定义常数，目标是找到E使得ψ(10a)=0。纸笔计算得E≈138.0214281，但代入代码后ψ(10a)结果为1.136e+142，解完全不符合预期。

关键错误点 #

1. 初始条件完全错误
   代码中y0 = [-10 * a, 1]违反了题目要求的ψ(-10a)=0，第一个元素应该是0而非-10*a。这个错误直接导致解从初始时刻就偏离正确轨道，最终指数爆炸。

2. 刚性问题的数值稳定性不足
   当E低于x=±10a处的势能时，方程属于刚性系统（解包含快速增长的指数项），默认的RK45方法难以稳定求解，建议改用适合刚性问题的方法（如Radau）。

修正后的代码 #

from scipy.constants import c, m_e, eV, e
from scipy.integrate import solve_ivp
import numpy as np

me_ev = m_e * c**2 / e        # 电子质量（能量单位）
l_bohr = 5.2918e-11           # 玻尔半径[m]
convert_length = 5067730.179  # [m]转自然单位长度[1/eV]
V0 = 50
a = 10e-11 * convert_length

def V(x):
    return V0 * ((x / a) ** 2)

def eqn2(x, y, energy):
    psi, v = y
    dpsidx = v
    dvdx = -2 * me_ev * psi * (energy - V(x))
    return dpsidx, dvdx

def solve2(energy):
    # 修正初始条件：ψ(-10a)=0，ψ'(-10a)=1
    y0 = [0, 1]
    x_range = [-10 * a, 10 * a]
    sol = solve_ivp(
        eqn2,
        x_range,
        y0,
        t_eval=np.linspace(x_range[0], x_range[1], 1000),
        args=(energy,),
        dense_output=True,
        method='Radau'  # 改用适合刚性问题的求解器
    )
    return sol.sol(x_range[1])[0]

# 测试修正后的代码
print(f"ψ(10a) = {solve2(138.0214281)}")

更优方案：使用边界值求解器solve_bvp #

由于这是典型的两点边值问题（两端ψ=0），solve_ivp（初值问题求解器）并非最优选择，solve_bvp专门处理这类问题，数值稳定性更强：

from scipy.integrate import solve_bvp

def bvp_eqn(x, y, energy):
    psi, v = y
    return np.vstack([v, -2 * me_ev * psi * (energy - V(x))])

def bvp_bc(ya, yb, energy):
    # 边界条件：ψ(-10a)=0，ψ(10a)=0
    return np.array([ya[0], yb[0]])

# 构造初始猜测解（基态波函数为高斯型）
x_guess = np.linspace(-10*a, 10*a, 100)
y_guess = np.zeros((2, x_guess.size))
y_guess[0] = np.exp(-(x_guess/a)**2)

# 求解边界值问题
sol_bvp = solve_bvp(
    lambda x, y: bvp_eqn(x, y, 138.0214281),
    lambda ya, yb: bvp_bc(ya, yb, 138.0214281),
    x_guess, y_guess
)

print(f"BVP求解结果：ψ(10a) = {sol_bvp.sol(10*a)[0]}")
print(f"求解状态：{'成功' if sol_bvp.success else '失败'}，信息：{sol_bvp.message}")

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