|
楼主 |
发表于 2017-9-18 20:56
|
显示全部楼层
%强迫振动达芬方程幅频响应曲线的解析求法
%参数设定
epsilon=0.2;%小参数
gammal=4;%非线性参数
F1=0.3;%激励幅值
zetal=0.25;%阻尼项
a=linspace(0.1,2.0,100);%计算响应幅值范围
for ii=1:1:length(a)
%根据(23)式,对应一个响应幅值,可能存在两个激励频率值;
%Omega1为较小的激励频率值
Omega1(ii)=(1+3*epsilon*gammal/8*a(ii)^2-sqrt((epsilon*F1/(2*a(ii)))^2-(epsilon*zetal)^2));
%Omega1对应的本征函数值
lamOmega11=sqrt(-((Omega1(ii)-1)/epsilon-3*gammal/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gammal/8*a(ii)^2))-zetal;
lamOmega12=-sqrt(-((Omega1(ii)-1)/epsilon-3*gammal/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gammal/8*a(ii)^2))-zetal;
%Omega2为较大的激励频率值
Omega2(ii)=(1+3*epsilon*gammal/8*a(ii)^2+sqrt((epsilon*F1/(2*a(ii)))^2-(epsilon*zetal)^2));
%Omega2对应的本征函数值
lamOmega21=sqrt(-((Omega2(ii)-1)/epsilon-3*gammal/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gammal/8*a(ii)^2))-zetal;
lamOmega22=-sqrt(-((Omega2(ii)-1)/epsilon-3*gammal/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gammal/8*a(ii)^2))-zetal;
if lamOmega11==conj(lamOmega21) %如果两个本征值相等退出计算(即此时计算到频响曲线最高点)
plot(Omega1(ii),a(ii),'k.','MarkerSize',15')
hold on
break
else
FG=((Omega1(ii)-1)/epsilon-3*gammal/8*a(ii)^2)*((Omega1(ii)-1)/epsilon-9*gammal/8*a(ii)^2)+zetal^2;
%对应上文中的(31)式
if FG<0
plot(Omega1(ii),a(ii),'color',[0.5 0.5 0.5],'MarkerSize',15) %不稳定点
else
plot(Omega1(ii),a(ii),'k.','MarkerSize',15) %稳定点
hold on
end
GF=((Omega2(ii)-1)/epsilon-3*gammal/8*a(ii)^2)*((Omega2(ii)-1)/epsilon-9*gammal/8*a(ii)^2)+zetal^2;
%对应上下文中的(31)式
if GF<0
plot(Omega2(ii),a(ii),'.','color',[0.5 0.5 0.5],'MarkerSize',15) %不稳定点
else
plot(Omega2(ii),a(ii),'k.','MarkerSize',15)%稳定点
hold on
end
xlabel('\Omega') %x轴标题
ylabel('a') %y轴标题
xlim([0.5,1.8]);%x轴范围
ylim([0,2]);%y轴范围
end %满足条件退出循环(计算得到幅频响应曲线顶点)
end %结束对不同a值得计算
|
评分
-
1
查看全部评分
-
|