胞映射周期确定问题求助

煜宸0922

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

水风而动

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

gghhjj

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

水风而动

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

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

gghhjj

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

1. % 本程序为用胞化积分轨迹法进行全局性态分析
   
2. 
   
3. % 初始化，划分胞化空间
   
4. 
   
5. clear
   
6. xcellnumber=10;
   
7. ycellnumber=10;
   
8. xleft=-5;
   
9. xright=5;
   
10. ytop=10;
    
11. ybottom=-10;
    
12. di=1;
    
13. width=(xright-xleft)/xcellnumber;
    
14. height=(ytop-ybottom)/ycellnumber;
    
15. center(1:xcellnumber/di,1)=[(xleft+width/2):di*width:(xright-width/2)]';
    
16. center(1:ycellnumber/di,2)=[(ybottom+height/2):di*height:(ytop-height/2)]';
    
17. Id=zeros(xcellnumber*ycellnumber,1);
    
18. Ib=zeros(xcellnumber*ycellnumber,1);
    
19. Y=zeros(xcellnumber*ycellnumber,2);
    
20. dpx=0.001;
    
21. dpy=0.001;
    
22. nm=10001;
    
23. ipsilong=0.0001;
    
24. row=0;
    
25. periodcells=0;
    
26. for i=1:10
    
27.     for k=1:10
    
28.         initvalue((i-1)*10+k,1)=center(i,1);
    
29.         initvalue((i-1)*10+k,2)=center(k,2);
    
30.     end
    
31. end
    
32. plot(initvalue(:,1),initvalue(:,2),'.');
    
33. % omiga=4.7547;
    
34. % B=58.4;
    
35. % G=311.5;
    
36. % eita=0.1;
    
37. % K=1;
    
38. % U=1e-6;
    
39. G=2;
    
40. omiga=3.9311;
    
41. B=2;
    
42. K=1;
    
43. U=1;
    
44. eita=0.1;
    
45. 
    
46. % 从每个胞的中心点开始进行映射。
    
47. 
    
48. for i=1:ycellnumber/di
    
49.    i;
    
50.    cput=cputime;
    
51.     for j=1:xcellnumber/di
    
52.          j;
    
53.         trajectory=zeros(1,2);            
    
54.         sequence=zeros;
    
55.         initial(1,1)=center(j,1);
    
56.         initial(1,2)=center(i,2);
    
57.         ci=(i-1)*di*xcellnumber+(j-1)*di+1;
    
58.         Ib(ci)=Ib(ci)+1;
    
59. 
    
60.       
    
61.         if Id(ci)>0
    
62.             if abs(Y(ci,1)-initial(1,1))<=dpx&abs(Y(ci,2)-initial(1,2))<=dpy
    
63.                 continue;
    
64.             else
    
65.                 Id(ci)=-Id(ci);
    
66.             end
    
67.         else               
    
68.             Id(ci)=-nm;
    
69.             Y(ci,:)=[initial(1,1),initial(1,2)];
    
70.         end
    
71.         trajectory(1,:)=initial(1,:);
    
72.         sequence(1)=ci;
    
73.         count1=0;
    
74.         count2=1;
    
75.         for count=2:2000
    
76.             if count==2000
    
77.                 Id(sequence)=-Id(sequence);
    
78.                 break
    
79.             end
    
80.                
    
81. %             [t,x]=ode45('duffing',[0 2*pi/omiga],[initial(1,1),initial(1,2)]);
    
82. %            
    
83.             opts=simset('Initialstate',[initial(1,1),initial(1,2)]);
    
84. %             [t,x,y]=sim('vanderpol',[2],opts);
    
85.             [t,x,y]=sim('duffingscm',[2*pi/omiga],opts);
    
86.             [m,e]=size(t);
    
87.             
    
88.                % 当映射的像点不在选定区域
    
89. 
    
90.             if x(m,1)>=xright|x(m,1)<xleft|x(m,2)<ybottom|x(m,2)>=ytop      
    
91.                 if count1~=count-1
    
92.                     k=1;
    
93.                     count1=0;
    
94.                 else
    
95.                     k=k+1;
    
96.                 end                              
    
97.                 if k>12
    
98.                     for index=1:count2
    
99.                         cj=sequence(index);
    
100.                         if Id(cj)==-nm
     
101.                             Id(cj)=nm+1;
     
102.                         else if Id(cj)<0
     
103.                                 Id(cj)=-Id(cj);
     
104.                             end
     
105.                         end
     
106.                     end
     
107.                    break
     
108.                 end
     
109.                 initial=x(m,:);
     
110.                 count1=count;
     
111.                 continue;
     
112.             end
     
113.             
     
114.             count2=count2+1;
     
115.             indexx=fix(abs(x(m,1)-xleft)/width)+1;
     
116.             indexy=fix(abs(x(m,2)-ybottom)/height)+1;
     
117.             
     
118.             ci=(indexy-1)*ycellnumber+indexx;
     
119.             sequence(count2)=ci;
     
120.             trajectory(count2,:)=[x(m,1),x(m,2)];
     
121.             Ib(ci)=Ib(ci)+1;
     
122.             
     
123.             
     
124.          % 定义四种状态
     
125.          
     
126.             if Id(ci)==0
     
127.                 state=1;
     
128.             else
     
129.                 if Id(ci)==-nm
     
130.                     state=2;
     
131.                 else
     
132.                     if Id(ci)>0
     
133.                         state=3;
     
134.                     else
     
135.                         state=4;
     
136.                     end
     
137.                 end
     
138.             end
     
139.                     
     
140.          % 对四个状态分别进行处理
     
141.             
     
142.             switch state
     
143.             case 1
     
144.                 Id(ci)=-nm;
     
145.                 Y(ci,:)=[x(m,1),x(m,2)];
     
146.                 initial=x(m,:);
     
147.             case 2
     
148.                 if sqrt(((Y(ci,1)-x(m,1))^2+(Y(ci,2)-x(m,2))^2))<=ipsilong
     
149.                     period=0;                 
     
150.                    for index=count2-1:-1:1
     
151.                         if sequence(index)==ci
     
152.                             period(1,1:(count2-index))=sequence(index+1:count2);
     
153.                             break
     
154.                         end
     
155.                     end
     
156.                     [r2,c2]=size(period);
     
157.                     isold=0;
     
158.                     C=find(Id<xcellnumber*ycellnumber);
     
159.                     Idd=Id(C);
     
160.                     for index=max(Idd(1,:)):-1:1
     
161.                         oldperiod=isoldperiod(period,periodcells);
     
162.                         if oldperiod
     
163.                             Id(sequence)=index;  
     
164.                             isold=1;
     
165.                             break
     
166.                         end
     
167.                     end
     
168.                     if isold~=1
     
169.                         periodcells(max(Idd)+1,1:c2)=period;
     
170.                         Id(sequence)=max(Idd)+1;
     
171.                         
     
172.                         break
     
173.                     end
     
174.                 else
     
175.                     Y(ci,:)=[x(m,1),x(m,2)];
     
176.                     initial=x(m,:);
     
177.                 end
     
178.             case 3
     
179.                 if abs(Y(ci,1)-x(m,1))<=dpx&abs(Y(ci,2)-x(m,2))<=dpy
     
180.                     for index=1:count2
     
181.                         cj=sequence(index);
     
182.                         if Id(cj)==-nm|Id(cj)==-Id(ci)
     
183.                             Id(cj)=Id(ci);
     
184.                         else
     
185.                             if Id(cj)<0&Id(cj)~=-nm&Id(cj)~=-Id(ci)
     
186.                                 Id(cj)=-Id(cj);
     
187.                             end
     
188.                         end
     
189.                     end
     
190.                     break
     
191.                 else
     
192.                     Id(ci)=-Id(ci);
     
193.                     initial=x(m,:);
     
194.                 end
     
195.             otherwise
     
196.                 if abs(Y(ci,1)-x(m,1))<=dpx&abs(Y(ci,2)-x(m,2))<=dpy
     
197.                     Id(ci)=-Id(ci);
     
198.                     for index=1:count2
     
199.                         cj=sequence(index);
     
200.                         if Id(cj)==-nm|Id(cj)==-Id(ci)
     
201.                             Id(cj)=Id(ci);
     
202.                         else
     
203.                             if Id(cj)<0&Id(cj)>-nm&Id(cj)~=-Id(ci)
     
204.                                 Id(cj)=-Id(cj);
     
205.                             end
     
206.                         end
     
207.                     end
     
208.                     break
     
209.                 else
     
210.                     period=0;
     
211.                     for index=count2-1:-1:1
     
212.                         if sequence(index)==ci
     
213.                             period(1,1:(count2-index))=sequence(index+1:count2);
     
214.                             Y1=trajectory(index,:);
     
215.                             break
     
216.                         end
     
217.                     end
     
218.                     if sqrt(((Y1(1,1)-x(m,1))^2+(Y1(1,2)-x(m,2))^2))<=ipsilong
     
219.                         [r2,c2]=size(period);                  
     
220.                         isold=0;
     
221.                         C=find(Id<xcellnumber*ycellnumber);
     
222.                         Idd=Id(C);
     
223.                         for index=max(Idd):-1:1
     
224.                             oldperiod=isoldperiod(period,periodcells(index,:));
     
225.                             if oldperiod
     
226.                                 Id(sequence)=index;  
     
227.                                 isold=1;
     
228.                                 break
     
229.                             end
     
230.                         end
     
231.                         if isold~=1
     
232.                             periodcells(max(Idd)+1,1:c2)=period;
     
233.                             Id(sequence)=max(Idd)+1;                                                   
     
234.                         end
     
235.                         break
     
236.                     end
     
237.                 end
     
238.                 initial=x(m,:);
     
239.             end
     
240.         end
     
241.     end
     
242.     et=(cputime-cput)/60;
     
243. end

复制代码

gghhjj

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

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

gghhjj

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

relax

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

relax

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

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

relax

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

lihaitao123

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

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

162.                       ...

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