关于lyapunov指数的问题

octopussheng

用“计算连续方程Lyapunov指数的程序”算了一个单自由度系统振动的lyapunov指数，
系统是一个二阶的微分方程，算出了有三个Lyapunov指数曲线，请问有没有人可以帮我解释一下，这三个是不是都是有效的，还是只有一个有用？？

这个在Lyapunov指数的定义中好像没有看到过啊！！

请帮忙讨论一下，谢谢！

方程如下：
x''=-A*omiga^2*sin(omiga*t)-B0-B1*x-B2*x^3-C1*x'-C2*x'^2

参数这里就不给出了，可以自己定义的！

gghhjj

三个Lyapunov指数？能否说明一下什么意思？

octopussheng

对 啊，我就是不明白用这个程序怎么算出三个Lyapunov指数的啊！！
请各位懂的朋友多多帮忙啊！！

VibInfo

原帖由 octopussheng 于 2007-4-25 07:39 发表
对 啊，我就是不明白用这个程序怎么算出三个Lyapunov指数的啊！！
请各位懂的朋友多多帮忙啊！！

那你就需要把程序给出，否则别人怎么知道你算得是什么

octopussheng

代码如下：

1. %
   
2. % 计算Rossler吸引子的Lyapunov指数
   
3. %  2004.10.20 zya
   
4. %
   
5. 
   
6. clear;
   
7. yinit = [1,1,0.5];
   
8. orthmatrix = [1 0 0;
   
9.               0 1 0;
   
10.               0 0 1];
    
11. a = 0.15;
    
12. b = 0.20;
    
13. c = 10.0;
    
14. y = zeros(12,1);
    
15. % 初始化输入
    
16. y(1:3) = yinit;
    
17. y(4:12) = orthmatrix;
    
18. 
    
19. tstart = 0; % 时间初始值
    
20. tstep = 1e-3; % 时间步长
    
21. wholetimes = 1e5; % 总的循环次数
    
22. steps = 10; % 每次演化的步数
    
23. iteratetimes = wholetimes/steps; % 演化的次数
    
24. 
    
25. mod = zeros(3,1);
    
26. lp = zeros(3,1);
    
27. % 初始化三个Lyapunov指数
    
28. Lyapunov1 = zeros(iteratetimes,1);
    
29. Lyapunov2 = zeros(iteratetimes,1);
    
30. Lyapunov3 = zeros(iteratetimes,1);
    
31. 
    
32. for i=1:iteratetimes
    
33.     tspan = tstart:tstep:(tstart + tstep*steps);   
    
34.     [T,Y] = ode45('Rossler_ly', tspan, y);
    
35.     % 取积分得到的最后一个时刻的值
    
36.     y = Y(size(Y,1),:);
    
37.     % 重新定义起始时刻
    
38.     tstart = tstart + tstep*steps;
    
39.     y0 = [y(4) y(7) y(10);
    
40.           y(5) y(8) y(11);
    
41.           y(6) y(9) y(12)];
    
42.     %正交化
    
43.     y0 = FourGS(y0);
    
44.     % 取三个向量的模
    
45.     mod(1) = sqrt(y0(:,1)'*y0(:,1));
    
46.     mod(2) = sqrt(y0(:,2)'*y0(:,2));
    
47.     mod(3) = sqrt(y0(:,3)'*y0(:,3));
    
48.     y0(:,1) = y0(:,1)/mod(1);
    
49.     y0(:,2) = y0(:,2)/mod(2);
    
50.     y0(:,3) = y0(:,3)/mod(3);
    
51.    
    
52.     lp = lp+log(abs(mod));
    
53.     %三个Lyapunov指数
    
54.     Lyapunov1(i) = lp(1)/(tstart);
    
55.     Lyapunov2(i) = lp(2)/(tstart);
    
56.     Lyapunov3(i) = lp(3)/(tstart);
    
57.    
    
58.     y(4:12) = y0';
    
59. end
    
60. % 作Lyapunov指数谱图
    
61. i = 1:iteratetimes;
    
62. plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)

复制代码

octopussheng

Rosser吸引子程序：

1. function dX = Rossler_ly(t,X)
   
2. %
   
3. %  Rossler吸引子，用来计算Lyapunov指数
   
4. %        a=0.15,b=0.20,c=10.0
   
5. %        dx/dt = -y-z,
   
6. %        dy/dt = x+ay,
   
7. %        dz/dt = b+z(x-c),
   
8. %  2004.10.20 zya
   
9. %
   
10. 
    
11. a = 0.15;
    
12. b = 0.20;
    
13. c = 10.0;
    
14. 
    
15. x=X(1); y=X(2); z=X(3);
    
16. % Y的三个列向量为相互正交的单位向量
    
17. Y = [X(4), X(7), X(10);
    
18.     X(5), X(8), X(11);
    
19.     X(6), X(9), X(12)];
    
20. % 输出向量的初始化，必不可少
    
21. dX = zeros(12,1);
    
22. 
    
23. % Rossler吸引子
    
24. dX(1) = -y-z;
    
25. dX(2) = x+a*y;
    
26. dX(3) = b+z*(x-c);
    
27. 
    
28. % Rossler吸引子的Jacobi矩阵
    
29. Jaco = [0 -1 -1;
    
30.         1 a  0;
    
31.         z 0  x-c];
    
32.         
    
33. dX(4:12) = Jaco*Y;

复制代码

octopussheng

还有中间用到的那个FourGS程序：

1. % 四维G-S变换
   
2. function A = FourGS(V)  % V 为4*4向量
   
3. v1 = V(:,1);
   
4. v2 = V(:,2);
   
5. v3 = V(:,3);
   
6. v4 = V(:,4);
   
7. a1 = zeros(4,1);
   
8. a2 = zeros(4,1);
   
9. a3 = zeros(4,1);
   
10. a4 = zeros(4,1);
    
11. a1 = v1;
    
12. a2 = v2-((a1'*v2)/(a1'*a1))*a1;
    
13. a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
    
14. a4 = v4-((a1'*v4)/(a1'*a1))*a1-((a2'*v4)/(a2'*a2))*a2-((a3'*v4)/(a3'*a3))*a3;
    
15. A = [a1,a2,a3,a4];
    

复制代码

其实这些都是论坛里面的！

shenyongjun

这个方程是应该有三个Lyapunov指数的！
虽然这个方程是二阶微分方程，但是含有强迫激励，因此在计算Lyapunov指数时，转换为状态空间应该是三维的，那么最后就应该有三个Lyapunov指数。

hohoo

由方程x''=-A*omiga^2*sin(omiga*t)-B0-B1*x-B2*x^3-C1*x'-C2*x'^2
改写为
dx=y;
dy=-A*omiga^2*sin(omiga*z)-B0-B1*x-B2*x^3-c1*y-c2*y^2;
dz=1;
对应的Jocobi矩阵也是3*3的
得到三个Lyapunov指数

net22446688

弱弱的问一句，什么是lyapunov指数？

gghhjj

原帖由 net22446688 于 2007-4-25 18:19 发表
弱弱的问一句，什么是lyapunov指数？

判断运动稳定性的方法，任何非线性动力学的书上都会介绍到，以下是转贴

1. The hardest thing to get right about Lyapunov exponents is the spelling of Lyapunov, which you will variously find as Liapunov, Lyapunof and even Liapunoff. Of course Lyapunov is really spelled in the Cyrillic(古代斯拉夫语) alphabet: (Lambda)(backwards r)(pi)(Y)(H)(0)(B). Now that there is an ANSI standard of transliteration for Cyrillic, we expect all references to converge on the version Lyapunov.
   
2. 
   
3. Lyapunov was born in Russia in 6 June 1857. He was greatly influenced by Chebyshev and was a student with Markov. He was also a passionate(充满热情的) man: Lyapunov shot himself the day his wife died. He died 3 Nov. 1918, three days later. According to the request on a note he left, Lyapunov was buried with his wife.
   
4. Lyapunov left us with more than just a simple note. He left a collection of papers on the equilibrium shape of rotating liquids, on probability, and on the stability of low-dimensional dynamical systems. It was from his dissertation(论文) that the notion of Lyapunov exponent emerged. Lyapunov was interested in showing how to discover if a solution to a dynamical system is stable or not for all times. The usual method of studying stability, i.e. linear stability, was not good enough, because if you waited long enough the small errors due to linearization would pile up and make the approximation invalid.
   
5. Lyapunov developed concepts (now called Lyapunov Stability) to overcome these difficulties. Lyapunov exponents measure the rate at which nearby orbits converge or diverge. There are as many Lyapunov exponents as there are dimensions in the state space of the system, but the largest is usually the most important. Roughly speaking the (maximal) Lyapunov exponent is the time constant, lambda, in the expression for the distance between two nearby orbits, exp(lambda * t).! If lambda is negative, then the orbits converge in time, and the dynamical system is insensitive to initial conditions.! However, if lambda is positive, then the distance between nearby orbits grows exponentially in time, and the system exhibits sensitive dependence on initial conditions.
   
6. 
   
7. There are basically two ways to compute Lyapunov exponents. In one way one chooses two nearby points, evolves them in time, measuring the growth rate of the distance between them. This is useful when one has a time series, but has the disadvantage that the growth rate is really not a local effect as the points separate. A better way is to measure the growth rate of tangent vectors to a given orbit.
   
8. More precisely, consider a map f in an m dimensional phase space, and its derivative matrix Df(x). Let v be a tangent vector at the point x. Then we define a function
   
9. L(x,v) = lim (1/n)ln |( Dnf (x)v )|
   
10.         n -> ￥
    
11. 
    
12. Now the Multiplicative Ergodic Theorem of Oseledec states that this limit exists for almost all points x and all tangent vectors v. There are at most m distinct values of L as we let v range over the tangent space. These are the Lyapunov exponents at x.

复制代码

gghhjj

原帖由 shenyongjun 于 2007-4-25 16:03 发表
这个方程是应该有三个Lyapunov指数的！
虽然这个方程是二阶微分方程，但是含有强迫激励，因此在计算Lyapunov指数时，转换为状态空间应该是三维的，那么最后就应该有三个Lyapunov指数。

1.     % 取三个向量的模
   
2. 
   
3.     mod(1) = sqrt(y0(:,1)'*y0(:,1));
   
4. 
   
5.     mod(2) = sqrt(y0(:,2)'*y0(:,2));
   
6. 
   
7.     mod(3) = sqrt(y0(:,3)'*y0(:,3));
   
8. 
   
9.     y0(:,1) = y0(:,1)/mod(1);
   
10. 
    
11.     y0(:,2) = y0(:,2)/mod(2);
    
12. 
    
13.     y0(:,3) = y0(:,3)/mod(3);
    
14. 
    
15.    
    
16. 
    
17.     lp = lp+log(abs(mod));
    
18. 
    
19.     %三个Lyapunov指数
    
20. 
    
21.     Lyapunov1(i) = lp(1)/(tstart);
    
22. 
    
23.     Lyapunov2(i) = lp(2)/(tstart);
    
24. 
    
25.     Lyapunov3(i) = lp(3)/(tstart);

复制代码

看这一段程序，已经说明的非常清楚了

octopussheng

感谢gghhjj的指教！

那这样的话，按照理解，是不是只有第一个Lyapunov指数的值才是对于判断系统稳定性有用的啊？？

gghhjj

原帖由 octopussheng 于 2007-4-26 08:01 发表
感谢gghhjj的指教！

那这样的话，按照理解，是不是只有第一个Lyapunov指数的值才是对于判断系统稳定性有用的啊？？

要看最大Lyapunov指数

octopussheng

谢谢gghhjj的指点！！
了解了！:lol