betaa=0.25;
F=1.093;
v=2/3;
Poin=inline(['[x(2);',...
'-2*betaa*x(2)-x(2).^2.*sin(x(1))+F*cos(v*t);',...
'v]'],...
't','x','flag','betaa','F','v');
% Poincare_section[绘制庞加莱截面图]
=ode45(Poin,,,[],betaa,F,v);
x(:,2)=mod(x(:,2),2*pi)-pi;
phi0=pi*2/3; % 选择phi=2*pi/3这个截面
for k=1:round(max(x(:,3))/2/pi);
d=x(:,3)-(k-1)*2*pi-phi0;
=sort(abs(d));
x1l=x(K(1),1);
x1r=x(K(2),1);
x2l=x(K(1),2);
x2r=x(K(2),2);
x3l=x(K(1),3);
x3r=x(K(2),3);
if abs(P(1))+abs(P(2))<3e-16;
X1(k)=x1l;
X2(k)=x2l;
else
Q=polyfit(,,1);
X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit(,,1);
X2(k)=polyval(Q,(k-1)*2*pi-phi0);
end
end
plot(X1,X2,'.');
xlabel('\theta','fontsize',14);
ylabel('d\theta/dt','fontsize',14);
我怎么用这个程序画的图点数都是一样多的。线性的也是三个点。不对了啊。
多自由度的话，是不是将x1,dx1/dt,t 三个作为一组数据来画啊？
x2,dx2/dt,t再作为一组。
这样的话，上面的x(:,3)就为t还是1呢？ 这个程序有这么麻烦吗？
rigid.m
function dx = rigid(t,x)
dx = zeros(2,1); % a column vector
betaa=0.25;
F=1.093;
v=2/3;
dx(1)= x(2);
dx(2)= -2*betaa*x(2)-x(2)^2*sin(x(1))+F*cos(v*t);
options = odeset('RelTol',1e-6,'AbsTol',);
= ode45(@rigid,,,options);
plot(Y(500:end,1),Y(500:end,2),'.');
axis(); 是看的论坛上的一个程序,也是觉得太麻烦.没弄清楚。 2种绘制poincare方法不一样。
siyanger在18楼贴出来的程序是“傻傻”的取一个截面，这是按照书上一般介绍的定义做的，因为不太可能刚好算出截面上的点，所以只好用拟合、插值的方法来得到截面上点的值；
19楼的gghhjj的方法直接算出了phi=2*pi/3面上点的值，不用拟合和插值，当然简单。
推荐用第2种方法。 原帖由 toes 于 2006-9-14 10:45 发表
2种绘制poincare方法不一样。
siyanger在18楼贴出来的程序是“傻傻”的取一个截面，这是按照书上一般介绍的定义做的，因为不太可能刚好算出截面上的点，所以只好用拟合、插值的方法来得到截面上点的值；
19楼的 ...
两种方法都存在很大的问题，不过没办法 原帖由 gghhjj 于 2006-9-15 01:08 发表
两种方法都存在很大的问题，不过没办法
有什么问题，说说撒。
询问关于poincare截面图的具体情况
看了半天上面关于poincare截面图的，还是知道得不具体，例如对于某个非线性方程，已知初始条件和参数，如何画？编程中各符号对应方程的那个参数？还请明示？ 谁能把2楼的程序详细地解释一下啊 原帖由 sczhang 于 2006-8-14 15:03 发表" phi0=pi*2/3; % 选择phi=2*pi/3这个截面 "请问这个截面是根据什么选择的？？？
一般是任意选取的
Matlab中怎样去掉坐标，而保留坐标轴
回复 #7 toes 的帖子
“频闪法”在哪本书上可以看到，还有关于“多频激励”有那些专著？ 哈哈，正要学习poincare方法 ，在这找到方向了谢谢了
新手
中华的程序好象画不出来截面图哦，大家有没有有还有没有谐振子的例子啊 18楼的步长取太长了吧？对计算精度有影响吧？[ 本帖最后由 tigerlun 于 2008-3-20 18:15 编辑 ]