对于duffing方程 幅频响应 这一经典问题的求助

zzqmessi

duffing方程 应该是很经典的问题了 试了很多方法 都不能很好得到其幅频响应 特指跳跃现象


第一
用ezplot 直接或者扭曲的骨干曲线


第二
用 谐波平衡法求解  当然用平均法 多尺度 都可以得到其关系式
但是再求解幅值时 利用的是盛金公式（见百度百科）
求出来的解跟实际情况有很大的差距
具体程序如下：
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);
w=0.7:0.01:1.4;
figure;plot(w,A1S,'+',w,A2S,'.',w,A3S,'*');


请高手指出其中的错误 谢谢

lihaitao123

我记得以前直接matlab中的“solve”来求解三次方程，只保留正实根就可以了！ 还是觉得盛金公式太书面化。

zzqmessi

lihaitao123 发表于 2014-4-18 15:36
我记得以前直接matlab中的“solve”来求解三次方程，只保留正实根就可以了！ 还是觉得盛金公式太书面化。

嗯 谢谢你 只留实根 选择自己需要的点

dum

可以用伽辽金方法吗，我觉得这个方法很简单啊，为啥没人用

zzqmessi

dum 发表于 2014-5-12 08:28
可以用伽辽金方法吗，我觉得这个方法很简单啊，为啥没人用

谢谢 用过  最后我明白了 需要自己人为的去选择 实数根