利用python自带的库sympy，求解不同阻尼比的振动方程表达式

#利用python自带的库sympy，求解不同阻尼比的振动方程表达式
import sympy #导入数学符号运算库，属于自带库
import numpy as np
import matplotlib.pyplot as plt #导入绘图模块库，是外部库
plt.rcParams['font.sans-serif'] = ['Microsoft Yahei']
plt.rcParams['axes.unicode_minus'] = False #绘图可以中文显示

t,x0,v0 = sympy.symbols('t,x0,v0',positive=True)    #定义变量和初始值，都为实数
omega_n = sympy.Symbol('omega',positive=True) 
zeta = sympy.Symbol('zeta',positive=True) 
#定义频率为实数和正数，可以把指数表达式化简为三角函数表达式;W为变量名，omega_n为变量
x = sympy.Function('x')     #定义位移函数
ode_SHV = sympy.diff(x(t),t,2)+2*zeta*omega_n*sympy.diff(x(t),t)+omega_n**2*x(t)      #定义微分方程
initcons = {x(0):x0,sympy.diff(x(t),t).subs(t,0):v0}   #定义初始条件，数据结构为字典型
ode_res0 = sympy.dsolve(ode_SHV)                #通解形式
ode_res1 = sympy.dsolve(ode_SHV,ics=initcons)   #代入初始条件确定常数C1,C2
xt_critical = sympy.limit(ode_res1.rhs, zeta, 1) 
#临界阻尼通过极限方式求解表达式，不能直接代入值，会出现分母为0

sympy.pprint(ode_res0)   #手写输出形式
print('\n')
sympy.pprint(xt_critical)   #手写输出形式
tt = np.linspace(0,15,500)
wn = 0.5*np.pi
x00,v00 = [1,0]

fig = plt.figure("有阻尼简谐振动",figsize=(12, 8))   #定义画布
ax =fig.add_subplot(1,1,1)                          #定义子图
plt.title("单自由度有阻尼简谐振动")
for zt in [0.1, 0.5, 1, 2, 5]:
    if zt==1:
        xt = sympy.lambdify(t,xt_critical.subs({omega_n:wn, x0:x00, v0:v00}))
    else:
        xt = sympy.lambdify(t,ode_res1.rhs.subs({omega_n:wn, x0:x00, v0:v00, zeta:zt}))
    ax.plot(tt, xt(tt).real, label=r"$\zeta = %.1f$" % zt)  #绘制响应曲线
ax.set_xlabel(r"$t$", fontsize=18)
ax.set_ylabel(r"$x(t)$", fontsize=18)
ax.legend()
plt.show()