python单自由度振动计算-数据可视化2

参考

单自由度振动理论推导参考：

结果

输出位移由3部分组成：

代码

import matplotlib.pyplot as plt
import numpy as np
from moviepy.video.io.bindings import mplfig_to_npimage
import moviepy.editor as mpy

from numpy import pi,sqrt,cos,sin,arctan,exp

##########begin:计算单自由度振动###############
k1=(2*pi)**2
m1=1
wn=sqrt(k1/m1)
s_=0.16
wd=sqrt(1-s_**2)*wn
c1=2*s_*wn*m1
t=np.array(np.linspace(0,4,1000,False))#0:0.002:2;
f0=2
f_x0=1*sin(2*pi*f0*t)
f1=f_x0
s=2*pi*f0/wn;
B0=1/k1#%B0=F0/k1 2019.3.31 by lijilin
B=B0/sqrt((1-s**2)**2+(2*s_*s)**2);
phai=arctan(2*s_*s/(1-s**2));
#%%%%%%%%%%%%%%%%%%
if s>1:
#    B=-B;
    phai=phai+pi;#%%%%当频率比s>1,响应滞后相位角为pi/2~pi!!!!!2019.3.31bylijilin
#%%%%%%%%%%%%%%%%%%%%%%
x0=0.02;
dx0=0.0;
x1=exp(-s_*wn*t)*(x0*cos(wd*t)+(dx0+s_*wn*x0)/wd*sin(wd*t))
x2=B*exp(-s_*wn*t)*(sin(phai)*cos(wd*t)+(s_*wn*sin(phai)-2*pi*f0*cos(phai))/wd*sin(wd*t))
x3=B*sin(2*pi*f0*t-phai)
x=x1+x2+x3
plt.figure
plt.subplot(3,1,1)
plt.plot(t,f1)
plt.subplot(3,1,2)
plt.plot(t,x)
plt.grid


# 用matplotlib绘制一个图形
plt.figure
duration = 4

fig_mpl, ax = plt.subplots(1,figsize=(5,3), facecolor='white')

fs=1000/4

xx=lambda d:t[:int(d*fs+1)]
zz=lambda d:x[:int(d*fs+1)]
zz1=lambda d:x1[:int(d*fs+1)]
zz2=lambda d:x2[:int(d*fs+1)]
zz3=lambda d:x3[:int(d*fs+1)]


ax.set_title("output")
ax.set_xlabel("t")
ax.set_ylim(-0.03,0.03)
ax.set_xlim(0,4)
line, = ax.plot(xx(0),zz(0), lw=3)
line1, = ax.plot(xx(0),zz1(0), lw=2)
line2, = ax.plot(xx(0),zz2(0), lw=2)
line3, = ax.plot(xx(0),zz3(0), lw=2)

# 用MoviePy制作动（为每个t更新曲面）。制作一个GIF

def make_frame_mpl(t):
    line.set_xdata(xx(t))
    line1.set_xdata(xx(t))
    line2.set_xdata(xx(t))
    line3.set_xdata(xx(t))
    line.set_ydata(zz(t) ) # 更新曲面
    line1.set_ydata(zz1(t) ) # 更新曲面
    line2.set_ydata(zz2(t) ) # 更新曲面
    line3.set_ydata(zz3(t) ) # 更新曲面
    
    return mplfig_to_npimage(fig_mpl) # 图形的RGB图像

animation =mpy.VideoClip(make_frame_mpl, duration=duration)
animation.write_gif("sin_mpl.gif", fps=20)