我正在尝试制作一个脚本来计算使用sympy的简单谐振子的特征状态/能级。
这是一个表现出周期性运动的系统，如弹簧上的质量。一般来说，简谐振子的运动方程可以写成:
F = -kx
地点:F
- 为作用在振荡器上的力。
- k为弹簧常数。
- x为离平衡位置的位移。
简谐振子的本征值是系统的可能能级。对于简谐振子，本征值由
给出:E_n = (n + 1/2) * h * nu
地点:
- n是一个整数，表示能量级别
- h普朗克常数
- ν为振荡器的频率，
这些特征值是等间隔的，形成一个离散谱，这是量子系统的一个特征属性。
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)
，其中代码定义位移(x)、弹簧常数(k)、质量(m)、普朗克常数(h)和频率(nu)作为符号变量，然后用这些变量定义简谐振子的运动方程。
然而，我得到这个错误，无法找到一个方法来调试它。
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
感谢您的帮助。
打印出表达式中的不同数量，以调试正在发生的事情。你的方程是
In [5]: eq
Out[5]: -k⋅x = 0
我猜那不是你想要的。原因如下:
In [6]: diff(x, t, 2)
Out[6]: 0
这是因为x作为一个符号不依赖于t，所以就wrt与t的微分而言，它是一个常数，其导数为零。您应该使x成为t的函数:
In [7]: x = Function('x')
In [10]: eq = Eq(-k*x(t), m*diff(x(t), t, 2))
In [11]: eq
Out[11]:
2
d
-k⋅x(t) = m⋅───(x(t))
2
dt
In [12]: dsolve(eq)
Out[12]:
_____ _____
╱ -k ╱ -k
-t⋅ ╱ ─── t⋅ ╱ ───
╲╱ m ╲╱ m
x(t) = C₁⋅ℯ + C₂⋅ℯ
如果你希望在那里看到sin和cos，那么问题是你没有说任何关于符号k和m:
In [13]: k, m = symbols('k, m', positive=True)
In [14]: eq = Eq(-k*x(t), m*diff(x(t), t, 2))
In [15]: dsolve(eq)
Out[15]:
⎛√k⋅t⎞ ⎛√k⋅t⎞
x(t) = C₁⋅sin⎜────⎟ + C₂⋅cos⎜────⎟
⎝ √m ⎠ ⎝ √m ⎠