马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function dbf3
m=20;h=0.3142;t0=0.0;ep=0.00000001;
x(1)=0.06594;x(2)=0.01;a=1;b=1;
c(1)=h/2;c(2)=c(1);c(5)=c(1);
c(3)=h;c(4)=h;s(1)=x(1);s(2)=x(2);
x1(1)=1.0;x1(2)=0.0;
x2(1)=0.0;x2(2)=1.0;
v=1;
while v==1
t=t0;
x(1)=s(1);x(2)=s(2);
iii=iii+1;
for i=1:m
t1=t;ts=t;
for ii=1:2
p(ii)=x(ii);w(ii)=x(ii);
end
for jj=1:4
f(1)=p(2);
f(2)=-0.2*p(2)-4.0*p(1)-p(1).^3+0.03*cos(t);
t=ts+c(jj);
for ii=1:2
p(ii)=c(jj).*f(ii)+w(ii);
x(ii)=c(jj+1).*f(ii)/3.0+x(ii);
end
end
t=t1;ts=t;
for ii=1:2
p(ii)=x1(ii);
w(ii)=x1(ii);
end
for jj=1:4
f(1)=p(2);
f(2)=-0.2*p(2)-(4.0+3.0*x(1).^2).*p(1);
t=ts+c(jj);
for ii=1:2
p(ii)=c(jj).*f(ii)+w(ii);
x1(ii)=c(jj+1).*f(ii)/3.0+x1(ii);
end
end
t=t1;ts=t;
for ii=1:2
p(ii)=x2(ii);
w(ii)=x2(ii);
end
for jj=1:4
f(1)=p(2);
f(2)=-0.2*p(2)-(4.0+3.0*x(1).^2).*p(1);
t=ts+c(jj);
for ii=1:2
p(ii)=c(jj).*f(ii)+w(ii);
x2(ii)=c(jj+1).*f(ii)/3.0+x2(ii);
end
end
b=1;N(a,b)=t;
b=b+1;N(a,b)=x(1);
b=b+1;N(a,b)=x(2);
a=a+1;
end
r(1)=s(1)-x(1);
r(2)=s(2)-x(2);
dr(1,1)=1.0-x1(1);
dr(1,2)=x2(1);
dr(2,1)=x1(2);
dr(2,2)=1.0-x2(2);
I=[1 0;0 1];
P=I-dr;
D=eig(P)
d=dr(1,1).*dr(2,2)-dr(1,2).*dr(2,1);
b=-r(2).*dr(1,2)+r(1).*dr(2,2);
e=-r(1).*dr(2,1)+r(2).*dr(1,1);
ds(1)=b/d;
ds(2)=e/d;
if abs(ds(1))<ep&abs(ds(2))<ep
break
else
s(1)=s(1)-ds(1);
s(2)=s(2)-ds(2);
end
end
N
P
D
这是打靶法求达芬系统的1T周期解的程序,我同时加了点同求系统的雅可比矩阵P和其特征值D,得到的结果是
P =
1.0e-005 *
0.132237285277448 0.049603766980269
-0.197319424225645 0.141870719549964
D =
1.0e-005 *
0.137054002413706 + 0.098815919648602i
0.137054002413706 - 0.098815919648602i
发现值都非常小,不知道合不合理,也不知道这么求对不对
另外有点疑惑的是蓝色标注行,dr(1,2)=x2(1);dr(2,1)=x1(2);是不是x2(1),x1(2)前应该加负号
b=-r(2).*dr(1,2)+r(1).*dr(2,2);e=-r(1).*dr(2,1)+r(2).*dr(1,1);是不是应该改为b=r(2).*dr(1,2)-r(1).*dr(2,2);e=r(1).*dr(2,1)-r(2).*dr(1,1);
请各位帮忙看下
|