nedu-fly
助理工程师 精华 0 积分 23 帖子 38 水位 76 技术分 0 | 关于定义法求解的程序,和matlab板块的“连续系统LE求解程序”差不多。以Rossler系统为例
Rossler系统微分方程定义程序 function dX = Rossler_ly(t,X) % Rossler吸引子,用来计算Lyapunov指数 % a=0.15,b=0.20,c=10.0 % dx/dt = -y-z, % dy/dt = x+ay, % dz/dt = b+z(x-c), a = 0.15; b = 0.20; c = 10.0; x=X(1); y=X(2); z=X(3); % Y的三个列向量为相互正交的单位向量 Y = [X(4), X(7), X(10); X(5), X(8), X(11); X(6), X(9), X(12)]; % 输出向量的初始化,必不可少 dX = zeros(12,1); % Rossler吸引子 dX(1) = -y-z; dX(2) = x+a*y; dX(3) = b+z*(x-c); % Rossler吸引子的Jacobi矩阵 Jaco = [0 -1 -1; 1 a 0; z 0 x-c]; dX(4:12) = Jaco*Y;
求解LE代码: % 计算Rossler吸引子的Lyapunov指数 clear; yinit = [1,1,1]; orthmatrix = [1 0 0; 0 1 0; 0 0 1]; a = 0.15; b = 0.20; c = 10.0; y = zeros(12,1); % 初始化输入 y(1:3) = yinit; y(4:12) = orthmatrix; tstart = 0; % 时间初始值 tstep = 1e-3; % 时间步长 wholetimes = 1e5; % 总的循环次数 steps = 10; % 每次演化的步数 iteratetimes = wholetimes/steps; % 演化的次数 mod = zeros(3,1); lp = zeros(3,1); % 初始化三个Lyapunov指数 Lyapunov1 = zeros(iteratetimes,1); Lyapunov2 = zeros(iteratetimes,1); Lyapunov3 = zeros(iteratetimes,1); for i=1:iteratetimes tspan = tstart:tsteptstart + tstep*steps); [T,Y] = ode45('Rossler_ly', tspan, y); % 取积分得到的最后一个时刻的值 y = Y(size(Y,1),; % 重新定义起始时刻 tstart = tstart + tstep*steps; y0 = [y(4) y(7) y(10); y(5) y(8) y(11); y(6) y(9) y(12)]; %正交化 y0 = ThreeGS(y0); % 取三个向量的模 mod(1) = sqrt(y0(:,1)'*y0(:,1)); mod(2) = sqrt(y0(:,2)'*y0(:,2)); mod(3) = sqrt(y0(:,3)'*y0(:,3)); y0(:,1) = y0(:,1)/mod(1); y0(:,2) = y0(:,2)/mod(2); y0(:,3) = y0(:,3)/mod(3); lp = lp+log(abs(mod)); %三个Lyapunov指数 Lyapunov1(i) = lp(1)/(tstart); Lyapunov2(i) = lp(2)/(tstart); Lyapunov3(i) = lp(3)/(tstart); y(4:12) = y0'; end % 作Lyapunov指数谱图 i = 1:iteratetimes; plot(i,Lyapunov1,i,Lyapunov2,i,Lyapunov3)
程序中用到的ThreeGS程序如下: %G-S正交化 function A = ThreeGS(V) % V 为3*3向量 v1 = V(:,1); v2 = V(:,2); v3 = V(:,3); a1 = zeros(3,1); a2 = zeros(3,1); a3 = zeros(3,1); a1 = v1; a2 = v2-((a1'*v2)/(a1'*a1))*a1; a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2; A = [a1,a2,a3];
计算得到的Rossler系统的LE为―――― 0.063231 0.092635 -9.8924 Wolf文章中计算得到的Rossler系统的LE为――――0.09 0 -9.77
需要注意的是――定义法求解的精度有限,对有些系统的计算往往出现计果和理论值有偏差的现象。
正交化程序可以根据上面的扩展到N*N向量,这里就不加以说明了,对matlab用户来说应该还是比较简单的!
|
|