数学与非线性科学 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

煜宸0922 发表于 2012-3-2 15:31

胞映射周期确定问题求助

各位前辈，我想问下，在简单胞映射中，如果映射后发现了周期运动，怎么确定前面出现过的这个胞的位置？

水风而动 发表于 2014-3-2 21:23

你好你也在做包映射么？？

gghhjj 发表于 2014-3-5 08:40

可以参考下面这篇文章的程序实现思路

水风而动 发表于 2014-3-23 10:00

gghhjj 发表于 2014-3-5 08:40
可以参考下面这篇文章的程序实现思路

能提供个代码参考下吗？？

gghhjj 发表于 2014-3-25 08:39

本程序为matlab程序，是通过网络收集的，并未经过验证，需要的请自己验证

% 本程序为用胞化积分轨迹法进行全局性态分析

% 初始化，划分胞化空间

clear
xcellnumber=10;
ycellnumber=10;
xleft=-5;
xright=5;
ytop=10;
ybottom=-10;
di=1;
width=(xright-xleft)/xcellnumber;
height=(ytop-ybottom)/ycellnumber;
center(1:xcellnumber/di,1)=[(xleft+width/2):di*width:(xright-width/2)]';
center(1:ycellnumber/di,2)=[(ybottom+height/2):di*height:(ytop-height/2)]';
Id=zeros(xcellnumber*ycellnumber,1);
Ib=zeros(xcellnumber*ycellnumber,1);
Y=zeros(xcellnumber*ycellnumber,2);
dpx=0.001;
dpy=0.001;
nm=10001;
ipsilong=0.0001;
row=0;
periodcells=0;
for i=1:10
for k=1:10
   initvalue((i-1)*10+k,1)=center(i,1);
   initvalue((i-1)*10+k,2)=center(k,2);
end
end
plot(initvalue(:,1),initvalue(:,2),'.');
% omiga=4.7547;
% B=58.4;
% G=311.5;
% eita=0.1;
% K=1;
% U=1e-6;
G=2;
omiga=3.9311;
B=2;
K=1;
U=1;
eita=0.1;

% 从每个胞的中心点开始进行映射。

for i=1:ycellnumber/di
i;
cput=cputime;
for j=1:xcellnumber/di
      j;
   trajectory=zeros(1,2);
   sequence=zeros;
   initial(1,1)=center(j,1);
   initial(1,2)=center(i,2);
   ci=(i-1)*di*xcellnumber+(j-1)*di+1;
   Ib(ci)=Ib(ci)+1;


   if Id(ci)>0
         if abs(Y(ci,1)-initial(1,1))<=dpx&abs(Y(ci,2)-initial(1,2))<=dpy
            continue;
         else
            Id(ci)=-Id(ci);
         end
   else
         Id(ci)=-nm;
         Y(ci,:)=;
   end
   trajectory(1,:)=initial(1,:);
   sequence(1)=ci;
   count1=0;
   count2=1;
   for count=2:2000
         if count==2000
            Id(sequence)=-Id(sequence);
            break
         end

%          =ode45('duffing',,);
%
         opts=simset('Initialstate',);
%          =sim('vanderpol',,opts);
         =sim('duffingscm',,opts);
         =size(t);

            % 当映射的像点不在选定区域

         if x(m,1)>=xright|x(m,1)<xleft|x(m,2)<ybottom|x(m,2)>=ytop
            if count1~=count-1
               k=1;
               count1=0;
            else
               k=k+1;
            end
            if k>12
               for index=1:count2
                     cj=sequence(index);
                     if Id(cj)==-nm
                        Id(cj)=nm+1;
                     else if Id(cj)<0
                           Id(cj)=-Id(cj);
                        end
                     end
               end
               break
            end
            initial=x(m,:);
            count1=count;
            continue;
         end

         count2=count2+1;
         indexx=fix(abs(x(m,1)-xleft)/width)+1;
         indexy=fix(abs(x(m,2)-ybottom)/height)+1;

         ci=(indexy-1)*ycellnumber+indexx;
         sequence(count2)=ci;
         trajectory(count2,:)=;
         Ib(ci)=Ib(ci)+1;


      % 定义四种状态

         if Id(ci)==0
            state=1;
         else
            if Id(ci)==-nm
               state=2;
            else
               if Id(ci)>0
                     state=3;
               else
                     state=4;
               end
            end
         end

      % 对四个状态分别进行处理

         switch state
         case 1
            Id(ci)=-nm;
            Y(ci,:)=;
            initial=x(m,:);
         case 2
            if sqrt(((Y(ci,1)-x(m,1))^2+(Y(ci,2)-x(m,2))^2))<=ipsilong
               period=0;
               for index=count2-1:-1:1
                     if sequence(index)==ci
                        period(1,1:(count2-index))=sequence(index+1:count2);
                        break
                     end
               end
               =size(period);
               isold=0;
               C=find(Id<xcellnumber*ycellnumber);
               Idd=Id(C);
               for index=max(Idd(1,:)):-1:1
                     oldperiod=isoldperiod(period,periodcells);
                     if oldperiod
                        Id(sequence)=index;
                        isold=1;
                        break
                     end
               end
               if isold~=1
                     periodcells(max(Idd)+1,1:c2)=period;
                     Id(sequence)=max(Idd)+1;

                     break
               end
            else
               Y(ci,:)=;
               initial=x(m,:);
            end
         case 3
            if abs(Y(ci,1)-x(m,1))<=dpx&abs(Y(ci,2)-x(m,2))<=dpy
               for index=1:count2
                     cj=sequence(index);
                     if Id(cj)==-nm|Id(cj)==-Id(ci)
                        Id(cj)=Id(ci);
                     else
                        if Id(cj)<0&Id(cj)~=-nm&Id(cj)~=-Id(ci)
                           Id(cj)=-Id(cj);
                        end
                     end
               end
               break
            else
               Id(ci)=-Id(ci);
               initial=x(m,:);
            end
         otherwise
            if abs(Y(ci,1)-x(m,1))<=dpx&abs(Y(ci,2)-x(m,2))<=dpy
               Id(ci)=-Id(ci);
               for index=1:count2
                     cj=sequence(index);
                     if Id(cj)==-nm|Id(cj)==-Id(ci)
                        Id(cj)=Id(ci);
                     else
                        if Id(cj)<0&Id(cj)>-nm&Id(cj)~=-Id(ci)
                           Id(cj)=-Id(cj);
                        end
                     end
               end
               break
            else
               period=0;
               for index=count2-1:-1:1
                     if sequence(index)==ci
                        period(1,1:(count2-index))=sequence(index+1:count2);
                        Y1=trajectory(index,:);
                        break
                     end
               end
               if sqrt(((Y1(1,1)-x(m,1))^2+(Y1(1,2)-x(m,2))^2))<=ipsilong
                     =size(period);
                     isold=0;
                     C=find(Id<xcellnumber*ycellnumber);
                     Idd=Id(C);
                     for index=max(Idd):-1:1
                        oldperiod=isoldperiod(period,periodcells(index,:));
                        if oldperiod
                           Id(sequence)=index;
                           isold=1;
                           break
                        end
                     end
                     if isold~=1
                        periodcells(max(Idd)+1,1:c2)=period;
                        Id(sequence)=max(Idd)+1;
                     end
                     break
               end
            end
            initial=x(m,:);
         end
   end
end
et=(cputime-cput)/60;
end

gghhjj 发表于 2014-3-25 08:41

水风而动发表于 2014-3-23 10:00
能提供个代码参考下吗？？

上面发的是matlab代码，另外论坛还有一个C语言的代码，见帖子
http://forum.vibunion.com/thread-60347-1-1.html

gghhjj 发表于 2014-3-25 08:43

再发一篇胞映射的参考资料

relax 发表于 2014-3-26 17:19

每个胞都有自己的群号，周期，及距离吸引子的步数。

relax 发表于 2014-3-26 17:28

gghhjj 发表于 2014-3-25 08:41
上面发的是matlab代码，另外论坛还有一个C语言的代码，见帖子
http://forum.vibunion.com/thread-60347-1- ...

这个程序缺头文件，要是有就好了。

relax 发表于 2014-3-26 17:30

每个胞的群数周期数距离吸引子的步数都是计算好了的。

lihaitao123 发表于 2014-4-27 11:43

gghhjj 发表于 2014-3-25 08:39
本程序为matlab程序，是通过网络收集的，并未经过验证，需要的请自己验证

161                   oldperiod=isoldperiod(period,periodcells);

162.                      if oldperiod

224.                         oldperiod=isoldperiod(period,periodcells(index,:));
不知道哪位能给解释一下.isoldperiod之前没有出现过. 还有if oldperiod
是不是不完整呢?

gghhjj 发表于 2014-5-16 17:15

lihaitao123 发表于 2014-4-27 11:43
161 oldperiod=isoldperiod(period,periodcells);

162. ...

这个应该是周期胞组检测的子程序，个人收集的程序中没有该函数
如果可能的话，只能你自己根据相关的原理补充了

页: [1]

声振论坛's Archiver

胞映射周期确定问题求助