使用 sympy 求解阻尼震动方程的解析解
|
810
|
0
|
未分类

695 字
|
6 分钟

阻尼微分方程的形式

考虑一维情形下的阻尼震动，考虑宏观低速情况，依据牛顿第二定律：$$F=ma=m\frac{\mathrm{d}^2x}{\mathrm{d}t^2}$$其中，$m$为物体质量，$x=x(t)$为物体位置矢量（下同），$t$为时间。以及一维情形下的阻尼振动，物体所受的合外力：$$F = -cv-kx=-c\frac{\mathrm{d}x}{\mathrm{d}t}-kx$$其中$c$为阻尼系数，$k$为回复系数。将上面两个方程联立，整理得到：$$m\frac{\mathrm{d}^2x}{\mathrm{d}t^2} = -c\frac{\mathrm{d}x}{\mathrm{d}t}-kx$$微分方程合写为：

$$F(t) =m\frac{\mathrm{d}^2x}{\mathrm{d}t^2} +c\frac{\mathrm{d}x}{\mathrm{d}t}+kx = 0$$

这里还需要强调，我们最终需要求解的是$x(t)$的表达式，不必拘泥于$F(t)$这一表示形式。对于阻尼振动，从物理上考虑，解具有的形式为：$$x(t) = Ae^{i(\omega_0 t+\phi)-\alpha t}+C$$其中，$\alpha$为阻尼比，整理得到的$Ae^{-\alpha t}$为振幅在时间上的衰减，即能量的耗散。$\omega_0$为振动角频率。$\phi$为振动初相。

考虑实际物理意义以及便于求解，可以做如下代换：$$\omega_0 = \sqrt{\frac{k}{m}}~~~~~~\zeta = \frac{c}{2\sqrt{mk}} $$其中，$\omega_0$为系统不受阻尼力时的振动角频率，称为自然角频率。$\zeta$为阻尼比，当$\zeta > 1$时，对应过阻尼振动，系统振动趋于平稳，即$\lim_{t\to +\infty}x(t) = 0$，$\zeta = 1$时为临界阻尼情况，系统达到平衡状态时不再振动。$0< \zeta < 1$时为欠阻尼振动，系统仍然保持振动状态，但角频率改变。

依据上述代换，微分方程的形式变为：$$\frac{\mathrm{d}^2 x}{\mathrm{d}t^2} + 2\zeta\omega_0\frac{\mathrm{d}x}{\mathrm{d}t}+\omega_0^2x = 0$$后面，我们将根据这个方程的形式来求解。

使用sympy求解阻尼振动方程

依据微分方程的形式，考虑初始条件$v_0 = 1$，$x_0 = 2$，$\omega_0=0.5\pi$。使用sympy建立阻尼振动方程的代码为：

import sympy
import numpy as np
import matplotlib.pyplot as plt

x = sympy.Function('x')         # 定义位移矢量函数
t = sympy.symbols('t')          # 定义时间变量

zeta, omega_0 = sympy.symbols('zeta, omega_0', real=True, positive=True)
# 定义阻尼系数和固有频率

v0, x0 = sympy.symbols('v0,x0') # 定义初始条件
initcons = {x(0): x0,                       # x(0) = x0
            x(t).diff(t).subs(t, 0): v0}    # (dx/dt)|t=0 = v0
# 方程初始条件，以字典形式存储

ode = sympy.diff(x(t),t,2)+2*zeta*omega_0*sympy.diff(x(t),t)+omega_0**2*x(t)
# 定义二阶常微分方程

ode_res1 = sympy.dsolve(ode, x(t), ics=initcons)    # 求解二阶常微分方程，非临界条件
xt_critical = sympy.limit(ode_res1.rhs, zeta, 1)    # 求解二阶常微分方程，临界条件

tt = np.linspace(0,10,1000)                         # 绘图用时间
# 考虑 omega_0 = 0.5*np.pi, x0 = 2, v0 = 1, zeta为变量
ode_res1 = ode_res1.subs({omega_0:0.5*np.pi, x0:2, v0:1})           # 非临界条件替换参数
xt_critical = xt_critical.subs({omega_0:0.5*np.pi, x0:2, v0:1})     # 临界条件替换参数

zts = [0.1, 0.5, 1.0, 2, 5]                 # 不同zeta值，用于绘图

plt.figure(dpi = 400)
for zt in zts:
    if zt == 1:
        xt = sympy.lambdify(t, xt_critical.subs({zeta: zt}), 'numpy')
        # 替换参数并实例化，得到临界条件函数表达式 xt(t)
    else:
        xt = sympy.lambdify(t, ode_res1.rhs.subs({zeta: zt}), 'numpy')
        # 替换参数并实例化，得到非临界条件函数表达式 xt(t)
    plt.plot(tt, xt(tt), label=r"$\zeta = {:.1f}$".format(zt),linewidth=0.7)

plt.legend()
plt.xlabel(r"$t/s$")
plt.ylabel(r"$x(t)$")
plt.tight_layout()

plt.savefig("阻尼震荡.png", dpi=400)
plt.show()

过程中，我们先使用sympy对方程形式进行了输入，并对方程的临界形式和非临界形式进行了求解，对求解结果实例化，并依据不同$\zeta$进行了代换。

上一篇