时滞系统的最大lyapunov指数估算？

朝阳似火

第二，如果不一样，那一般应用什么方法去求时滞系统的最大lyapunov指数？

springsigs

两者没有关系，你给出的是运动微分方程，一般的处理思路是现对微分方程进行求解，求的方程的时间序列，然后再用你所说的wolf法等方法对相关序列进行分析求解lyapunov指数

朝阳似火

你有什么好的带有程序的参考资料吗？

朝阳似火

我把写的程序发上来，但是有问题，请问怎么修改？

1. function dx=loren(t,x,z)
   
2. %global r k4 b d1 a A B u mu1 C d2 tau
   
3. xlag=z(:,1);
   
4. s=zeros(4,1);
   
5. %% 模型形式
   
6. s(1)=50+1000.*x(2).*x(4)-0.000098.*x(3).*x(1)-0.011.*x(1);
   
7. s(2)=0.000098.*x(3).*x(1)-0.1.*x(2)-1200.*x(2).*x(4);
   
8. s(3)=200.*x(2)-0.67.*x(3);
   
9. s(4)=exp(-0.5).*0.05.*xlag(2).*xlag(4)-0.3.*x(4);

复制代码

1. 
   
2. function lambda_1=lyapunov_wolf(data,N,m,tau,p)
   
3. %  该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法
   
4. %  m: 嵌入维数
   
5. %  tau:时间延迟
   
6. %  data:时间序列
   
7. %  N:时间序列长度
   
8. %  lambda_1:返回最大lyapunov指数值
   
9. %**************************************************************************
   
10. %ode计算整数阶系统的时间序列
    
11. %******************************************************************
    
12. Delays=[1];
    
13. delt_t1 = 0.001;
    
14. t1 = 0:delt_t1:60;  
    
15. y1=dde23('loren',Delays,[537.35;16.31;12.17;0.01],t1);
    
16. xx1 = y1(:,1)';  
    
17. x1 = spline(y1.x, xx1, t1);  
    
18. data= x1(20000:10:60000);                                %采样
    
19. N=length(data);
    
20. m=4;
    
21. tau=3;
    
22. %*****************************************************
    
23. % FFT计算平均周期  
    
24. %**********************************************************
    
25. x=data;  
    
26. xPower=abs(fft(x)).^2;
    
27. NN=length(xPower);  
    
28. xPower(1)=[];                                           %去除直流分量
    
29. NN=floor(NN/2);  
    
30. xPower=xPower(1:NN);
    
31. freq=(1:NN)/NN*0.5;
    
32. [mP,index]=max(xPower);
    
33. P=index;  
    
34. %*************************************************************
    
35. min_point=1  ; %&&要求最少搜索到的点数
    
36. MAX_CISHU=5 ;  %&&最大增加搜索范围次数
    
37. %FLYINGHAWK
    
38. %   求最大、最小和平均相点距离
    
39.     max_d = 0;                                         %最大相点距离
    
40.     min_d = 1.0e+100;                                  %最小相点距离
    
41.     avg_dd = 0;
    
42.     Y=reconstitution(data,N,m,tau);                    %相空间重构
    
43.     M=N-(m-1)*tau;                                     %重构相空间中相点的个数
    
44.     for i = 1 : (M-1)
    
45.         for j = i+1 : M
    
46.             d = 0;
    
47.             for k = 1 : m
    
48.                 d = d + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
    
49.             end
    
50.             d = sqrt(d);
    
51.             if max_d < d
    
52.                max_d = d;
    
53.             end
    
54.             if min_d > d
    
55.                min_d = d;
    
56.             end
    
57.             avg_dd = avg_dd + d;
    
58.         end
    
59.     end
    
60.     avg_d = 2*avg_dd/(M*(M-1));                %平均相点距离
    
61.      
    
62.     dlt_eps = (avg_d - min_d) * 0.02 ;         %若在min_eps～max_eps中找不到演化相点时，对max_eps的放宽幅度
    
63.     min_eps = min_d + dlt_eps / 2 ;            %演化相点与当前相点距离的最小限
    
64.     max_eps = min_d + 2 * dlt_eps  ;           %&&演化相点与当前相点距离的最大限
    
65.      
    
66. %     从P+1～M-1个相点中找与第一个相点最近的相点位置(Loc_DK)及其最短距离DK
    
67.     DK = 1.0e+100;                             %第i个相点到其最近距离点的距离
    
68.     Loc_DK = 2;                                %第i个相点对应的最近距离点的下标
    
69.     for i = (P+1):(M-1)                        %限制短暂分离，从点P+1开始搜索
    
70.         d = 0;
    
71.         for k = 1 : m
    
72.             d = d + (Y(k,i)-Y(k,1))*(Y(k,i)-Y(k,1));
    
73.         end
    
74.         d = sqrt(d);
    
75.         if (d < DK) & (d > min_eps)  
    
76.            DK = d;
    
77.            Loc_DK = i;
    
78.         end
    
79.     end
    
80. %     以下计算各相点对应的李氏数保存到lmd()数组中
    
81. %     i 为相点序号，从1到(M-1)，也是i-1点的演化点；Loc_DK为相点i-1对应最短距离的相点位置，DK为其对应的最短距离
    
82. %     Loc_DK+1为Loc_DK的演化点，DK1为i点到Loc_DK+1点的距离，称为演化距离
    
83. %     前i个log2（DK1/DK）的累计和用于求i点的lambda值
    
84.     sum_lmd = 0 ;                              % 存放前i个log2（DK1/DK）的累计和
    
85.     for i = 2 : (M-1)                          % 计算演化距离      
    
86.         DK1 = 0;
    
87.         for k = 1 : m
    
88.             DK1 = DK1 + (Y(k,i)-Y(k,Loc_DK+1))*(Y(k,i)-Y(k,Loc_DK+1));
    
89.         end
    
90.         DK1 = sqrt(DK1);
    
91.         old_Loc_DK = Loc_DK ;                  % 保存原最近位置相点
    
92.         old_DK=DK;
    
93. 
    
94. %     计算前i个log2（DK1/DK）的累计和以及保存i点的李氏指数
    
95.         if (DK1 ~= 0)&( DK ~= 0)
    
96.            sum_lmd = sum_lmd + log(DK1/DK) /log(2);
    
97.         end
    
98.         lmd(i-1) = sum_lmd/(i-1);
    
99. %     以下寻找i点的最短距离：要求距离在指定距离范围内尽量短，与DK1的角度最小
    
100.         point_num = 0  ; % &&在指定距离范围内找到的候选相点的个数
     
101.         cos_sita = 0  ; %&&夹角余弦的比较初值 ——要求一定是锐角
     
102.         zjfwcs=0     ;%&&增加范围次数
     
103.          while (point_num == 0)
     
104.            % * 搜索相点
     
105.             for j = 1 : (M-1)
     
106.                 if abs(j-i) <=(P-1)      %&&候选点距当前点太近，跳过！
     
107.                    continue;      
     
108.                 end
     
109.                  
     
110.                 %*计算候选点与当前点的距离
     
111.                 dnew = 0;
     
112.                 for k = 1 : m
     
113.                    dnew = dnew + (Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));
     
114.                 end
     
115.                 dnew = sqrt(dnew);
     
116.                  
     
117.                 if (dnew < min_eps)|( dnew > max_eps )   %&&不在距离范围，跳过！
     
118.                   continue;              
     
119.                 end
     
120.                                 
     
121.                 %*计算夹角余弦及比较
     
122.                 DOT = 0;
     
123.                 for k = 1 : m
     
124.                     DOT = DOT+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,old_Loc_DK+1));
     
125.                 end
     
126.                 CTH = DOT/(dnew*DK1);
     
127.                  
     
128.                 if acos(CTH) > (3.14151926/4)      %&&不是小于45度的角，跳过！
     
129.                   continue;
     
130.                 end
     
131.                  
     
132.                 if CTH > cos_sita   %&&新夹角小于过去已找到的相点的夹角，保留
     
133.                     cos_sita = CTH;
     
134.                     Loc_DK = j;
     
135.                     DK = dnew;
     
136.                 end
     
137. 
     
138.                 point_num = point_num +1;
     
139.                  
     
140.             end         
     
141.          
     
142.             if point_num <= min_point
     
143.                max_eps = max_eps + dlt_eps;
     
144.                zjfwcs =zjfwcs +1;
     
145.                if zjfwcs > MAX_CISHU    %&&超过最大放宽次数，改找最近的点
     
146.                    DK = 1.0e+100;
     
147.                    for ii = 1 : (M-1)
     
148.                       if abs(i-ii) <= (P-1)      %&&候选点距当前点太近，跳过！
     
149.                        continue;      
     
150.                       end
     
151.                       d = 0;
     
152.                       for k = 1 : m
     
153.                           d = d + (Y(k,i)-Y(k,ii))*(Y(k,i)-Y(k,ii));
     
154.                       end
     
155.                       d = sqrt(d);
     
156.          
     
157.                       if (d < DK) & (d > min_eps)  
     
158.                          DK = d;
     
159.                          Loc_DK = ii;
     
160.                       end
     
161.                    end
     
162.                    break;  
     
163.                end
     
164.                point_num = 0          ;     %&&扩大距离范围后重新搜索
     
165.                cos_sita = 0;
     
166.             end
     
167.         end
     
168.    end
     
169. 
     
170. %取平均得到最大李雅普诺夫
     
171. lambda_1=sum(lmd)/length(lmd);
     
172. plot(lmd)

复制代码

1. function X=reconstitution(data,N,m,tau)
   
2. %该函数用来重构相空间
   
3. % m      为嵌入空间维数
   
4. % tau    为时间延迟
   
5. % data   为输入时间序列
   
6. % N      为时间序列长度
   
7. % X       为输出,是m*n维矩阵
   
8. M=N-(m-1)*tau;     %相空间中点的个数
   
9. for j=1:M          %相空间重构
   
10.   for i=1:m
    
11.       X(i,j)=data((i-1)*tau+j);
    
12.   end
    
13. end

复制代码

暂时运行结果图：

朝阳似火

结果还是我自己搞出来啦。

朝阳似火

1. wolf计算的最大lyapunov指数值，请问下图的横坐标是表示什么？

2. 下图t1的初值与末值应怎样取最好？这个是否求数值解的时间的初值与末值有关？

3. 为什么下图这个只取y1(:,1)'，而不取y1(:,2)'，y1(:,3)'？

ranjie

可以分享一下吗，