马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
有没有做伪不动点追踪求系统多周期解的前辈帮忙看看
clc
clear all
global Num; global eta11;global eta22;global eta33; global eta13; global eta23;global b;
global k11; global k22; global Fb1; global Fb2; global Fm; global weh; global Fah1;global emsl;
Num=3; eta11=0.00;eta22=0.00;eta13=0.0125;eta23=0.0125;eta33=0.05;k11=1.25;k22=1.25;
Fb1=0;Fb2=0;Fah1=0.05;emsl=0.2;b=1;Fm=0.1;
weh=1.498;
options=odeset('RelTol',1e-8,'AbsTol',1e-9);
det_num=100;
h=2*pi/weh;t1=0;t0=0; h1=h/det_num;K_X=[];
X0=zeros(2*Num,1);
% chuzhi;
for i=1:1000
[t,X]=ode45(@three_degree_gear, [t0,t0+h], X0,options);
t0=t0+h; X0=X(end,:)';
if i>900
K_X=[K_X X0];
end
end
t0=0 ;t1=0;
chuzhi;
% XX=budongdian2(t1,X0,h);
XX_list=[];DDK=[];KKK=[];
% X0=XX;
for i=1:5*det_num
t1=t0+i*h1;i
XX=budongdian2(t0,X0,t1);[t,X]=ode45(@three_degree_gear, [t0,t1], XX,options);
DK=norm(X(end,:)'-XX)
X0=XX;
DDK=[DDK DK];XX_list=[XX_list XX];
if mod(i,det_num)==0
KKK=[KKK [XX;DK]];
end
end
[t,X]=ode45(@three_degree_gear, [0,50*h], XX,options);
***************************************************************************
function XX=budongdian2(t,X0,h)
global Num;
options=odeset('RelTol',1e-8,'AbsTol',1e-9);
X_E=[];X_D=[];
% 不动点迭代法流程
K=1;v=1;t1=t;t2=t1;
while v==1
while K<=2*Num+1
[t,X]=ode45(@three_degree_gear, [t1,t1+h], X0,options);
X1=X(end,:)';
% t1=t1+h;
ek=norm((X1-X0)/X1,2);
if ek<1e-9
XX=X1; v=0
break;
else
X_E=[X_E X1];X_D=[X_D [X0;1]];X0=X1;
end
K=K+1;
end
HC=X_E*pinv(X_D);H=HC(:,1:end-1);C=HC(:,end);
E=eye(2*Num);XX=pinv(E-H)*C;v=0;X0=XX;
end
eb=1;kk=1;
while eb>1e-7;
X2=X0;K=1;v=1; t1=t2;
while v==1
while K<=2*Num+1
[t,X]=ode45(@three_degree_gear, [t1,t1+h], X0,options);
X1=X(end,:)';
ek=norm((X1-X0)/X1,2);
if ek<1e-8
XX=X1; v=0
break;
else
X_E=[X_E X1];X_D=[X_D [X0;1]];X0=X1;
end
K=K+1;
end
HC=X_E*pinv(X_D);H=HC(:,1:end-1);C=HC(:,end);
E=eye(2*Num);XX=pinv(E-H)*C; v=0;
end
eb=norm((X2-XX)/XX,2);kk=kk+1;
if kk>200
kk
break
end
X0=XX;
end
做不出正确结果,为什么? |