fork download
copy 
  1. function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
  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://w...content-available-to-author-only...u.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://w...content-available-to-author-only...u.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.  
  80.  
  81. %
  82. %       n=number of nonlinear odes
  83. %       n2=n*(n+1)=total number of odes
  84. %
  85.  
  86. n1=n; n2=n1*(n1+1);
  87.  
  88. %  Number of steps
  89.  
  90. nit = round((tend-tstart)/stept);
  91.  
  92. % Memory allocation 
  93.  
  94. y=zeros(n2,1); cum=zeros(n1,1); y0=y;
  95. gsc=cum; znorm=cum;
  96.  
  97. % Initial values
  98.  
  99. y(1:n)=ystart(:);
  100.  
  101. for i=1:n1 y((n1+1)*i)=1.0; end;
  102.  
  103. t=tstart;
  104.  
  105. % Main loop
  106.  
  107. for ITERLYAP=1:nit
  108.  
  109. % Solutuion of extended ODE system 
  110.  
  111.   [T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);  
  112.  
  113.   t=t+stept;
  114.   y=Y(size(Y,1),:);
  115.  
  116.   for i=1:n1 
  117.       for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;
  118.   end;
  119.  
  120. %
  121. %       construct new orthonormal basis by gram-schmidt
  122. %
  123.  
  124.   znorm(1)=0.0;
  125.   for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;
  126.  
  127.   znorm(1)=sqrt(znorm(1));
  128.  
  129.   for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;
  130.  
  131.   for j=2:n1
  132.       for k=1:(j-1)
  133.           gsc(k)=0.0;
  134.           for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;
  135.       end;
  136.  
  137.       for k=1:n1
  138.           for l=1:(j-1)
  139.               y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);
  140.           end;
  141.       end;
  142.  
  143.       znorm(j)=0.0;
  144.       for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end;
  145.       znorm(j)=sqrt(znorm(j));
  146.  
  147.       for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;
  148.   end;
  149.  
  150. %
  151. %       update running vector magnitudes
  152. %
  153.  
  154.   for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;
  155.  
  156. %
  157. %       normalize exponent
  158. %
  159.  
  160.   for k=1:n1 
  161.       lp(k)=cum(k)/(t-tstart); 
  162.   end;
  163.  
  164. % Output modification
  165.  
  166.   if ITERLYAP==1
  167.      Lexp=lp;
  168.      Texp=t;
  169.   else
  170.      Lexp=[Lexp; lp];
  171.      Texp=[Texp; t];
  172.   end;
  173.  
  174.   if (mod(ITERLYAP,ioutp)==0)
  175.      fprintf('t=%6.4f',t);
  176.      for k=1:n1 fprintf(' %10.6f',lp(k)); end;
  177.      fprintf('\n');
  178.   end;
  179.  
  180.   for i=1:n1 
  181.       for j=1:n1
  182.           y(n1*j+i)=y0(n1*i+j);
  183.       end;
  184.   end;
  185.  
  186. end;
Success #stdin #stdout #stderr 0.17s 65056KB
comments ()
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
parse error:

  function body open at end of input
https://ideone.com/KJhpp3
language:
Octave (octave 4.4.1)
created:
11 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!