Sympy求解简谐振子本征态报错：无法自动检测函数

简谐振子基础 #

简谐振子是一类做周期性运动的系统（比如弹簧振子），经典运动方程由胡克定律给出：

F = -kx

其中：

• F：振子所受的力
• k：弹簧常数
• x：振子偏离平衡位置的位移

量子简谐振子的本征值（即系统可能的能级）公式为：

Eₙ = (n + 1/2) * h * ν

其中：

• n：表示能级的非负整数（n=0,1,2,...）
• h：普朗克常数
• ν：振子的固有频率

这些能级呈等间距分布，构成离散谱，是量子系统的典型特征。

问题代码与报错 #

我尝试用Sympy编写脚本计算上述本征值，但运行报错，代码如下：

import sympy as sym

# Define the variables and the equation of motion
x, k, m, h, nu,t = sym.symbols('x k m h nu t')
eq = sym.Eq(-k * x, m * sym.diff(x,t, 2))

# Solve the equation of motion for the displacement x(t)
sol = sym.dsolve(eq)

# Compute the eigenvalues of the oscillator
E = sym.solve(sym.det(sol.lhs - nu * sym.eye(2)), nu)

# Print the eigenvalues
print(E)

运行时抛出如下错误：

ValueError                                Traceback (most recent call last)
Input In [18], in <cell line: 8>()
      5 eq = sym.Eq(-k * x, m * sym.diff(x,t, 2))
      7 # Solve the equation of motion for the displacement x(t)
----> 8 sol = sym.dsolve(eq)
     10 # Compute the eigenvalues of the oscillator
     11 E = sym.solve(sym.det(sol.lhs - nu * sym.eye(2)), nu)

File /opt/python/3.9/envs/pulsar/lib/python3.10/site-packages/sympy/solvers/ode/ode.py:605, in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    602 given_hint = hint  # hint given by the user
    604 # See the docstring of _desolve for more details.
--> 605 hints = _desolve(eq, func=func,
    606     hint=hint, simplify=True, xi=xi, eta=eta, type='ode', ics=ics,
    607     x0=x0, n=n, **kwargs)
    608 eq = hints.pop('eq', eq)
    609 all_ = hints.pop('all', False)

File /opt/python/3.9/envs/pulsar/lib/python3.10/site-packages/sympy/solvers/deutils.py:180, in _desolve(eq, func, hint, ics, simplify, prep, **kwargs)
    178 # preprocess the equation and find func if not given
    179 if prep or func is None:
--> 180     eq, func = _preprocess(eq, func)
    181     prep = False
    183 # type is an argument passed by the solve functions in ode and pde.py
    184 # that identifies whether the function caller is an ordinary
    185 # or partial differential equation. Accordingly corresponding
    186 # changes are made in the function.

File /opt/python/3.9/envs/pulsar/lib/python3.10/site-packages/sympy/solvers/deutils.py:82, in _preprocess(expr, func, hint)
     80     funcs = set().union(*[d.atoms(AppliedUndef) for d in derivs])
     81     if len(funcs) != 1:
--> 82         raise ValueError('The function cannot be '
     83             'automatically detected for %s.' % expr)
     84     func = funcs.pop()
     85 fvars = set(func.args)

ValueError: The function cannot be automatically detected for -k*x

错误原因与修复方案 #

错误根源 #

1. 函数定义问题：当前代码里的x是普通符号，不是关于t的函数，Sympy无法识别diff(x,t,2)中的待解函数，导致dsolve无法自动检测函数，抛出错误。
2. 逻辑偏差：代码试图通过解经典运动方程来求量子能级，这是方法错误——经典运动方程的解是位移随时间的变化，和量子本征值的求解逻辑完全不同。量子简谐振子的能级需要通过求解薛定谔方程得到。

修正后的代码 #

方案1：直接利用能级公式计算

import sympy as sym

# 定义符号
n, h, nu = sym.symbols('n h nu', integer=True, nonnegative=True)

# 量子简谐振子能级公式
E_n = (n + sym.Rational(1, 2)) * h * nu

# 打印前5个能级（n=0到4）
for level in range(5):
    print(f"n={level}: E = {E_n.subs(n, level)}")

方案2：求解薛定谔方程推导能级

import sympy as sym

# 定义符号与波函数
x, m, h_bar, E = sym.symbols('x m h_bar E')
psi = sym.Function('psi')(x)

# 定义角频率ω（ω=√(k/m)）
omega = sym.symbols('omega')
# 量子简谐振子的薛定谔方程
schrodinger_eq = sym.Eq(
    -h_bar**2/(2*m) * sym.diff(psi, x, 2) + sym.Rational(1,2)*m*omega**2*x**2*psi,
    E*psi
)

# 求解方程（Sympy会返回包含能级条件的解）
solution = sym.dsolve(schrodinger_eq, psi)
# 直接生成标准能级表达式
n = sym.symbols('n', integer=True, nonnegative=True)
E_n = (n + sym.Rational(1,2)) * h_bar * omega
print("量子简谐振子能级：", E_n)

运行说明 #

• 方案1直接利用已知的能级公式计算，简单高效，适合快速得到结果。
• 方案2通过求解薛定谔方程，更贴近量子力学的理论推导，能同时得到波函数形式与能级。

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