查看: 38980|回复: 70

[共享资源] matlab编写的Lyapunov指数计算程序汇总

  [复制链接]
发表于 2006-9-2 06:07 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
申明:以下各程序为个人在网络上收集的Lyapunov指数计算程序,未经过验证,不保证程序的正确性和计算结果的正确性,请大家见谅,也欢迎大家探讨!


计算连续方程Lyapunov指数的程序,比较好用的

其中给出了两个例子:
计算Rossler吸引子的Lyapunov指数
计算超混沌Rossler吸引子的Lyapunov指数

http://www.pudn.com/downloads39/sourcecode/math/detail132231.html

点评

赞成: 4.0
赞成: 4
  发表于 2014-3-27 18:52

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2006-9-2 06:08 | 显示全部楼层
基于Volterra滤波器混沌时间序列多步预测

----------------------------------------------
参考文献:
1、张家树.混沌时间序列的Volterra自适应预测.物理学报.2000.03
2、Scott C.Douglas, Teresa H.-Y. Meng, Normalized Data Nonlinearities for LMS Adaptation. IEEE Trans.Sign.Proc. Vol.42 1994

----------------------------------------------
文件说明:
1、original_MultiStepPred_main.m 程序主文件,直接运行此文件即可
2、original_train.m                 训练函数
3、original_test.m                  测试函数
4、LorenzData.dll                  产生Lorenz离散序列
5、normalize_1.m                 归一化
6、PhaSpaRecon.m                 相空间重构
7、PhaSpa2VoltCoef.m                 构造 Volterra 自适应 FIR 滤波器的输入信号矢量 Un
8、TrainTestSample_2.m                 将特征矩阵前 train_num 个为训练样本,其余为测试样本
9、FIR_NLMS.dll                         NLMS自适应算法

http://www.pudn.com/downloads45/sourcecode/math/detail150408.html
 楼主| 发表于 2006-9-2 06:09 | 显示全部楼层
recnstitution重构相空间,在非线性系统分析中有重要的作用

  1. function [Texp,Lexp]=lyapunov(n,tstart,stept,tend,ystart,ioutp);
  2. global DS;
  3. global P;
  4. global calculation_progress first_call;
  5. global driver_window;
  6. global TRJ_bufer Time_bufer bufer_i;

  7. %
  8. %    Lyapunov exponent calcullation for ODE-system.
  9. %
  10. %    The alogrithm employed in this m-file for determining Lyapunov
  11. %    exponents was proposed in
  12. %
  13. %         A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,
  14. %        "Determining Lyapunov Exponents from a Time Series," Physica D,
  15. %        Vol. 16, pp. 285-317, 1985.
  16. %
  17. %    For integrating ODE system can be used any MATLAB ODE-suite methods.
  18. % This function is a part of MATDS program - toolbox for dynamical system investigation
  19. %    See:    http://www.math.rsu.ru/mexmat/kvm/matds/
  20. %
  21. %    Input parameters:
  22. %      n - number of equation
  23. %      rhs_ext_fcn - handle of function with right hand side of extended ODE-system.
  24. %              This function must include RHS of ODE-system coupled with
  25. %              variational equation (n items of linearized systems, see Example).                  
  26. %      fcn_integrator - handle of ODE integrator function, for example: @ode45                  
  27. %      tstart - start values of independent value (time t)
  28. %      stept - step on t-variable for Gram-Schmidt renormalization procedure.
  29. %      tend - finish value of time
  30. %      ystart - start point of trajectory of ODE system.
  31. %      ioutp - step of print to MATLAB main window. ioutp==0 - no print,
  32. %              if ioutp>0 then each ioutp-th point will be print.
  33. %
  34. %    Output parameters:
  35. %      Texp - time values
  36. %      Lexp - Lyapunov exponents to each time value.
  37. %
  38. %    Users have to write their own ODE functions for their specified
  39. %    systems and use handle of this function as rhs_ext_fcn - parameter.      
  40. %
  41. %    Example. Lorenz system:
  42. %               dx/dt = sigma*(y - x)     = f1
  43. %               dy/dt = r*x - y - x*z = f2
  44. %               dz/dt = x*y - b*z     = f3
  45. %
  46. %    The Jacobian of system:
  47. %        | -sigma  sigma  0 |
  48. %    J = |   r-z    -1   -x |
  49. %        |    y      x   -b |
  50. %
  51. %    Then, the variational equation has a form:
  52. %
  53. %    F = J*Y
  54. %    where Y is a square matrix with the same dimension as J.
  55. %    Corresponding m-file:
  56. %        function f=lorenz_ext(t,X)
  57. %         SIGMA = 10; R = 28; BETA = 8/3;
  58. %         x=X(1); y=X(2); z=X(3);
  59. %
  60. %         Y= [X(4), X(7), X(10);
  61. %             X(5), X(8), X(11);
  62. %             X(6), X(9), X(12)];
  63. %         f=zeros(9,1);
  64. %         f(1)=SIGMA*(y-x); f(2)=-x*z+R*x-y; f(3)=x*y-BETA*z;
  65. %
  66. %         Jac=[-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA];
  67. %  
  68. %         f(4:12)=Jac*Y;
  69. %
  70. %    Run Lyapunov exponent calculation:
  71. %     
  72. %    [T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[0 1 0],10);   
  73. %   
  74. %    See files: lorenz_ext, run_lyap.   
  75. %  
  76. % --------------------------------------------------------------------
  77. % Copyright (C) 2004, Govorukhin V.N.
  78. % This file is intended for use with MATLAB and was produced for MATDS-program
  79. % http://www.math.rsu.ru/mexmat/kvm/matds/
  80. % lyapunov.m is free software. lyapunov.m is distributed in the hope that it
  81. % will be useful, but WITHOUT ANY WARRANTY.
  82. %



  83. %
  84. %       n=number of nonlinear odes
  85. %       n2=n*(n+1)=total number of odes
  86. %

  87. options = odeset('RelTol',DS(1).rel_error,'AbsTol',DS(1).abs_error,'MaxStep',DS(1).max_step,...
  88.                  'OutputFcn',@odeoutp,'Refine',0,'InitialStep',0.001);
  89.             

  90. n_exp = DS(1).n_lyapunov;
  91. n1=n; n2=n1*(n_exp+1);
  92. neq=n2;

  93. %  Number of steps

  94. nit = round((tend-tstart)/stept)+1;

  95. % Memory allocation

  96. y=zeros(n2,1);
  97. cum=zeros(n2,1);
  98. y0=y;
  99. gsc=cum;
  100. znorm=cum;

  101. % Initial values

  102. y(1:n)=ystart(:);

  103. for i=1:n_exp y((n1+1)*i)=1.0; end;

  104. t=tstart;


  105. Fig_Lyap = figure;
  106. set(Fig_Lyap,'Name','Lyapunov exponents','NumberTitle','off');
  107. set(Fig_Lyap,'CloseRequestFcn','');
  108. hold on;
  109. box on;
  110. timeplot = tstart+(tend-tstart)/10;
  111. axis([tstart timeplot -1 1]);
  112. title('Dynamics of Lyapunov exponents');
  113. xlabel('t');
  114. ylabel('Lyapunov exponents');
  115. Fig_Lyap_Axes = findobj(Fig_Lyap,'Type','axes');

  116. for i=1:n_exp
  117.        PlotLyap{i}=plot(tstart,0);        
  118. end;

  119.        uu=findobj(Fig_Lyap,'Type','line');
  120.        for i=1:size(uu,1)
  121.            set(uu(i),'EraseMode','none') ;
  122.            set(uu(i),'XData',[],'YData',[]);
  123.            set(uu(i),'Color',[0 0 rand]);
  124.        end

  125. ITERLYAP = 0;

  126. % Main loop
  127. calculation_progress = 1;

  128. while t<tend
  129.     tt = t + stept;
  130.     ITERLYAP = ITERLYAP+1;
  131.     if tt>tend, tt = tend; end;
  132. % Solutuion of extended ODE system

  133. %  [T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);     
  134.    while calculation_progress == 1
  135.        [T,Y] = integrator(DS(1).method_int,@ode_lin,[t tt],y,options,P,n,neq,n_exp);
  136.        first_call = 0;
  137.        if calculation_progress == 99, break; end;
  138.        if ( T(size(T,1))<tt ) & (calculation_progress~=0)
  139.          y=Y(size(Y,1),:);
  140.          y(1,1:n)=TRJ_bufer(bufer_i,1:n);
  141.          t = Time_bufer(bufer_i);
  142.          calculation_progress = 1;
  143.        else
  144.          calculation_progress = 0;
  145.        end;
  146.    end;

  147.    if (calculation_progress == 99)
  148.        break;
  149.    else
  150.        calculation_progress = 1;
  151.    end;
  152.    
  153.   t=tt;
  154.   y=Y(size(Y,1),:);
  155.   
  156.   first_call = 0;

  157. %
  158. %       construct new orthonormal basis by gram-schmidt
  159. %

  160.   znorm(1)=0.0;
  161.   for j=1:n1 znorm(1)=znorm(1)+y(n1+j)^2; end;

  162.   znorm(1)=sqrt(znorm(1));

  163.   for j=1:n1 y(n1+j)=y(n1+j)/znorm(1); end;

  164.   for j=2:n_exp
  165.       for k=1:(j-1)
  166.           gsc(k)=0.0;
  167.           for l=1:n1 gsc(k)=gsc(k)+y(n1*j+l)*y(n1*k+l); end;
  168.       end;

  169.       for k=1:n1
  170.           for l=1:(j-1)
  171.               y(n1*j+k)=y(n1*j+k)-gsc(l)*y(n1*l+k);
  172.           end;
  173.       end;

  174.       znorm(j)=0.0;
  175.       for k=1:n1 znorm(j)=znorm(j)+y(n1*j+k)^2; end;
  176.       znorm(j)=sqrt(znorm(j));

  177.       for k=1:n1 y(n1*j+k)=y(n1*j+k)/znorm(j); end;
  178.   end;

  179. %
  180. %       update running vector magnitudes
  181. %

  182.   for k=1:n_exp cum(k)=cum(k)+log(znorm(k)); end;

  183. %
  184. %       normalize exponent
  185. %

  186.   rescale = 0;
  187.   u1 =1.e10;
  188.   u2 =-1.e10;

  189.   for k=1:n_exp
  190.       lp(k)=cum(k)/(t-tstart);

  191. %  Plot

  192.       Xd=get(PlotLyap{k},'Xdata');
  193.       Yd=get(PlotLyap{k},'Ydata');

  194.       if timeplot<t
  195.             u1=min(u1,min(Yd));
  196.             u2=max(u2,max(Yd));
  197.       end;
  198.       Xd=[Xd t]; Yd=[Yd lp(k)];
  199.       set(PlotLyap{k},'Xdata',Xd,'Ydata',Yd);
  200.   end;
  201.   if timeplot<t
  202.      timeplot = timeplot+(tend-tstart)/20;
  203.      figure(Fig_Lyap);
  204.      axis([tstart timeplot u1 u2]); end;

  205.   drawnow;



  206. % Output modification

  207.   if ITERLYAP==1
  208.      Lexp=lp;
  209.      Texp=t;
  210.   else
  211.      Lexp=[Lexp; lp];
  212.      Texp=[Texp; t];
  213.   end;

  214.   
  215.   if (mod(ITERLYAP,ioutp)==0)
  216.      for k=1:n_exp
  217.          txtstring{k}=['\lambda_' int2str(k) '=' num2str(lp(k))];
  218.      end
  219.      legend(Fig_Lyap_Axes,txtstring,3);
  220.   end;


  221. end;

  222.    ss=warndlg('Attention! Plot of lyapunov exponents will be closed!','Press OK to continue!');
  223.    uiwait(ss);
  224.    delete(Fig_Lyap);
  225.    
  226.      fprintf('\n \n Results of Lyapunov exponents calculation: \n t=%6.4f',t);
  227.      for k=1:n_exp fprintf(' L%d=%f; ',k,lp(k)); end;
  228.      fprintf('\n');
  229.    
  230.          if ~isempty(driver_window)
  231.              if ishandle(driver_window)
  232.                  delete(driver_window);
  233.                  driver_window = [];
  234.              end;         
  235.          end;
  236.          
  237.    calculation_progress = 0;

  238.    update_ds;
  239.    
复制代码
 楼主| 发表于 2006-9-2 06:11 | 显示全部楼层
 楼主| 发表于 2006-9-2 06:12 | 显示全部楼层
给出了分形计算的源代码的matlab程序,可以迅速帮助大家进行分形的分析与计算,参数容易设置

  1. function [Texp,Lexp]=new1lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp,d);
  2. %
  3. %    Lyapunov exponent calcullation for ODE-system.
  4. %
  5. %    The alogrithm employed in this m-file for determining Lyapunov
  6. %    exponents was proposed in
  7. %
  8. %         A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,
  9. %        "Determining Lyapunov Exponents from a Time Series," Physica D,
  10. %        Vol. 16, pp. 285-317, 1985.
  11. %
  12. %    For integrating ODE system can be used any MATLAB ODE-suite methods.
  13. % This function is a part of MATDS program - toolbox for dynamical system investigation
  14. %    See:    http://www.math.rsu.ru/mexmat/kvm/matds/
  15. %
  16. %    Input parameters:
  17. %      n - number of equation
  18. %      rhs_ext_fcn - handle of function with right hand side of extended ODE-system.
  19. %              This function must include RHS of ODE-system coupled with
  20. %              variational equation (n items of linearized systems, see Example).                  
  21. %      fcn_integrator - handle of ODE integrator function, for example: @ode45                  
  22. %      tstart - start values of independent value (time t)
  23. %      stept - step on t-variable for Gram-Schmidt renormalization procedure.
  24. %      tend - finish value of time
  25. %      ystart - start point of trajectory of ODE system.
  26. %      ioutp - step of print to MATLAB main window. ioutp==0 - no print,
  27. %              if ioutp>0 then each ioutp-th point will be print.
  28. %
  29. %    Output parameters:
  30. %      Texp - time values
  31. %      Lexp - Lyapunov exponents to each time value.
  32. %
  33. %    Users have to write their own ODE functions for their specified
  34. %    systems and use handle of this function as rhs_ext_fcn - parameter.      
  35. %
  36. %    Example. Lorenz system:
  37. %               dx/dt = sigma*(y - x)     = f1
  38. %               dy/dt = r*x - y - x*z = f2
  39. %               dz/dt = x*y - b*z     = f3
  40. %
  41. %    The Jacobian of system:
  42. %        | -sigma  sigma  0 |
  43. %    J = |   r-z    -1   -x |
  44. %        |    y      x   -b |
  45. %
  46. %    Then, the variational equation has a form:
  47. %
  48. %    F = J*Y
  49. %    where Y is a square matrix with the same dimension as J.
  50. %    Corresponding m-file:
  51. %        function f=lorenz_ext(t,X)
  52. %         SIGMA = 10; R = 28; BETA = 8/3;
  53. %         x=X(1); y=X(2); z=X(3);
  54. %
  55. %         Y= [X(4), X(7), X(10);
  56. %             X(5), X(8), X(11);
  57. %             X(6), X(9), X(12)];
  58. %         f=zeros(9,1);
  59. %         f(1)=SIGMA*(y-x); f(2)=-x*z+R*x-y; f(3)=x*y-BETA*z;
  60. %
  61. %         Jac=[-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA];
  62. %  
  63. %         f(4:12)=Jac*Y;
  64. %
  65. %    Run Lyapunov exponent calculation:
  66. %     
  67. %    [T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[0 1 0],10);   
  68. %   
  69. %    See files: lorenz_ext, run_lyap.   
  70. %  
  71. % --------------------------------------------------------------------
  72. % Copyright (C) 2004, Govorukhin V.N.
  73. % This file is intended for use with MATLAB and was produced for MATDS-program
  74. % http://www.math.rsu.ru/mexmat/kvm/matds/
  75. % lyapunov.m is free software. lyapunov.m is distributed in the hope that it
  76. % will be useful, but WITHOUT ANY WARRANTY.
  77. %



  78. %
  79. %       n=number of nonlinear odes
  80. %       n2=n*(n+1)=total number of odes
  81. %

  82. n1=n; n2=n1*(n1+1);

  83. %  Number of steps

  84. nit = round((tend-tstart)/stept);

  85. % Memory allocation

  86. y=zeros(n2,1); cum=zeros(n1,1); y0=y;
  87. gsc=cum; znorm=cum;

  88. % Initial values

  89. y(1:n)=ystart(:);

  90. for i=1:n1 y((n1+1)*i)=1.0; end;

  91. t=tstart;

  92. % Main loop

  93. for ITERLYAP=1:nit

  94. % Solutuion of extended ODE system

  95.   [T,Y] = unit(t,stept,y,d);  
  96.   
  97.   t=t+stept;
  98.   y=Y(size(Y,1),:);

  99.   for i=1:n1
  100.       for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;
  101.   end;

  102. %
  103. %       construct new orthonormal basis by gram-schmidt
  104. %

  105.   znorm(1)=0.0;
  106.   for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;

  107.   znorm(1)=sqrt(znorm(1));

  108.   for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;

  109.   for j=2:n1
  110.       for k=1:(j-1)
  111.           gsc(k)=0.0;
  112.           for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;
  113.       end;

  114.       for k=1:n1
  115.           for l=1:(j-1)
  116.               y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);
  117.           end;
  118.       end;

  119.       znorm(j)=0.0;
  120.       for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end;
  121.       znorm(j)=sqrt(znorm(j));

  122.       for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;
  123.   end;

  124. %
  125. %       update running vector magnitudes
  126. %

  127.   for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;

  128. %
  129. %       normalize exponent
  130. %

  131.   for k=1:n1
  132.       lp(k)=cum(k)/(t-tstart);
  133.   end;

  134. % Output modification

  135.   if ITERLYAP==1
  136.      Lexp=lp;
  137.      Texp=t;
  138.   else
  139.      Lexp=lp;
  140.      Texp=t;
  141.   end;


  142.   for i=1:n1
  143.       for j=1:n1
  144.           y(n1*j+i)=y0(n1*i+j);
  145.       end;
  146.   end;

  147. end;
复制代码
 楼主| 发表于 2006-9-2 06:14 | 显示全部楼层
小数据量法计算 Lyapunov 指数的 Matlab 程序 - (mex 函数,超快)

----------------------------------------------
参考文献:
Michael T.Rosenstein,
"A practical method for calculating largest lyapunov exponents from small sets",
Physica D,1993,65: 117-134

文件说明:
----------------------------------------------
1、LargestLyapunov_example1.m         程序主文件1,直接运行此文件即可
2、LargestLyapunov_example2.m        程序主文件2,直接运行此文件即可
3、LorenzData.dll                产生 Lorenz 离散数据
4、PhaSpaRecon.m                相空间重构
5、Lyapunov_luzhenbo.dll        Lyapunov 计算主函数
6、lyapunov_buffer.dll                Lyapunov 计算缓存


http://www.pudn.com/downloads63/sourcecode/math/detail221870.html
 楼主| 发表于 2006-9-2 06:14 | 显示全部楼层
计算lyapunov指数的程序

  1. program lylorenz
  2.         parameter(n=3,m=12,st=100)
  3.         integer::i,j,k
  4.     real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3
  5.         y(1)=10.
  6.         y(2)=1.
  7.         y(3)=0.
  8.         a=10.
  9.     b=8./3.
  10.     r=28.
  11.         t=0.
  12.         h=0.01

  13. !!!!!initial conditions
  14. do i=n+1,m
  15.         y(i)=0.
  16. end do
  17. do  i=1,n
  18.         y((n+1)*i)=1.
  19.     s(i)=0
  20. end do

  21. open(1,file='lorenz1.dat')
  22. open(2,file='ly lorenz.dat')
  23. do 100 k=1,st   !!!!!!!!st iterations
  24. call rgkt(m,h,t,y,f,yc,y1)

  25.     !!!!normarize vector 123
  26.         z=0.   
  27.     do i=1,n
  28.          do j=1,n
  29.         z(i)=z(i)+y(n*j+i)**2
  30.         enddo
  31.         if(z(i)>0.)z(i)=sqrt(z(i))
  32.         do j=1,n
  33.     y(n*j+i)=y(n*j+i)/z(i)
  34.         enddo
  35.     end do
  36.    !!!!generate gsr coefficient
  37.    k1=0.
  38.    k2=0.
  39.    k3=0.
  40. do i=1,n       
  41.         k1=k1+y(3*i+1)*y(3*i+2)
  42.         k2=k2+y(3*i+3)*y(3*i+2)
  43.         k3=k3+y(3*i+1)*y(3*i+3)
  44. end do
  45.    !!!!conduct new vector2 and 3
  46. do i=1,n
  47.         y(3*i+2)=y(3*i+2)-k1*y(3*i+1)       
  48.     y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)
  49. end do
  50.    
  51.     !!!generate new vector2 and 3,normarize them
  52. do i=2,n
  53.     z(i)=0.
  54.     do j=2,n
  55.         z(i)=z(i)+y(n*j+i)**2
  56.         enddo
  57.     if(z(i)>0.)z(i)=sqrt(z(i))
  58.         do j=2,n
  59.     y(n*j+i)=y(n*j+i)/z(i)
  60.         end do
  61. end do

  62.     !!!!!!!update lyapunov exponent
  63. do i=1,n
  64.    if(z(i)>0)s(i)=s(i)+log(z(i))
  65. enddo

  66. 100 continue
  67. do i=1,n
  68.   s(i)=s(i)/(1.*st*h*1000.)
  69. write(2,*)s(i)
  70. enddo
  71. end

  72. !!!!!!!!!!!!!!!!!!!!!!!!!
  73. subroutine rgkt(m,h,t,y,f,yc,y1)

  74. real y(m),f(m),y1(m),yc(m),a,b,r
  75. integer::i,j
  76. do j=1,1000
  77.         call df(m,t,y,f)
  78.         t=t+h/2.0
  79.         do i=1,m
  80.                 yc(i)=y(i)+h*f(i)/2.0
  81.                 y1(i)=y(i)+h*f(i)/6.0               
  82.                 end do                       
  83.         call df(m,t,yc,f)
  84.         do i=1,m
  85.                 yc(i)=y(i)+h*f(i)/2.0
  86.                 y1(i)=y1(i)+h*f(i)/3.0               
  87.         end do       
  88.         call df(m,t,yc,f)
  89.         t=t+h/2.0
  90.         do i=1,m
  91.                 yc(i)=y(i)+h*f(i)
  92.                 y1(i)=y1(i)+h*f(i)/3.0
  93.         end do
  94.         call df(m,t,yc,f)
  95.         do i=1,m
  96.                 y(i)=y1(i)+h*f(i)/6.0
  97.                 end do
  98. if(j>500)write(1,*)t,y(1),y(2),y(3)
  99. end do
  100. return
  101. end
  102. !!!!!!!!!!!!!!!!!!!!!!!!
  103. subroutine df(m,t,y,f)
  104. real y(m),a,b,r,f(m)
  105. common a,b,r
  106.     a=10.
  107.     b=8./3.
  108.     r=28.
  109.         f(1)=a*(y(2)-y(1))
  110.         f(2)=y(1)*(r-y(3))-y(2)
  111.         f(3)=y(1)*y(2)-b*y(3)

  112.     do  i=0,2
  113.         f(4+i)=a*y(7+i)-y(4+i)
  114.         f(7+i)=y(4+i)*(r-y(3))-y(7+i)-y(1)*y(10+i)
  115.         f(10+i)=y(2)*y(4+i)-b*y(10+i)+y(1)*y(7+i)
  116.     enddo
  117.         return
  118. end
复制代码
 楼主| 发表于 2006-9-2 06:16 | 显示全部楼层
C-C方法计算时间延迟和嵌入维数计算Lyapunov指数计算关联维数混沌时间序列预测

C-C方法计算时间延迟和嵌入维数
主程序:C_CMethod.m, C_CMethod_independent.m
子函数:correlation_integral.m(计算关联积分)
        disjoint.m(将时间序列拆分成t个不相关的子序列)
        heaviside.m(计算时间序列的海维赛函数值)
参考文献Nonlinear dynamics, delay times, and embedding windows。
计算Lyapunov指数:
largest_lyapunov_exponent.m(用吕金虎的方法计算最大Lyapunov指数)
参考文献:基于Lyapunov指数改进算法的边坡位移预测。
lyapunov_wolf.m(用wolf方法计算最大Lyapunov指数)
计算关联维数:G_P.m(G-P算法)

混沌时间序列预测
主函数
MainPre_by_jiaquanyijie_1.m(该程序用加权一阶局域法对数据进行进行一步预测)
MainPre_by_jiaquanyijie_n.m(该程序用加权一阶局域法对数据进行进行n步预测)
MainPre_by_Lya_1.m(基于最大Lyapunov指数的一步预测)
MainPre_by_Lya_n.m(基于最大Lyapunov指数的n步预测)
nearest_point.m(计算最后一个相点的最近相点的位置及最短距离)
子函数
jiaquanyijie.m(该函数用加权一阶局域法(xx)、零级近似(yy)和基于零级近似的加权一阶局域法(zz)对时间数据进行一步预测)
pre_by_lya.m(基于最大Lyapunov指数的预测方法)
pre_by_lya_new.m(改进的基于最大Lyapunov指数的预测方法)

如果需要该工具包,大家到http://www.pudn.com/downloads25/ ... rs/detail82624.html下载
 楼主| 发表于 2006-9-2 06:18 | 显示全部楼层
求取lyapunov指数的小数据量方法,采用混合编程

lylorenz.f90的程序如下:

program lylorenz
        parameter(n=3,m=12,st=100)
        integer::i,j,k
    real y(m),z(n),s(n),yc(m),h,y1(m),a,b,r,f(m),k1,k2,k3
        y(1)=10.
        y(2)=1.
        y(3)=0.
        a=10.
    b=8./3.
    r=28.
        t=0.
        h=0.01

!!!!!initial conditions
do i=n+1,m
        y(i)=0.
end do
do  i=1,n
        y((n+1)*i)=1.
    s(i)=0
end do

open(1,file='lorenz1.dat')
open(2,file='ly lorenz.dat')
do 100 k=1,st   !!!!!!!!st iterations
call rgkt(m,h,t,y,f,yc,y1)

    !!!!normarize vector 123
        z=0.   
    do i=1,n
         do j=1,n
        z(i)=z(i)+y(n*j+i)**2
        enddo
        if(z(i)>0.)z(i)=sqrt(z(i))
        do j=1,n
    y(n*j+i)=y(n*j+i)/z(i)
        enddo
    end do
   !!!!generate gsr coefficient
   k1=0.
   k2=0.
   k3=0.
do i=1,n       
        k1=k1+y(3*i+1)*y(3*i+2)
        k2=k2+y(3*i+3)*y(3*i+2)
        k3=k3+y(3*i+1)*y(3*i+3)
end do
   !!!!conduct new vector2 and 3
do i=1,n
        y(3*i+2)=y(3*i+2)-k1*y(3*i+1)       
    y(3*i+3)=y(3*i+3)-k2*y(3*i+2)-k3*y(3*i+1)
end do
   
    !!!generate new vector2 and 3,normarize them
do i=2,n
    z(i)=0.
    do j=2,n
        z(i)=z(i)+y(n*j+i)**2
        enddo
    if(z(i)>0.)z(i)=sqrt(z(i))
        do j=2,n
    y(n*j+i)=y(n*j+i)/z(i)
        end do
end do

    !!!!!!!update lyapunov exponent
do i=1,n
   if(z(i)>0)s(i)=s(i)+log(z(i))
enddo

100 continue
do i=1,n
  s(i)=s(i)/(1.*st*h*1000.)
write(2,*)s(i)
enddo
end

!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine rgkt(m,h,t,y,f,yc,y1)

real y(m),f(m),y1(m),yc(m),a,b,r
integer::i,j
do j=1,1000
        call df(m,t,y,f)
        t=t+h/2.0
        do i=1,m
                yc(i)=y(i)+h*f(i)/2.0
                y1(i)=y(i)+h*f(i)/6.0               
                end do                       
        call df(m,t,yc,f)
        do i=1,m
                yc(i)=y(i)+h*f(i)/2.0
                y1(i)=y1(i)+h*f(i)/3.0               
        end do       
        call df(m,t,yc,f)
        t=t+h/2.0
        do i=1,m
                yc(i)=y(i)+h*f(i)
                y1(i)=y1(i)+h*f(i)/3.0
        end do
        call df(m,t,yc,f)
        do i=1,m
                y(i)=y1(i)+h*f(i)/6.0
                end do
if(j>500)write(1,*)t,y(1),y(2),y(3)
end do
return
end
!!!!!!!!!!!!!!!!!!!!!!!!
subroutine df(m,t,y,f)
real y(m),a,b,r,f(m)
common a,b,r
    a=10.
    b=8./3.
    r=28.
        f(1)=a*(y(2)-y(1))
        f(2)=y(1)*(r-y(3))-y(2)
        f(3)=y(1)*y(2)-b*y(3)

    do  i=0,2
        f(4+i)=a*y(7+i)-y(4+i)
        f(7+i)=y(4+i)*(r-y(3))-y(7+i)-y(1)*y(10+i)
        f(10+i)=y(2)*y(4+i)-b*y(10+i)+y(1)*y(7+i)
    enddo
        return
end

评分

1

查看全部评分

 楼主| 发表于 2006-9-2 06:21 | 显示全部楼层
计算各种混沌系统李雅普洛夫指数的MATLAB 源程序

http://www.pudn.com/downloads50/sourcecode/math/detail172856.html

评分

1

查看全部评分

发表于 2006-9-13 15:09 | 显示全部楼层
正好帮上忙了
让我看看程序的正确性怎么样
先表示感谢
这个工作量是很大的,辛苦了
发表于 2006-9-20 20:45 | 显示全部楼层
想问一下对于不连续的非线性方程,有没有计算程序,
我知道有文章讨论这个问题
P.Muller,Chaos Solitions Fractals 5(1995),1671
但还找不到程序
发表于 2006-10-7 18:22 | 显示全部楼层
离散问题的论坛也有好几个程序,自己搜索一下吧

[ 本帖最后由 ChaChing 于 2010-5-11 10:36 编辑 ]
发表于 2006-10-9 20:46 | 显示全部楼层

应用中出现问题,向大侠求救ing

各位好,我在用楼主所提供的C-C方法求时间延迟时遇到了下面的错误提示,请各位帮小弟诊断一下,谢谢!
t =
     1
t =
     2
??? Index exceeds matrix dimensions.

有朋友试过这个C-C这个程序吗?怎么我试了下出现问题的 ,

[ 本帖最后由 ChaChing 于 2010-5-11 10:35 编辑 ]
发表于 2006-10-11 09:32 | 显示全部楼层
不知道你给的data是什么样的?我这里没有这样的序列

从错误提示上看估计是因为提供序列有问题造成的

[ 本帖最后由 ChaChing 于 2010-5-11 10:36 编辑 ]
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表