在进行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:nx(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:nx2=1-a(j)*x^2+y; y=0.3*x; x=x2; %首先进行序列迭代endif x(end)>-100 & x(end)<100 %若不发散，再计算指数x=0.2; y=0.3;JJ=eye(2);%2阶单位矩阵for i=1:nx2=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;endL=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:nx2=1-a(j)*x^2+y; y=0.3*x; x=x2; %首先进行序列迭代endif x(end)>-100 & x(end)<100 %若不发散，再计算指数x=0.2; y=0.3;JJ=eye(2);%2阶单位矩阵for i=1:nx2=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)));endendLE = 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指数谱：

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

定义指令时“控制器”，“链接”和“编译”函数之间的区别
本文翻译自:Difference between the 'controller', 'link' and 'compile' functions when defining a directive ...
lyapunov指数 matlab计算_Matlab学习笔记1——B站台大课
1.Matlab调出某一视窗菜单栏layout里面找 2.Matlab作为Calculator (1)ans 是计算结果 (2)计算规则:先*/再+-,合理使用(),从内往外写 (3)一些特殊输入 ...
在画电路图时，想问下几种地之间的区别？ power-GND singal-GND GND
在电路中,任何信号或者电源都是有回路的,这个回路的终点就是地. 电源地,一般是指在DCDC输入端的地,这个地有可能噪声比较大,或者别的情况不能和DCDC的输出端共地而区分的. 数字地模拟地故名思意,就 ...
人体运动学非线性分析（二）—Lyapunov指数
基于轨迹的Lyapunov指数计算原理及在人体运动学中的应用 Lyapunov基本概念 Lyapunov指数能够通过测量相邻轨迹收敛或分离的速率来衡量混沌的状态,Lyapunov指数越大,表明系统局部 ...
坐标转换参数求取的若干问题
摘要:不同坐标系下的测量成果需要相互统一,统一的方法是进行坐标转换,而坐标转换的实质是转换参数的求取.本文针对坐标转换中转换参数的初值计算问题.转换参数的精度问题及转换参数的求取方法三个方面进行讨论 ...
MySQL学习笔记（五）并发时经典常见的死锁原因及解决方法
MySQL学习笔记(五)并发时经典常见的死锁原因及解决方法参考文章: (1)MySQL学习笔记(五)并发时经典常见的死锁原因及解决方法 (2)https://www.cnblogs.com/tiny ...
图形处理（九）点云重建（下）法矢求取、有向距离场等值面提取
一.相关理论在点云重建算法中,点云法矢的求取是非常重要的一步.之前导师让我做点云尖锐特征边重建时,发了一堆异向法矢求取的英文paper,当时我很迷糊,就问了老师,让我做点云重建,为什么发的文献资料都 ...
matlab lyapunov指数,lyapunov指数matlab
(1)定义法 Lyapunov 指数的计算方法定义法求解 Lyapunov 指数.JPG 关于定义法求解的程序,与 matlab 板块的"连续系统 LE 求解程序"差不多.以 R ...
N维离散混沌映射Lyapunov指数的计算
本文代码均以Python实现首先给出按照定义法进行求解的程序,这里选取Henon映射为示例 #-*-coding:utf-8-*- ''' 多变量非线性方程求解 ''' from sympy imp ...

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

1.准备工作

1.1 henon映射

动力学系统

代码实现

1.2 jacobi矩阵

1.3 lyapunov指数

2.Jacobi法

代码实现

3.qr法

代码实现

lyapunov指数谱：

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

代码实现

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

最新文章

热门文章