本帖最后由 thdl520 于 2014-11-18 22:21 编辑
求雪缘和无水两位元老 以及IHB方面的专家们给予解惑啊
下面是我编的非对称分段线性IHB法程序 貌似只能得到稳定的解(可以通过改变初始条件来得到),但同一频率下的不稳定解却怎么也得不到该怎么样做啊
clear
clc
% 参数
ks=0.05;F=0.4;
m=1; % order of subharmonic
eps=0.6;ac=-0.0083;
w0=0.98;
% initial parameter
A0=[ 1 1 0 0 0 0 0];
syms t
CS=[1 cos(t) cos(2*t) cos(3*t) sin(t) sin(2*t) sin(3*t)];
CS1=diff(CS,t);CS2=diff(CS1,t);
% 确定解位置
t0=0;t1=2*pi;
for k=1:30
a=[];e=[];
for i=t0:0.1:t1 %限定区间,步长为0.1
[x,fval]=fsolve(@(t)A0(1)+A0(2)*cos(t)+A0(3)*cos(2*t)+A0(4)*cos(3*t)+A0(5)*sin(t)+A0(6)*sin(2*t)+A0(7)*sin(3*t),i,optimset('Display','off')); %求出一个根
flag=0; %检查是否已经存在在根集中
for j=1:length(a)
if abs(x-a(j))<=1e-3
flag=1; %已经存在
break;
end
end
if flag==0 %没有存在
a=[a,x];e=[e,fval];
end
end
% 选取xo和x1之间的解,去掉极值点
for j=length(a):-1:1
if a(j)>t1||a(j)<t0||abs(e(j))>1e-3
a(j)=[];
end
end
% 排序
a=sort(a);
% linear term integration
rl=inline((w0^2*A0*CS2'+2*ks*m*w0*A0*CS1'-F*m^2*cos(m*t))*CS');
Rl=quadv(rl,0,2*pi);
cl=inline(2*ks*m*w0*CS'*CS1+w0^2*CS'*CS2);
Cl=quadv(cl,0,2*pi);
% nonlinear term integration
a=[0 a 2*pi]
n=length(a);
x0=A0*CS';
Rnl=0;
Cnl=0;
rnl1=inline(x0*CS');
rnl2=inline((eps*x0+(1-eps)*ac)*CS');
cnl=inline(1*CS'*CS);
cn2=inline(eps*CS'*CS);
for i=1:n-1
tc=(a(i)+a(i+1))/2;
if (A0(1)+A0(2)*cos(tc)+A0(3)*cos(2*tc)+A0(4)*cos(3*tc)+A0(5)*sin(tc)+A0(6)*sin(2*tc)+A0(7)*sin(3*tc))-ac>0
Rnl=quadv(rnl1,a(i),a(i+1))+Rnl;
Cnl=quadv(cnl,a(i),a(i+1))+Cnl;
else
Rnl=quadv(rnl2,a(i),a(i+1))+Rnl;
Cnl=quadv(cn2,a(i),a(i+1))+Cnl;
end
end
R=Rl+Rnl
C=Cl+Cnl;
dA=-C\R;
A0=A0+dA'
nm=norm(R)
if nm<1e-6
break
end
end
|