poincare截面程序求助

zhailiangjun

请各位大虾帮忙看一下我写的这个计算poincare截面图的程序究竟出了什么问题呐。就是画不出图来，自己查了很久了，可还是找不到错误，也不知道该怎样改，小弟先谢谢了。

function dxdt=zhengzefun(t,x)
global w;
global f;
global j;
w=1;
f=0.1;
j=4.5;
dxdt=zeros(4,1);
%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2
dxdt(1)=-w*x(2)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(4)-f*x(2)^2*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dxdt(2)=w*x(1)-x(1)*x(2)*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dxdt(3)=-w*x(4)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(2);
dxdt(4)=w*x(3);

上面是一个哈密顿正则运动方程，我试图画出它的poincare截面图。所以参考本论坛上许多类似的帖子，然后编出了下面的程序。

%poincare
%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2 lamt=f
global w,
global f,
global j
w=1;
f=0.5;
j=4.5;
[t,x]=ode113('zhengzefun',[0:6.3021e-004:630.211],[2.79429;-0.4326;0.1253;-1.6656]);
for i=1:1:(length(x)-1)
    if x(i,4)<0
        if  x(i+1,4)>0;
            if x(i,3)>0
        qq(i)=x(i,1);
        pp(i)=x(i,2);
            end
        end
    end
   if  x(i,4)>0
       if x(i+1,4)<0;
           if x(i,3)>0
         qq(i)=x(i,1);
         pp(i)=x(i,2);
           end
       end
end
if x(i,4)=0
   if x(i,3)>0
   qq(i)=x(i,1);
   pp(i)=x(i,2);
   end
end
end
plot(qq(:),pp(:),'.')

我是想画出x（4）=0的面上，当x（3）>0时的x（1），x（2）的截面图，x（1）与x（2）的数值校正没有写出，但是每次都得出了很少的点。改变了积分步长和积分时间还是不行。是不是我对程序理解的不对呐，请教大家了，多谢多谢。

[ 本帖最后由 eight 于 2008-3-6 16:43 编辑 ]

咕噜噜

%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2
dxdt(1)=-w*x(2)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(4)-f*x(2)^2*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dxdt(2)=w*x(1)-x(1)*x(2)*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dxdt(3)=-w*x(4)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(2);
dxdt(4)=w*x(3);
这里不对，搜索论坛看看别人怎么写的

zhailiangjun

function dx=zhengzefun(t,x)
global w;
global f;
global j;
w=1;
f=0.1;
j=4.5;
%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2
dx=[-w*x(2)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(4)-f*x(2)^2*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
w*x(1)-x(1)*x(2)*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;    -w*x(4)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(2);   w*x(3);]
你看写成这样子就行了吗
或者
function dx=zhengzefun(t,x)
global w;
global f;
global j;
w=1;
f=0.1;
j=4.5;
dx=zeros(4,1);
%x(1)=p1 x(2)=q1 x(3)=p2 x(4)=q2
dx(1)=-w*x(2)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(4)-f*x(2)^2*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dx(2)=w*x(1)-x(1)*x(2)*x(4)*((4*j-x(2)^2-x(1)^2)/4/j)^(-1/2)/j;
dx(3)=-w*x(4)-2*f*((4*j-x(2)^2-x(1)^2)/4/j)^(1/2)*x(2);
dx(4)=w*x(3);
这样对吗

院长姐姐，到底是怎么不对啊，改成上面的样子还是不对啊。:'(  :'(

[ 本帖最后由 ChaChing 于 2010-4-25 16:16 编辑 ]