0
点赞
收藏
分享

微信扫一扫

lyapunov指数求取时运用qr法与jacobi法之间的区别与联系【基于matlab的动力学模型学习笔记_9】

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

目录

1.准备工作

1.1 henon映射

动力学系统

代码实现

 1.2 jacobi矩阵

1.3 lyapunov指数

2.Jacobi法

 代码实现 

3.qr法

 代码实现

 lypunov指数谱:

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

代码实现


1.准备工作

1.1 henon映射

动力学系统

Henon映射为二维映射,其动力学方程通常有以下两种写法:

 和

 式中a、b为参数,0<b≤1我们通常运用式(1)表达形式。

代码实现

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

混沌图像:

图1 .henon映射 

 1.2 jacobi矩阵

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

1.3 lyapunov指数

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

首先令:

Jk的m个复根可表示为:

 离散动力学系统的第i个lypunov指数可表示为:

2.Jacobi法

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

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

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

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

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

故我们对(4)式中的特征值,进行取模排序,依次从小到大排序为:

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

 代码实现 

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

 lyapunov指数谱:

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

(henon映射最大Lyapunov指数谱

 jacobi法求得的lypunov指数谱与混沌的对应:

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

3.qr法

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

 代码实现

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

 lyapunov指数谱:

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

(henon映射Lyapunov指数谱,2维系统,一共两条lyapunov)

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

在jacobi法代码基础上,我们稍作调整也可以得qr法对应的lyapunov指数谱:

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

②引入qr函数

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

代码实现

 

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

 lyapunov指数谱:

 

 

举报

相关推荐

0 条评论