使用solve_ivp求解TISE后调用fsolve寻找基态能量时遇数组形状不匹配错误的解决方案咨询

你好，我看了你的代码和遇到的错误，这个数组形状不匹配的问题主要是因为fsolve在迭代过程中可能会传递数组类型的energy参数给你的solve函数，而你的微分方程求解逻辑是基于标量能量设计的，导致内部运算出现形状冲突。另外，你已经导入了更适合标量根查找的root_scalar，用它替代fsolve会更稳妥，也能直接避开这个问题。

错误原因详解 #

fsolve是为多变量优化设计的工具，它在计算数值导数时可能会传递数组作为输入参数。而你的solve函数和eqn微分方程都是针对标量能量编写的，当energy变成数组时，eqn里的-2 * me_ev * energy * y[0]会生成数组结果，和原本应该是标量的y[0]运算后，返回的导数结果形状就会和solve_ivp期望的长度为2的一维数组不一致，从而抛出形状不匹配的错误。

修正后的代码及说明 #

from scipy.constants import c, m_e, e
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import root_scalar

# 全局定义常量，避免重复计算，提升可读性
me_ev = m_e * c**2 / e          # 电子质量，单位[eV]
l_bohr = 5.2918e-11             # 玻尔半径，单位[m]
convert_length = 5067730.179    # 长度单位转换：[m] → [1/eV]
l_ev = l_bohr * convert_length  # 转换后的区间长度L，单位[1/eV]

# 一维TISE微分方程：y = [phi, dphi/dx]，返回[dphi/dx, d²phi/dx²]
def eqn(x, y, energy):
    f0 = y[1]  # dphi/dx
    f1 = -2 * me_ev * energy * y[0]  # d²phi/dx²
    return np.array([f0, f1])

# 给定能量，求解TISE并返回phi(L)的标量值
def solve(energy):
    # 初始条件：phi(0)=0，dphi(0)/dx=1
    init_cond = np.array([0.0, 1.0])
    # 启用dense_output获取连续解，更高效地计算L点的phi值
    sol = solve_ivp(eqn, (0., l_ev), init_cond, args=(energy,), dense_output=True)
    # 直接计算L位置的phi值
    phi_L = sol.sol(l_ev)[0]
    return phi_L

# 用root_scalar找基态能量：目标是让solve(energy)=0
# 选择brentq方法，需要给出包含根的大致区间（基于你提供的134初始值设置）
result = root_scalar(solve, bracket=[100, 150], method='brentq')

if result.converged:
    print(f"找到基态能量：{result.root:.4f} eV")
else:
    print("根查找未收敛，请调整搜索区间或尝试其他方法")

额外优化点 #

1. 把重复定义的常量移到全局，避免每次调用solve时重复计算，提升运行效率。
2. 用solve_ivp的dense_output=True生成连续解，直接计算L点的phi值，比生成采样点再取最后一个值更准确高效。
3. root_scalar的brentq方法稳定性强，只需要提供包含根的区间即可，你可以根据需求调整区间范围来寻找更高能级。
4. 简化了solve函数的参数传递，直接在内部调用eqn，逻辑更清晰。

这样修改后，就不会再出现数组形状不匹配的问题，根查找的稳定性也会更好。

备注：内容来源于stack exchange，提问作者Emily