求助 duffing 方程 频域分析时 跳跃的程序

zzqmessi

现在只是用ezplot画出图形，但是得不到跳跃的情况 ，请大家指教。 谢谢

zzqmessi

还有一个问题就是 我怎么获得ezplot 里面所有的值呢

zzqmessi

我用谐波平衡法 理论求解 和实际差距太大

zzqmessi

附上程序
clc
clear all
n=0;
m=1;k=1;kn=0.0005;C=0.02;F=1;
for w=0.7:0.01:1.4
    n=n+1;
a=9*kn^2/16;b=3*kn*(k-m*w^2)/2;c=(k-m*w^2)^2+C^2*w^2;d=-F^2;
X=b^2-3*a*c;
Y=b*c-9*a*d;
Z=c^2-3*b*d;
delt=Y^2-4*X*Z;
delt1(n)=Y^2-4*X*Z;
if X==0&&Y==0
    A1(n)=-b/(3*a);
    A2(n)=-c/(b);
    A3(n)=-3*d/(c);
else if delt>0
        W1=X*b+3*a*(-Y+abs(sqrt(Y^2-4*X*Z)))/2;
        W2=X*b+3*a*(-Y-abs(sqrt(Y^2-4*X*Z)))/2;
        %W1=X*b+3*a*(-Y+sqrt(Y^2-4*X*Z))/2;
       % W2=X*b+3*a*(-Y-sqrt(Y^2-4*X*Z))/2;
        A1(n)=(-b-(W1^(1/3)+W2^(1/3)))/(3*a);
        A2(n)=(-b+1/2*(W1^(1/3)+W2^(1/3))+sqrt(3)/2*(W1^(1/3)-W2^(1/3))*i)/(3*a);
        A3(n)=(-b+1/2*(W1^(1/3)+W2^(1/3))-sqrt(3)/2*(W1^(1/3)-W2^(1/3))*i)/(3*a);
    else if delt==0
            K=Y/X;
            A1(n)=-b/a+K;
            A2(n)=-K/2;
            A3(n)=-K/2;
        else if delt<0
                T=(2*X*b-3*a*Y)/(2*sqrt(X^3));
                sita=acos(T);
                A1(n)=(-b-2*sqrt(X)*cos(sita/3))/(3*a);
                A2(n)=(-b+sqrt(X)*(cos(sita/3)+sqrt(3)*sin(sita/3)))/(3*a);
                A3(n)=(-b+sqrt(X)*(cos(sita/3)-sqrt(3)*sin(sita/3)))/(3*a);
            end
        end
    end
end
end
A1=sqrt(A1);
A2=sqrt(A2);
A3=sqrt(A3);
A1S=abs(A1);
A2S=abs(A2);
A3S=abs(A3);

w=0.7:0.01:1.4;
figure;plot(w,A1S,'+',w,A2S,'.',w,A3S,'*');