求解这种图像MATLAB如何画出来的！！！

PrettyJerry

系统的方程为：


系统的图像为：


图中F的大小分别为：0.002,0.01,0.05,0.1,0.2
求大神解救，真情感谢！！！

linlin820

这个是典型的非线性系统跳跃现象，可以参考下面的文献

PrettyJerry

linlin820 发表于 2017-9-13 15:45
这个是典型的非线性系统跳跃现象，可以参考下面的文献

多谢！

PrettyJerry

linlin820 发表于 2017-9-13 15:45
这个是典型的非线性系统跳跃现象，可以参考下面的文献

你好，文献里面2.1程序我运行出来是这种。和结果差太多。这是什么情况？？？

MVH

PrettyJerry 发表于 2017-9-17 18:28
你好，文献里面2.1程序我运行出来是这种。和结果差太多。这是什么情况？？？

是不是绘图所选取的参数域不全造成的
这种图应该是分段做出来的

PrettyJerry

MVH 发表于 2017-9-18 11:17
是不是绘图所选取的参数域不全造成的
这种图应该是分段做出来的

谢谢帮助，是我编程有问题。修改后已经能画出来了。

陌影

可否把代码分享一下 楼主

PrettyJerry

陌影 发表于 2017-9-18 16:05
可否把代码分享一下 楼主

%强迫振动达芬方程幅频响应曲线的解析求法
%参数设定
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值得计算

sizhiyuan2006

终于解决了，谢谢分享