Rossler吸引子的M文件怎么写？

playfish

继续实验Huang的论文。前天的duffing方程已经成功。原来是一不小心的失误。但huang的论文里最后一个实验就不知道该怎么写了。Rossler吸引子的方程是：
dX/dY=-（Y+Z）
dY/dt=X+aY
dZ/dT=b+Z(X-r)
其中a=0.2，b=0.2，r=3.5，X(0)=-2；Y(0)=3；Z(0)=0；
不知道怎么编程序求解这个方程？

playfish

没人知道吗？

Vickyvictoria

1. % Author: Thomas Lee
   
2. % E-mail: lixf1979@126.com % Corresponding: School of Mathematics, Physics
   
3. % and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China
   
4. function dX = Rossler_ly(t,X)
   
5. %  Rossler吸引子，用来计算Lyapunov指数
   
6. %         a=0.15,b=0.20,c=10.0
   
7. %         dx/dt = -y-z,
   
8. %         dy/dt = x+ay,
   
9. %         dz/dt = b+z(x-c),
   
10. a = 0.15;
    
11. b = 0.20;
    
12. c = 10.0;
    
13. x=X(1); y=X(2); z=X(3);
    
14. % Y的三个列向量为相互正交的单位向量
    
15. Y = [X(4), X(7), X(10);
    
16.      X(5), X(8), X(11);
    
17.      X(6), X(9), X(12)];
    
18. % 输出向量的初始化，必不可少
    
19. dX = zeros(12,1);
    
20. % Rossler吸引子
    
21. dX(1) = -y-z;
    
22. dX(2) = x+a*y;
    
23. dX(3) = b+z*(x-c);
    
24. % Rossler吸引子的Jacobi矩阵
    
25. Jaco = [0 -1 -1;
    
26.          1 a  0;
    
27.          z 0  x-c];
    
28. dX(4:12) = Jaco*Y;

复制代码

求解LE代码：

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

复制代码

程序中用到的ThreeGS程序如下：

1. %G-S正交化
   
2. function A = ThreeGS(V)  % V 为3*3向量
   
3. v1 = V(:,1);
   
4. v2 = V(:,2);
   
5. v3 = V(:,3);
   
6. a1 = zeros(3,1);
   
7. a2 = zeros(3,1);
   
8. a3 = zeros(3,1);
   
9. a1 = v1;
   
10. a2 = v2-((a1'*v2)/(a1'*a1))*a1;
    
11. a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
    
12. A = [a1,a2,a3];

复制代码

计算得到的Rossler系统的LE为————  0.063231  0.092635  -9.8924

杰式悦敏

学习一下 谢谢分享