向各位大侠求助庞加莱投影图的画法

煜宸0922

各位大侠，小弟刚入非线性不久，做两自由度的间隙碰撞系统，相图已经画出，但庞加莱投影图十分不会，好长时间没进展，请各位大侠指教，不胜感激。

kangarooli

回复 1 # 煜宸0922 的帖子

方便的话你可以把系统发上来看一下，那样可能会更针对些，poincare图我现在参考的是论坛里频闪法思想，你可以搜搜看，好多的，看看有没有适合你的

煜宸0922

回复 2 # kangarooli 的帖子

我要用的是截面法，系统的程序三部分：
1，系统方程的程序：
function dxdt=fun1(t,x)
global w
q1=0;
q2=w;   
q3=0.05;
q4=0;
q5=4;
q6=4;
q7=3;
dxdt=zeros(4,1);
dxdt(1)=x(3);
dxdt(2)=x(4);
dxdt(3)=(1-q1)*sin(q2*t+q4)-2*q3*x(3)+2*q3*x(4)-x(1)+x(2);
dxdt(4)=(q1*sin(q2*t+q4)+2*q3*x(3)-2*q3*(1+q5)*x(4)+x(1)-(1+q6)*x(2))/q7;
end
2.龙格库塔法：
function  x_new=rk144(x,h,T)
k1=h*feval('fun1',T,x(:,1));
k2=h*feval('fun1',T+h/2,x(:,1)+k1/2);
k3=h*feval('fun1',T+h/2,x(:,1)+k2/2);
k4=h*feval('fun1',T+h,x(:,1)+k3);
x_new=x(:,1)+(k1+2*k2+2*k3+k4)/6;
end
3.主程序
function [T,x,V]=changerk344(N)
x(:,1)=[0;0;0;0;];      % 初始状态
h=0.01;                % 初始步长,最大步长
AbsTol=1e-5;         % 容差
p=4;                        
T=zeros(1,N);          % 时间向量
T(1,1)=0;              % 时间起点
p21=2^p-1;
V =[];
tic
for i=1:N-1
     if abs(x(1,i))<0.09             % 定步长
        h=0.01;      
        x(:,i+1)=rk144(x(:,i),h,T(1,i));
        T(1,i+1) = T(1,i) + h;          % 更新时间
    else                                 % 变步长
        x1=rk144(x(:,i),h,T(1,i));
        x2=rk144(x(:,i),h/2,T(1,i));
        a=abs(x1-x2)/p21;
        while a(1,1)>AbsTol                % 选择合适步长
            h=h/10;
            x1=rk144(x(:,i),h,T(1,i));
            x2=rk144(x(:,i),h/2,T(1,i));
            a=abs(x1-x2)/p21;
        end
        x(:,i+1)=x1;      % 根据选择的步长更新状态
       if   abs(x(1,i+1))>=0.1
           V=x(3,i+1);
           x(3,i+1)=-0.8*x(3,i+1);
           x(1,i+1)=sign(x(1,i+1))*min(0.1,abs(x(1,i+1)));  %限制位置在0.1处
       end
    T(1,i+1) = T(1,i) + h;   
      end
     h=0.01;
end
toc
T = T(:);   % 将时间转换为列向量
x = x(:,1:N);
选的截面是x(1)=0.1,1的速度等于碰后的速度

kangarooli

回复 3 # 煜宸0922 的帖子

截面法没研究过，不过你程序既然有了，那现在问题是什么呢，出来的映射图对应不上还是怎么的

煜宸0922

回复 4 # kangarooli 的帖子

程序只能画出来相图，庞加莱截面图我不会画呀

kangarooli

回复 5 # 煜宸0922 的帖子

我运行了下你的程序老是有错误，看到截面我还以为你上面的程序是画截面图的呢，这么说你是根据荣格库塔思想自己编程序画的相图啊，我看你的模型也不是太复杂，怎么不用matlab自带的ode45什么的解呢，可能我还没看明白你的问题

煜宸0922

回复 6 # kangarooli 的帖子

程序是对的，我这边运行没有问题，那三个程序要分开存储，然后调用。我们要求要用自己编的龙格库塔法。模型就是两自由度含间隙的碰撞振动系统，一个二阶常微分方程组。化成四状态方程。其中x1是质块1的位移，x3是质块1的速度，x2是质块2的位移，x4是它的速度。相图已经可以画出，但我要继续画庞加莱投影图和分岔图，但我试了试都不对，所以才求助

煜宸0922

回复 7 # 煜宸0922 的帖子

主程序中用的变步长，另加碰撞条件

kangarooli

回复 7 # 煜宸0922 的帖子

程序我是分开存储的，但老是提示N未定义，可能是我这里的问题，模型基本都是一个类型的，惭愧啊，以前都用现成的，还没做过变步长的呢，期待高手给你解答

煜宸0922

回复 9 # kangarooli 的帖子

n是点数，可以直接赋值的。

liliangbiao

频闪法不适合做含有间隙的碰撞动力系统，需要在碰撞面上建立poincare section. 建议看看罗冠炜的书籍和论文。

kezairenjian

liliangbiao 发表于 2011-6-9 15:44

                               
登录/注册后可看大图

频闪法不适合做含有间隙的碰撞动力系统，需要在碰撞面上建立poincare section. 建议看看罗冠炜的书籍和论文 ...

关于间隙振动系统的庞加莱截面法，是否可以不建立poincare setion，而直接在相图上找到截面呢？

liliangbiao

建议看看罗冠炜的书籍和论文,我只能告诉你这么多。呵呵。至于poincare section 那是必须的！因为碰撞面就是速度相反的时刻。那么这个时刻就是所有轨线都被Poincare垂直相交的平面，难道你的相图上看不出来吗？

meiyongyuandeze

可以截取碰撞面来最为庞加莱面吧，当初谢建华老师就是这么给我们讲的！

liliangbiao

呵呵，你还上过谢教授的课，那么你应该会了，不妨给他说说！