在进行lyapunov指数的求取时&＃xff0c;需要知道离散动力学系统对应Jacobi矩阵的特征值&＃xff0c;qr法与Jacobi法都可以求解矩阵特征值&＃xff0c;其中qr法求解的是矩阵所有特征值&＃xff0c;而Jacobi法求解的是矩阵的最大特征值。本文以二维Henon映射为例&＃xff0c;分别展示两种方法在求解时的区别与联系。

目录

1.准备工作

1.1 henon映射

动力学系统

代码实现

 1.2 jacobi矩阵

1.3 lyapunov指数

2.Jacobi法

 代码实现 

3.qr法

 代码实现

 lypunov指数谱&＃xff1a;

4.qr法与Jacobi法之间代码的联系

代码实现

1.准备工作

1.1 henon映射

动力学系统

Henon映射为二维映射&＃xff0c;其动力学方程通常有以下两种写法&＃xff1a;

 和

 式中a、b为参数&＃xff0c;0

代码实现

clc;clear all;close all
%初值设置
x0&＃61;0.2;y0&＃61;0.3;
%参数设置
b&＃61;0.3;n&＃61;800;t&＃61;600;
hold on
for a &＃61; 0:0.01:1.4 %参数a
x(1)&＃61;1&＃43;y0-a*x0^2;
y(1)&＃61;b*x0;
for i &＃61;2:nx(i)&＃61;1&＃43;y(i-1)-a*x(i-1)^2;y(i)&＃61;b*x(i-1);
end
xn&＃61;x(t&＃43;1:n);%取80次迭代之后的数据
pause(0.1);
plot(a,xn,&＃39;r.&＃39;,&＃39;Markersize&＃39;,2);
xlabel(&＃39;a&＃39;);ylabel(&＃39;x&＃39;);
set(gca,&＃39;XLim&＃39;,[0 1.4]);
set(gca,&＃39;XTick&＃39;,[0:0.2:1.4]);
set(gca,&＃39;YLim&＃39;,[-1.5 1.5]);
set(gca,&＃39;YTick&＃39;,[-1.5:0.5:1.5]);
end


混沌图像&＃xff1a;

图1 .henon映射 

 1.2 jacobi矩阵

 求Jacobi矩阵&＃xff0c;实际上就是对动力学系统方程f关于x求偏导。

1.3 lyapunov指数

在式(3)jacobi矩阵的基础上有

首先令&＃xff1a;

Jk的m个复根可表示为&＃xff1a;

 离散动力学系统的第i个lypunov指数可表示为&＃xff1a;

2.Jacobi法

Jacobi法求特征值&＃xff0c;是求的最大特征值。故在(4)式的基础上&＃xff0c;我们应将特征值进行排序&＃xff0c;然后取出最大的特征值即可&＃xff0c;但是&＃xff1a;

①在求lyapunov指数是特征值取模&＃xff1b;

②对于复特征值而言不取模的话并无意义&＃xff1b;

③对于单纯数大&＃xff0c;但数值并不大的特征值而言&＃xff0c;其对应求得的lyapunov指数并没有十分特殊的意义&＃xff1b;

④特征值模最大&＃xff0c;可以对应求得系统的最大lyapunov指数&＃xff0c;对于henon映射这种2维系统而言&＃xff0c;对最大lyapunov指数进行观察就可以了解系统混沌与否。

故我们对(4)式中的特征值&＃xff0c;进行取模排序&＃xff0c;依次从小到大排序为&＃xff1a;

 那么利用jacobi法求系统特征值后所对应的lyapunov指数可表示为&＃xff1a;

 代码实现 

clc, clear
a&＃61;0:0.01:1.6; n&＃61;500; S&＃61;[];
for j&＃61;1:length(a)x&＃61;0.2; y&＃61;0.3; %x,y的初值for i&＃61;1:nx2&＃61;1-a(j)*x^2&＃43;y; y&＃61;0.3*x; x&＃61;x2; %首先进行序列迭代endif x(end)>-100 & x(end)<100 %若不发散&＃xff0c;再计算指数x&＃61;0.2; y&＃61;0.3;JJ&＃61;eye(2);%2阶单位矩阵for i&＃61;1:nx2&＃61;1-a(j)*x^2&＃43;y; y&＃61;0.3*x; x&＃61;x2;J&＃61;[-2*a(j)*x, 1; 0.3, 0];%jacobi矩阵JJ&＃61;JJ*J;endL&＃61;eigs(JJ,1); %求模最大的特征值,并返回1个最大特征值S&＃61;[S,[a(j);log(abs(L))/n]]; %把a值及指数值保存end
end
plot(S(1,:),S(2,:),&＃39;.-&＃39;) %画出a对应的最大Lyapunov指数
hold on, plot([a(1),a(end)],[0,0],&＃39;k&＃39;)
xlabel(&＃39;\it a&＃39;), ylabel(&＃39;最大Lyapunov指数&＃39;)


 lyapunov指数谱&＃xff1a;

图2 .jacobi法求得的henon映射lypunov指数谱

&＃xff08;henon映射最大Lyapunov指数谱&＃xff09;

 jacobi法求得的lypunov指数谱与混沌的对应&＃xff1a;

图3 .jacobi法求得的henon映射与最大lyapunov指数之间的对应 

3.qr法

qr法求的是系统所有的特征值&＃xff0c;即把式(4)中所有特征值&＃xff0c;进行取模后全部转化为lyapunov。由于henon映射为2维系统&＃xff0c;故其对应的特征值最多就2个&＃xff1a;

 代码实现

clc;close all;clear all
N &＃61; 1000;
a &＃61; (0:0.001:1.4)&＃39;;
b &＃61; 0.3;
na &＃61; length(a);
LE1 &＃61; zeros(na,1);
LE2 &＃61; zeros(na,1);
x &＃61; 0.2; y &＃61; 0.3;
for i&＃61;1:na LCEvector &＃61; zeros(2,1); Q &＃61; eye(2); for j&＃61;1:N xprev &＃61; x; yprev &＃61; y; x &＃61; 1-a(i)*xprev*xprev&＃43;yprev; y &＃61; b*xprev; Ji &＃61; [-2*a(i)*x,1;b 0];B &＃61; Ji*Q;[Q,R] &＃61; qr(B); %R得B矩阵的上三角矩阵LCEvector &＃61; LCEvector&＃43;log(diag(abs(R))); end LE &＃61; LCEvector/N; LE1(i) &＃61; LE(1); LE2(i) &＃61; LE(2);
end
figure(1)
plot(a,LE1,&＃39;g&＃39;,&＃39;linewidth&＃39;,1) ;
hold on
plot(a,LE2,&＃39;b&＃39;,&＃39;linewidth&＃39;,1);
set(gca,&＃39;XLim&＃39;,[0 1.4]);set(gca,&＃39;YLim&＃39;,[-2 1]);
legend(&＃39;\lambda1&＃39;,&＃39;\lambda2&＃39;);
hold on
plot([0 1.4],[0 0],&＃39;k--&＃39;,&＃39;linewidth&＃39;,1)
hold off xlabel(&＃39;a&＃39;);ylabel(&＃39;LE&＃39;);set(gca,&＃39;fontsize&＃39;,10)


 lyapunov指数谱&＃xff1a;

图4 .qr法求得的henon映射lyapunov指数谱

&＃xff08;henon映射Lyapunov指数谱,2维系统&＃xff0c;一共两条lyapunov&＃xff09;

4.qr法与Jacobi法之间代码的联系

在jacobi法代码基础上&＃xff0c;我们稍作调整也可以得qr法对应的lyapunov指数谱&＃xff1a;

①去掉eigs()这个求最大特征值函数

②引入qr函数

③jacobi法每次对应1个特征值&＃xff0c;qr法此处对应两个&＃xff0c;因此为了累加计算的方便将lyapunov指数迭代放到循环里面。

代码实现

clc, clear
a&＃61;0:0.001:1.4; n&＃61;500;
for j&＃61;1:length(a)LE &＃61; zeros(2,1); x&＃61;0.2; y&＃61;0.3; %x,y的初值for i&＃61;1:nx2&＃61;1-a(j)*x^2&＃43;y; y&＃61;0.3*x; x&＃61;x2; %首先进行序列迭代endif x(end)>-100 & x(end)<100 %若不发散&＃xff0c;再计算指数x&＃61;0.2; y&＃61;0.3;JJ&＃61;eye(2);%2阶单位矩阵for i&＃61;1:nx2&＃61;1-a(j)*x^2&＃43;y; y&＃61;0.3*x; x&＃61;x2;J&＃61;[-2*a(j)*x, 1; 0.3, 0];%jacobi矩阵JJ&＃61;J*JJ;[JJ,L]&＃61;qr(JJ);LE &＃61; LE&＃43;log(diag(abs(L)));endendLE &＃61; LE/n;LE1(j) &＃61; LE(1);LE2(j) &＃61; LE(2);
end
figure(1)
plot(a,LE1,&＃39;g&＃39;,&＃39;linewidth&＃39;,1) ;
hold on
plot(a,LE2,&＃39;b&＃39;,&＃39;linewidth&＃39;,1);
set(gca,&＃39;XLim&＃39;,[0 1.4]);set(gca,&＃39;YLim&＃39;,[-2 1]);
legend(&＃39;\lambda1&＃39;,&＃39;\lambda2&＃39;);
hold on
plot([0 1.4],[0 0],&＃39;k--&＃39;,&＃39;linewidth&＃39;,1)
hold off xlabel(&＃39;a&＃39;);ylabel(&＃39;LE&＃39;);set(gca,&＃39;fontsize&＃39;,10)


 lyapunov指数谱&＃xff1a;