摆轮运动的振幅（角度）和角速度随时间变化的曲线的matlab程序

[程序]ltg1main.m如下

%摆轮在平衡点从静止开始的受迫振动曲线的主程序(用解析式)

clear,b=0.1;W=2;tm=60;%清除变量,阻尼系数参考值,驱动力参考值,最大时间

ltg1fun(b,W,tm)%调用函数文件

ltg1fun(b,1/2,tm)%调用函数文件

ltg1fun(b,1,tm)%调用函数文件

W=sqrt(1-b^2);%稳态振幅最大的圆频率

ltg1fun(b,W,tm)%调用函数文件函数文件

ltg1fun.m

%物体在平衡点从静止开始的受迫振动曲线的函数文件(用解析式)

function fun(b,W,tm)

if W==1 W=1-eps;end%如果为1则改小一点

W1=sqrt(1-b^2);%约化阻尼圆频率

t=linspace(0,tm,1000);%时间向量

TH1=sqrt(W1^2*(W^2-1)^2+b^2*(W^2+1)^2)/W1/((W^2-1)^2+4*b^2*W^2);%阻尼振幅

PHI1=atan2(b*(W^2+1),W1*(W^2-1));%阻尼振动初相

th1=TH1*exp(-b*t).*cos(W1*t+PHI1);%阻尼振动角度

TH2=1/sqrt((W^2-1)^2+4*b^2*W^2);%等幅振动振幅

PHI2=atan2(-2*b*W,1-W^2);%等幅振动初相

th2=TH2*cos(W*t+PHI2);%等幅振动角度

figure,plot(t,th1,'--b',t,th2,'-.g',t,th1+th2,'r','LineWidth',2)%创建图形窗口,画曲线

legend('\it\theta\rm_1/\it\theta\rm_0','\it\theta\rm_2/\it\theta\rm_0',...

    '\it\theta\rm/\it\theta\rm_0')%加图例

txt=['\it\beta/\omega\rm_0=',num2str(b)];%阻尼因子字符串

txt=[txt,',\it\Omega/\omega\rm_0=',num2str(W)];%连接驱动力圆频率字符串

grid on,fs=14;%加网格,字体大小

title(['摆轮振动的角度(',txt,')'],'FontSize',fs)%标题

xlabel('\it\omega\rm_0\itt','FontSize',fs)%标记横坐标

hold on,plot([0,tm],[TH2,TH2],':','LineWidth',2)%保持属性,画水平线

text(tm,TH2,num2str(TH2))%标记稳态角度

w1=-b*th1-W1*TH1*exp(-b*t).*sin(W1*t+PHI1);%暂态角速度

w2=-W*TH2*sin(W*t+PHI2);%稳态角速度

figure,plot(t,w1,'--',t,w2,'-.',t,w1+w2,'LineWidth',2)%创建图形窗口,画曲线

legend('\it\omega\rm_1/\it\omega\rm_0\it\theta\rm_0',...

   '\it\omega\rm_2/\it\omega\rm_0\it\theta\rm_0',...

   '\it\omega\rm/\it\omega\rm_0\it\theta\rm_0')%加图例

title(['摆轮振动的角速度(',txt,')'],'FontSize',fs)%标题

xlabel('\it\omega\rm_0\itt','FontSize',fs),grid on%标记横坐标,加网格

W2=W*TH2;%稳态角速度幅度

hold on,plot([0,tm],[W2,W2],':','LineWidth',2)%保持属性,画水平线

text(tm,W2,num2str(W2))%标记稳态角速度幅度

不同阻力矩下受迫振动的幅频特性与相频特性matlab程序

%受迫振动的稳态位移振幅,初相和速度振幅曲线

clear,b=1:-0.2:0;n=length(b);%清除变量,约化阻尼因子向量,向量的个数

wm=2.5;w=0:0.05:wm;w(w==1)=1+eps;%最大约化驱动角频率,约化驱动角频率,为1者取大一点

[B,W]=meshgrid(b,w);%阻尼因子和驱动角频率矩阵

TH2=1./sqrt((1-W.^2).^2+(2*B.*W).^2);figure%约化位移振幅,创建图形窗口

plot(w,TH2(:,1),'o-',w,TH2(:,2),'d-',w,TH2(:,3),'s-',w,TH2(:,4),'p-',...

   w,TH2(:,5),'h-',w,TH2(:,6),'<-')      %画振幅曲线族

h=legend([repmat('\it\beta/\omega\rm_0=',n,1),num2str(b')]);%加图例

fs=14;set(h,'FontSize',fs),grid on%字体大小,放大图例,加网格

[m,i]=max(TH2);hold on,stem(w(i),m,'--')%求峰值和下标,保持图像,画杆图

bm=linspace(0.05,sqrt(2)/2);%峰值曲线的约化阻尼因子

omegam=sqrt(1-2*bm.^2);%峰值曲线的驱动力约化角频率

thm=1./(2*bm.*sqrt(1-bm.^2));%峰值曲线的振幅

plot(omegam,thm,'r--','LineWidth',2)%画峰值曲线

xlabel('\it\Omega/\omega\rm_0','FontSize',fs)%标记横坐标

ylabel('\it\Theta\rm_2/\it\theta\rm_0','FontSize',fs),axis([0,wm,0,6])%标记纵坐标,曲线范围

title('受迫振动的幅频响应曲线','FontSize',fs)%标题

PHI2=atan2(-2*B.*W,1-W.^2)*180/pi;figure%初相,创建图形窗口

plot(w,PHI2(:,1),'o-',w,PHI2(:,2),'d-',w,PHI2(:,3),'s-',w,PHI2(:,4),'p-',...

   w,PHI2(:,5),'h-',w,PHI2(:,6),'<-')  %画振幅曲线族

h=legend([repmat('\it\beta/\omega\rm_0=',n,1),num2str(b')]);%加图例

set(h,'FontSize',fs),grid on%放大图例,加网格

thm=sqrt(1-2*bm.^2);thm(imag(thm)~=0)=NaN;%计算分子,复数改为非数

phim=-atan(thm./bm)*180/pi;%位移共振初相

hold on,plot(omegam,phim,'r--','LineWidth',2)%保持图像,画峰值曲线

xlabel('\it\Omega/\omega\rm_0','FontSize',fs)%标记横坐标

ylabel('\it\Phi\rm_2/\circ','FontSize',fs)%标记纵坐标

title('受迫振动的相频响应曲线','FontSize',fs)%标题

V=W./sqrt((1-W.^2).^2+(2*B.*W).^2);figure%约化速度振幅,创建图形窗口

plot(w,V(:,1),'o-',w,V(:,2),'d-',w,V(:,3),'s-',w,V(:,4),'p-',...

   w,V(:,5),'h-',w,V(:,6),'<-')      %画振幅曲线族

h=legend([repmat('\it\beta/\omega\rm_0=',n,1),num2str(b')]);%加图例

set(h,'FontSize',fs),grid on%放大图例,加网格

hold on,plot([1;1],[0;6],'r--','LineWidth',2)%保持图像,画竖线

xlabel('\it\Omega/\omega\rm_0','FontSize',fs)%标记横坐标

ylabel('\it\omega\rm_m/\it\omega\rm_0\it\theta\rm_0','FontSize',fs)%标记纵坐标

title('受迫振动的稳态角速度振幅曲线','FontSize',fs)%标题

axis([0,wm,0,6])%曲线范围

摆轮运动的振幅随时间变化的曲线动画matlab程序

[程序]ltg1main.m

%摆轮在平衡点从静止开始的受迫振动曲线的主程序(用解析式)

clear,tm=80;%清除变量,最大时间

w0=3.572;W=3.566;a=0.8381;b=0.05438;%自由振动圆频率,驱动力圆频率,角加速度,阻尼系数

ltg3fun(w0,W,b,a,tm)%调用函数文件

ltg3fun.m

%摆轮在平衡点从静止开始的受迫振动曲线的函数文件(用解析式)

function fun(w0,W,b,a,tm)

if W==1 W=1-eps;end%如果为1则改小一点

W1=sqrt(w0^2-b^2);%阻尼圆频率

t=linspace(0,tm,20000);%时间向量

TH1=a*sqrt(W1^2*(W^2-w0^2)^2+b^2*(W^2+w0^2)^2)/W1/((W^2-w0^2)^2+4*b^2*W^2);%阻尼振幅

PHI1=atan2(b*(W^2+w0^2),W1*(W^2-w0^2));%阻尼振动初相

th1=TH1*exp(-b*t).*cos(W1*t+PHI1);%阻尼振动角度

TH2=a/sqrt((W^2-w0^2)^2+4*b^2*W^2);%等幅振动振幅

PHI2=atan2(-2*b*W,w0^2-W^2);%等幅振动初相

th2=TH2*cos(W*t+PHI2);%等幅振动角度

figure,plot(t,th1,'--b',t,th2,'-.g','LineWidth',2)%创建图形窗口,画曲线

txt=['\it\beta\rm=',num2str(b)];%阻尼因子字符串

txt=[txt,',\itm\rm=',num2str(a)];%连接角加速度字符串

txt=[txt,',\it\omega\rm_0=',num2str(w0)];%连接自由振动圆频率字符串

txt=[txt,',\it\omega\rm=',num2str(W)];%连接驱动力圆频率字符串

grid on,fs=14;%加网格,字体大小

title(['摆轮振动的角度(',txt,')'],'FontSize',fs)%标题

xlabel('\itt','FontSize',fs)%标记横坐标

hold on,plot([0,tm],[TH2,TH2],':','LineWidth',2)%保持属性,画水平线

text(tm,TH2,num2str(TH2))%标记稳态角度

pause,th=th1+th2;%暂停,振动叠加

hold on,comet(t,th)%保持属性,画彗星线

plot(t,th,'r','LineWidth',2)%补画曲线