SymPy求解振荡运动方程报错：仅支持单变量函数问题求助

以下是用于求解振荡系统运动方程的代码：

#import libraries
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from sympy import *
import sympy as sp
from sympy.plotting import plot as symplot

#initialize symbols
delta1,delta2,mass,youngs,width,thickness,length = sp.symbols("delta_1 delta_2 m E b t L") #spring
permeability,magnetization,radius,height = sp.symbols("mu_0 M R h") #magnetic replusion
rho,dragconst = sp.symbols("rho Cd")

#EOM
EOM=Eq(mass*Derivative(delta1, dt,2)-((delta1**2)*youngs*width*thickness**3)/4*length**3-((pi*permeability*(magnetization**2)*radius**4*((delta2-delta1)**4)+(10*(delta2-delta1)**2)*thickness**2+(4*(delta2-delta1)**3*thickness+12*(delta2-delta1)*thickness**3+4*thickness**4))/((4*(delta2-delta1)**2)*(((delta2-delta1)+2*thickness)**2)*((delta2-delta1)+thickness)**2))+Rational(-1, 2)*rho*(Derivative(delta1, dt,1)**2)*dragconst*width*length,0)
display(EOM)
dsolve(EOM,delta1)

运行后输出的运动方程与错误信息：

−𝐶𝑑𝐿𝑏𝜌(𝑑𝑑𝑡𝛿1)22−𝐸𝐿3𝑏𝛿21𝑡34+𝑚𝑑2𝑑𝑡2𝛿1−𝜋𝑀2𝑅4𝜇0(−𝛿1+𝛿2)4+4𝑡4+𝑡3(−12𝛿1+12𝛿2)+10𝑡2(−𝛿1+𝛿2)2+4𝑡(−𝛿1+𝛿2)34(−𝛿1+𝛿2)2(−𝛿1+𝛿2+𝑡)2(−𝛿1+𝛿2+2𝑡)2=0
 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_24300/1246088876.py in <module>
     14 EOM=Eq(mass*Derivative(delta1, dt,2)-((delta1**2)*youngs*width*thickness**3)/4*length**3-((pi*permeability*(magnetization**2)*radius**4*((delta2-delta1)**4)+(10*(delta2-delta1)**2)*thickness**2+(4*(delta2-delta1)**3*thickness+12*(delta2-delta1)*thickness**3+4*thickness**4))/((4*(delta2-delta1)**2)*(((delta2-delta1)+2*thickness)**2)*((delta2-delta1)+thickness)**2))+Rational(-1, 2)*rho*(Derivative(delta1, dt,1)**2)*dragconst*width*length,0)
     15 display(EOM)
---> 16 dsolve(EOM,delta1)
     17 
     18 

~\anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    602 
    603         # See the docstring of _desolve for more details.
--> 604         hints = _desolve(eq, func=func,
    605             hint=hint, simplify=True, xi=xi, eta=eta, type='ode', ics=ics,
    606             x0=x0, n=n, **kwargs)

~\anaconda3\lib\site-packages\sympy\solvers\deutils.py in _desolve(eq, func, hint, ics, simplify, prep, **kwargs)
    207     # recursive calls.
    208     if kwargs.get('classify', True):
--> 209         hints = classifier(eq, func, dict=True, ics=ics, xi=xi, eta=eta,
    210         n=terms, x0=x0, hint=hint, prep=prep)
    211 

~\anaconda3\lib\site-packages\sympy\solvers\ode\ode.py in classify_ode(eq, func, dict, ics, prep, xi, eta, n, **kwargs)
    940 
    941     if func and len(func.args) != 1:
--> 942         raise ValueError("dsolve() and classify_ode() only "
    943         "work with functions of one variable, not %s" % func)
    944 

ValueError: dsolve() and classify_ode() only work with functions of one variable, not delta_1

错误原因 #

核心问题是你将delta1定义为普通符号，而非关于时间t的单变量函数。SymPy的dsolve要求待求解的未知量必须是单变量函数（此处应为t的函数），而非独立符号。此外原方程中弹簧项的括号位置错误，会导致物理量纲不符，且存在符号命名冲突（厚度符号t与时间变量重名）。

修改步骤 #

• 先定义时间符号t，再将delta1声明为t的函数：delta1 = sp.Function('delta_1')(t)
• 修正厚度的符号名，避免和时间t冲突（如改为t_s）
• 修正弹簧项的括号，确保分母是4*length**3（原代码中误写为除以4后乘以length³）
• 导数写法可简化为delta1.diff(t, 2)（二阶导数）和delta1.diff(t)（一阶导数），更直观

修正后的完整代码 #

#import libraries
from sympy.interactive import printing
printing.init_printing(use_latex=True)
import sympy as sp
from sympy.plotting import plot as symplot

#initialize symbols
t = sp.symbols('t')  # 先定义时间变量
delta2,mass,youngs,width,thickness,length = sp.symbols("delta_2 m E b t_s L")  # 厚度符号改为t_s，避免和时间t冲突
permeability,magnetization,radius,height = sp.symbols("mu_0 M R h")
rho,dragconst = sp.symbols("rho Cd")

# 将delta1定义为t的函数
delta1 = sp.Function('delta_1')(t)

# 修正后的EOM
EOM = sp.Eq(
    mass * delta1.diff(t, 2) 
    - (youngs * width * thickness**3 * delta1**2) / (4 * length**3)  # 修正分母括号
    - (
        sp.pi * permeability * magnetization**2 * radius**4 * (delta2 - delta1)**4 
        + 10 * (delta2 - delta1)**2 * thickness**2 
        + 4*(delta2-delta1)**3*thickness + 12*(delta2-delta1)*thickness**3 + 4*thickness**4
    ) / (
        4 * (delta2 - delta1)**2 * ((delta2 - delta1) + 2*thickness)**2 * ((delta2 - delta1) + thickness)**2
    )
    + sp.Rational(-1, 2) * rho * (delta1.diff(t))**2 * dragconst * width * length,
    0
)
display(EOM)

# 尝试求解
try:
    solution = sp.dsolve(EOM, delta1)
    display(solution)
except Exception as e:
    print(f"求解失败：{str(e)}")

补充说明 #

这个方程是非线性二阶常微分方程，包含δ₁的平方项、速度平方项（阻尼）和复杂的磁斥力非线性项，SymPy通常无法给出解析解。如果需要求解，建议：

• 对系统进行线性化近似（如在平衡位置附近展开），再尝试解析求解
• 使用数值解法，将方程转换为一阶方程组后，用SciPy的scipy.integrate.solve_ivp进行数值积分

内容的提问来源于stack exchange，提问作者This is Mark