使用 SymPY dsolve 求解临界阻尼振荡器的微分方程时出现的问题

Question

我一直在尝试编写一个程序，该程序将为谐振子取一组值和初始条件，求解微分方程并显示 phi(t) 的函数。

import sympy as sp
from matplotlib import pyplot as plt
import numpy as np

t = sp.Symbol("t")
phi = sp.Function('phi')

#Values of the oscilator
beta = 0.4
omega_0 = 0.4
f0 = 0
omega = 0

#Initial conditions
phi_0 = 3.14/4
v_phi_0 = 0

#phi''(t) + 2*beta*phi'(t) + omega_0^2*phi(t) - f0*sin(omega*t) = 0
eq = sp.diff(phi(t),t,2)+2*beta*sp.diff(phi(t),t)+omega_0**2*phi(t)-f0*sp.sin(omega*t)
print(sp.dsolve(eq))
equation = sp.dsolve(eq,ics={phi(0):phi_0,sp.diff(phi(t),t).subs(t,0):v_phi_0})
print(equation)
wzor = sp.lambdify(t,equation.rhs,"numpy")
xvals = np.arange(0,30,.1)
yvals = wzor(xvals)

# make figure
fig, ax = plt.subplots(1,1)
ax.plot(xvals, yvals)
ax.set_xlabel('t')
ax.set_ylabel('phi(t)')
plt.show()

一切似乎都很好，但是当beta和omega_0的值相等时（所谓的临界阻尼）和低于1, 程序好像解方程有误

Eq(phi(t), (C1*sin(3.65002571393103e-9*t) + C2*sin(3.6500259065127e-9*t) + C3*cos(3.65002571393103e-9*t) + C4*cos(3.6500259065127e-9*t))*exp(-0.4*t))

其他情况，包括beta和omega_0等于且大于1时，解微分方程后只有C1和C22个常数，但在本例中有 4 个，这是不应该发生的。因此，程序无法继续绘制绘图，因为这 2 个额外的常量我无法求解。

ax.plot(xvals, yvals)

TypeError: can't convert expression to float

编辑：另外，当施加驱动力时，f0 和 omega 与 0 不同，出现以下错误

NotImplementedError: Cannot find 2 solutions to the homogeneous equation necessary to apply undetermined coefficients to 0.16*phi(t) - 5*sin(0.4*t) + 0.8*Derivative(phi(t), t) + Derivative(phi(t), (t, 2)) (number of terms != order)

这可能与主要问题有关，但我认为提一下是件好事。

编辑 2：当 beta 和 omega_0 相等且大于 3，但也是十进制时，其他问题开始出现。当它们是整数时没有问题。

Answer 1

符号计算和浮点常量不能很好地结合在一起。一种简单的解决方法是使用 sympy Rational 数据类型

将所有常量替换为有理数

#Values of the oscilator
beta = sp.Rational(2,5)
omega_0 = sp.Rational(2,5)
f0 = 0
omega = 0

#Initial conditions
phi_0 = sp.Rational(314,400)
v_phi_0 = 0

仅进行这些更改，输出就正确了

Eq(phi(t), (C1 + C2*t)*exp(-2*t/5))
Eq(phi(t), (157*t/500 + 157/200)*exp(-2*t/5))

情节生成为

使用 SymPY dsolve 求解临界阻尼振荡器的微分方程时出现的问题

Problem while using SymPY's dsolve to solve differential equation of criticaly damped oscilator

python

sympy

ode