Lyapunov函数的一些构造方法
-Lyapunov指数的计算方法
【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。
（1）定义法定义法求解Lyapunov指数.JPG关于定义法求解的程序，和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:iteratetimestspan = tstart:tstep:(tstart + 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.8924Wolf文章中计算得到的Rossler系统的LE为————0.09 0 -9.77需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
lyapunov函数
lyapunov函数Lyapunov函数是一种用于研究系统稳定性的重要工具，它可以用来检验系统耐受外界干扰的能力，以及系统发生振荡现象的可能性。
它是由俄罗斯数学家安德烈利亚普诺夫（Andrey Lyapunov）最早提出的，是一种重要的动态系统的稳定性理论。
Lyapunov函数也称为质能函数（或者拉普拉斯函数），它是一个定义在特定空间中的实值函数，它能够浓缩动态系统状态并同时反映系统的稳定性。
简而言之，Lyapunov函数有助于发现动态系统中某个状态点是否是稳定的。
Lyapunov函数的定义为：对于一个系统，它可以根据任意一个参数来评估系统状态点的稳定性。
当系统从状态点A移动到状态点B 时，系统的Lyapunov函数能够帮助我们了解更多的系统的行为：第一，该函数能够衡量不同状态点之间的影响，即该函数能够测量系统从A点变为B点的过程中，A点对B点的影响有多大。
第二，它能够表明系统能否在未来某个时刻稳定地维持本身的状态点。
由于Lyapunov函数可以衡量系统状态点的稳定性，因此可以用来实施控制策略，来防止系统振荡。
Lyapunov函数可以作为一个智能控制系统中的一个重要分量，它能够有效地检测和对应外部的环境因素，进而把外部的环境因素转换为控制指令，以便调节系统的状态。
Lyapunov函数在动态系统建模和分析方面，具有无可比拟的效用。
除了上述提到的用于检验系统稳定性以外，它还可以用来检查系统控制功能、建立系统模型和探索各种可能性，以此探究系统的行为特征，它对数学原理研究和应用研究都是非常重要的。
Lyapunov函数的研究仍在不断发展，越来越多的研究者将Lyapunov函数与其他技术融合在一起，以便更好地理解系统的行为，提高控制策略的可靠性，更好地探索异常情况的发生可能性，以及更多的分析细节。
总之，Lyapunov函数是一种重要的动态系统理论，它能够帮助我们检验系统的稳定性，为智能控制系统提供重要参考，在模型建立和控制策略制定等方面具有重要意义，而且它也在不断发展，以适应不断变化的环境和情况。
常微分方程中李雅普诺夫函数构造方法
武夷学院学报JOURNAL OF WUYI UNIVERSITY第40卷第3期2021年3月Vol.40 No.3Mar. 2021常微分方程中李雅普诺夫函数构造方法陆求赐】，张宋传2,王学彬2(1.武夷学院人文与教师教育学院，福建武夷山354300； 2.武夷学院数计学院，福建武夷山354300)扌商 要：李雅普诺夫直接法,就是在不求方程组解的情况下，构造一个李雅普诺夫函数,通过微分方程组所计算岀来的全导数的符号性质，来判断微分方程组零解的稳定性，但至今仍没有一种统一的方法来构造李雅普诺夫函数。
通过 对一个含参数例子的分析,介绍几种常见且适用的李雅普诺夫函数的构造方法:首次积分法、能量函数法、分离变量法、待定系数法和二次型矩阵法等。
关键词：常微分方程;李雅普诺夫函数;零解的稳定性;构造方法中图分类号：O175.26 文献标识码:A 考虑如下自治的微分方程组牯(x )(1)这里X =(x 1,x 2，…x ”)表示n ("逸2)维向量,i =1,2，…,n ,假设f (0)=0(0表示向量),且f (x )在某邻域G ：||x ||臆A (A 为正常数)内有连续的偏导数，从而方程组(1)的由其初值条件x (t °)=x 。
所决定的解在邻域G 内存在且唯一°李雅普诺夫第二方法:通过方程组(1)来构造一个特殊的函数V (x ),并假设V (x )关于所有变元的偏导数存在且连续，将V (x )对变量t 求导，并将方程组(1)代入得：d V § dV dx , § 坠卩 f(2)dt =移 dx, dt=移 dx ,)这个式(2)就是函数V (x )沿方程组(1)求得的全导数，现通过式(2)的符号来判定方程组(1)零解的稳定 性,这就是李雅普诺夫第二方法的思路,这个函数收稿日期：2020-07-21基金项目：武夷学院校科研基金项目(XL201408)；福建省教育厅科技项目(JA15512、JAT160519)；福建省自然科学基金(2016J01682)；高级引进人才 科研启动基金(YJ201802)资助°作者简介:陆求赐(1975-)，男,汉族,副教授，主要从事基础数学教学和微分方程的研究°文章编号：1674-2109(2021)03-0016-05V (x )就称为李雅普诺夫函数，简称V 函数［1-2］°李雅普诺夫函数作为常微分方程的重要内容之一,学好它可以为进入常微分方程领域的研究打好基础。
Lyapunov函数的一些构造方法
dL 并且 1 dt
y v 0 x x* , * * , w w* . y v M
OK
Lyapunov 函数的构造
三、Lyapunov函数形式的推广
* y u y u x* x*v* v u v* 其中 F * du * du du. * * x y v u u ky u x
推广 形式：
u w* w* u du,
w
w w*
u w* u w du, * du, w f u
*
w( w* ) du, f u
w
w*
f (u ) f ( w* ) du. u
Lyapunov 函数的构造
四、一类SI模型的全局稳定性
Lyapunov 函数的一些构造方法
李建全
空军工程大学(西安)
E-mail: jianq_li@
Lyapunov 函数的构造
一、引言 二、带有免疫的病毒动力学模型 三、Lyapunov函数形式的推广 四、一类SI模型的全局稳定性
Lyapunov 函数的构造
一、引言
全局渐近稳定性 (GAS)是动力系统研究的重要 问题。 对于高维自治微分系统 ，证明GAS最常用的方法是 Lyapunov - LaSalle定理， 即若自治系统 * x' f x, x Rn
函数 2、 情形f w 1. 定义Lyapunov
*
dL 则 1 dt
* x x p 0 * * x 2 w w * x x M c
2
x* xy*v y v* x v 3 x x* yv* y *v 0,
Liapunov函数的构造
Liapunov 函数的构造摘要：Liapunov 函数是一种判定微分方程零解稳定性的重要方法，所以本文首先介绍了Liapunov 函数以及判断微分方程的稳定性定理，然后着重介绍了Liapunov 函数的几种构造方法，包括常系数线性系统的巴尔巴欣公式、线性类比法.通过这两种构造方法，我们将初步了解Liapunov 函数的构造在判断微分方程零解稳定性中的重要作用. 关键词：Liapunov 函数；零解；稳定性.引言在常微分方程中，稳定性理论研究是很重要的一部分，即研究当时间趋于无穷时，其解的形态将会怎样变化，他在自然科学、工程力学、环境生态、社会经济等方面有着重要的应用。
在本章第一节中介绍了稳定性的相关定义，也介绍了对于可以求得微分方程的解析解时，如何利用定义判断其零解的稳定性。
但是在实际问题中提出的微分方程往往是很复杂的，无法求得其解析解，这就需要从方程本身来判断零解的稳定性Liapunov 直接方法就是求解这一问题的有效途径。
本文先引入Liapunov 函数，即V 函数的定义，以及Liapunov 稳定性的定理，然后介绍几种构造Liapunov 函数方法。
1 Liapunov 稳定性的定理1.1 V 函数设函数(x)V 在nR 中原点的某邻域U 中有定义,(x)V 在U 中连续可微，且满足(0)0V =定义 1.1若除原点外对所有x U ∈均有(x)0((x)0)V V ><，则称(x)V 为正定函数（负定函数）；若除原点外对所有x U ∈均有(x)0((x)0)V V ≥≤，则称(x)V 为半正定函数（半负定函数）；若在U 中原点的任一邻域内(x)V 既可以取正，也可以取负，则称(x)V 为变号函数.例如，222123(x)x x x V =++是3R 中的正定函数，但在4R 中确实半正定函数，而222123(x)x x x V =+-是3R 中的变号函数.一般(x)V 函数的符号判断十分困难，通常把(x)V 在原点展开为Taylor 级数12(x)V (x)V (x)V (x),m m m V ++=+++其中12V (x),V (x),V (x)m m m ++分C 是x 的m 次，1m +，2m +齐次函数，根据(x)V 展开式中的最低次项的系数，通常就可以判断(x)V 在原点邻域内的符号.因为再原点附近其他项都可以视为第一项的高阶无穷小. 1.2 Liapunov 稳定性定理设n 维自治微分方程(t)(x),(0)0dx f f dt==（1.1） 的解为12(t)(x (t),x (t),x (t))Tn x = ，为了研究方程（1.1）零解的稳定性，考察随时间变化时(x(t))V 的变化情况，将(x(t))V 视为t 的复合函数，关于t 求导可得11(x(t))(x)(x)n n k k k k k kdx dV V Vf V dt x dt x ==∂∂==∂∂∑∑ （1.2） 式（1.2）称为函数(x)V 沿着方程（1.1）轨线的全导数.介绍了(x)V 以及其全导数后，接下来简单介绍下Liapunov 稳定性理论的几个定理. 定理1.1若有原点的邻域U 和一个正定（负定）函数(x)V ，使得其全导数(x)V是半负定（半正定），则称系统（1.1）的零解是稳定的；特别地，当(x)V是负定（正定）时，系统（1.1）的零解是渐进稳定的.定理1.2设在原点的邻域U 内有函数(x)V ，它沿着方程（1.1）的轨线的全导数(x)V是正定（负定）的，而(x)V 本身不是半负定（半正定）的，则方程（1.1）的零解是不稳定的. 这两个定理是直接通过构造Liapunov 函数，来判断方程的零解是否稳定的，定理在本章第2节已经详细的证明过，这里不再做证明.2 Liapunov 函数的构造第一节中所讨论的两个定理都是一个函数稳定的充分条件，即存在一个()V x ，和它的全导数满足定理1.1时，系统的零解是稳定的. 满足定理1.2的条件时，系统的零解是不稳定的. 在使用Liapunov 函数判定稳定性时应当注意，当找不到满足稳定性定理的条件的函数()V x 时，并无法断言此系统的零解是不稳定的,并且构造的Liapunov 函数不同时，判断零解是否渐进稳定以及吸引域的大小也会有些差异.再利用Liapunov 方法判断系统零解稳定性时，需要明确满足一定条件的Liapunov 函数是否存在，即当系统的零解有某种稳定性时，满足这个稳定性定理的()V x 是否存在，这就是上述定理1.1和定理1.2的逆命题.是否成立. 2.1 Liapunov 函数的存在性考虑微分方程组d (t,x),(t,0)0xf f dt==（2.1）记(){}0,|,G t x t t x h =≥≤，设(t,)f x 在G 连续，关于x 满足Liapunov 条件.令00(t)(t,t ,)x x φ=是方程组（2.1）满足00(t )x φ=的解.定理2.1 若方程组2.1的零解是稳定的，则有正定函数(t,)V x ，使得其全导数(t,)V x是半负定的.证 首先根据方程组（2.1）的零解构造出正定函数(t,)V x ，在验证(t,)V x是半负定.取02(t t )0(t,)(1)(t ,t,)V x e x φ--=+020,()min(t ,t,)t t x y hW x y φ≥≤≤=其中0(t ,t,)x φ表示方程组（2.1）在t 时刻过x 的解在0t 时刻的位置坐标. 显然有02(t t )0(t,0)(1)(t ,t,0)0V e φ--=+=因为方程组（2.1）的零解是稳定的，所以对任意的0ε>，有0δ>，使得当x δ<时，有00(t,t ,),x t t φε<≥（2.2）于是当x h ε≤≤时，对所有的0t t ≥有0(t ,t,)0x φδ≥>.否则就有0x ，0x h ε≤≤，以及10t t ≥，使得010(t ,t ,)x φδ<.由式（2.2）得1001010(t ,t ,(t ,t ,x )),t t φφε<≥又因为1001001000(t ,t ,(t ,t ,x ))(t ,t ,x )x ,x φφφε==≥与式（2.2）矛盾. 由此得(x)W 是正定函数，显然有2010(t,x)(t ,t ,x )(x)V W φ≥≥即(t,x)V 是正定函数，另一方面002(t t )000002(t t )(t,(t,t ,x ))(1)(t ,t,(t,t ,x ))(1)V e ex φφφ----=+=+所以有02(t t )000(t,(t,t ,x ))0dV e x dtφ--=-<即(t,x)V是半负定的.例1研究下述微分方程组零解的稳定性.1221dx x dtdx xdt⎧=-⎪⎪⎨⎪=⎪⎩ 满足初值问题10102020(t ),(t )x x x x ==.解容易解得满足上述初值问题的解为101020100200201020100200(t,t ,,)cos(t t )sin(t t ),(t,t ,,)sin(t t )cos(t t ).x x x x x x x x φφ=---⎧⎨=-+-⎩ 所以有定理2.1，可以取函数0(t t )221210122012(t,,)(1e )((t ,t,,)(t ,t,,))V x x x x x x φφ--=++0(t t )2212(1e )()x x --=++显然0(t t )2222121212(t,,)(1e )()V x x x x x x --=++≥+即通过这种方法构造的12(t,,)V x x 是正定函数.并求得其全导数为：0(t t )221212(t,,)e ()V x x x x --=-+即12(t,,)V x x 正定，12(t,,)V x x半负定，所以上述方程组的零解是稳定的. 2.2 常系数线性系统的巴尔巴欣公式对于常系数线性微分方程组d ,(),,n ij n n xAx A a x R dt⨯==∈（2.3） 这里可以假设该线性自治系统的Liapunov 函数为一个二次型，不妨设其为()TV x x Bx =，则他沿系统（2.3）的全导数为()()(A B BA).TT T d V x x Bx x x dt==+显然其全导数依然是一个二次型.这样，可以从()V x的二次型出发，利用A B BA C T+=进而来确定二次型()V x ，再根据()V x 和()V x的符号来判断方程组（2.3）零解的稳定性.下面就利用这一思想介绍二阶方程组的巴尔巴欣公式.对于二维微分方程组11111222211222dx a x a x dtdx a x a xdt⎧=+⎪⎪⎨⎪=+⎪⎩（2.4）可以给出一个二次型221111212222()2W x w x w x x w x =++进而来求这样的二次型221111212222()2V x v x v x x v x =++使得()V x 沿着方程组（2.4）解曲线的全导数满足(2.4)()2()dV x W x dt =这里的()W x 是选取的，因此111222,,w w w 可视为已知量，利用()W x ，()V x 的定义和上述全导数的等式可以求得未知量111222,,v v v 满足下述方程组：111121121112111122122122121212222222()a v a v w a v a a v a v w a v a v w +=⎧⎪+++=⎨⎪+=⎩当112211221221()()0a a a a a a ∆=+-≠时，可以求解出111222,,v v v ，即可得出()V x2211221111211212112221221222201()20x x x x w a a V x w a a a a w a a =-∆+通过这个过程，可知构造出了一个Liapunov 函数()V x ，并求出其沿着自治系统（2.4）解曲线的全导数()2()V x W x =，从而可以根据()V x 和()W x 的符号来判断系统（2.4）零解的稳定性.例2 讨论微分方程组1121124426dx x x dtdx x xdt⎧=-+⎪⎪⎨⎪=-⎪⎩ 解取221212(,)()W x x x x =-+，利用上述的巴尔巴欣公式可得2211220214201()16004102146x x x x V x --=---()222221122121211(786)(32)2020x x x x x x x x =++==+++ 显然()V x 是正定的，()2()V x W x =是负定的，故上述方程组的零解是渐进稳定的. 2.3 线性类比法线性类比法是将一些非线性系统当作线性系统，用类比的方法构造出需要的Liapunov 函数.例3设1()f x 连续可导，(0)0f =，讨论微分方程组111222211222()dx f x a x dtdx a x a xdt⎧=+⎪⎪⎨⎪=+⎪⎩（2.5） 零解的稳定性.解当1111()f x a x =时，上述非线性系统（2.5）就是线性系统11111222211222dx a x a x dtdx a x a xdt⎧=+⎪⎪⎨⎪=+⎪⎩（2.6） 线性系统（2.6）的特征方程为2112211221221()0a a a a a a λλ-++-=容易看出，当1122112212210,0a a a a a a +<->时，特征方程有两个根都有负实部，线性系统（2.6）的零解是渐进稳定的.非线性系统（2.5）的零解稳定性无法用特征根的方法来判断，但可以用类似于线性系统的Liapunov 函数去判断其稳定性.事实上，对线性系统（2.6）取2121122112212211(,)()()V x x a a a a a a x =+-则12(,)V x x是半负定的函数，利用巴尔巴欣公式得221211221221122112211(,)()()22V x x a a a a x a x a x =-+-12(,)V x x 是正定函数，所以线性系统（2.6）的零解是稳定的.对比非线性系统（2.5）和线性系统（2.6），线性系统（2.6）中的11a 相当于非线性系统（2.5）中的11()f x x .所以很自然地猜想，当112222122111()()0,0f x f x a a a a x x +<->（2.7）成立时，非线性系统（2.5）的零解很可能是稳定的，注意到12112212211112212211101()()2x a a a a x a a a a x dx -=-⎰由此类比构造与线性系统类似的V 函数121122212211122112201()1(,)()()2x f x V x x a a a x dx a x a x x =-+-⎰ 当式（2.7）的条件成立时，12(,)V x x 正定，计算全导数得12221122111122(2.5)(,)(())(())V x x a f x a a x f x a x =-++22112222112212211222()((())())a x a x a f x a x a a x a x -+-+21122122122111()()()()0f x f x a a a a x x x =-+≤ 所以，非线性系统（2.5）的零解是稳定的.3 总结本文中我们仅仅给出了两种构造Liapunov 函数的方法，也很容易看出即便是这两种发放也都有其局限性. 仅仅对于常系数线性系统才可以使用巴尔巴欣公式，并且当微分方程组的阶数太大时，利用巴尔巴欣公式将会产生巨大的计算量.线性类比法需要我们先找到非线性系统对应的线性系统的Liapunov 函数，然后类比构造出非线性系统的Liapunov 函数，但是在实际问题中，往往很难进行类比，更多的要靠个人的经验总结.在教材上还介绍了能量函数法、分离变量法、变梯度法，微分方程组的形式结构不同，就要用不同的方法，遗憾的是至今仍未形成构造Liapunov 函数的通用方法，找到通用方法最大的困难是依赖于一个未知的V 函数的存在，然而构造这种V 函数又没有一般的规律可循，这是一种技巧性问题，更大的可能只能靠个人的经验.但是，我相信随着学科的发展，数学家们会继续研究构造Liapunov 函数的一般方法，Liapunov 函数的构造定会取得重大突破.4 参考文献[1]马知恩，周义仓，李承治. 常微分方程定性与稳定性方法[M]. 北京:科学出版社,2015.6. [2]张志芬，丁同仁，黄文灶，董镇喜. 微分方程定性理论[M]. 北京: 科学出版社,1985. [3]王高雄，周之铭，朱思铭等. 常微分方程[M]. 北京：高等教育出版社,2006.7.。
离散周期lyapunov方程和离散周期riccati方程的迭代算法
摘要作为线性时变系统的最简单形式，线性周期系统由于其广泛的应用，一直是学者们研究的热点。
线性周期系统，是一类系数矩阵带有周期性的线性系统，在各个领域中都有着广泛的应用。
为了研究离散周期系统的稳定性问题，离散周期Lyapunov方程的求解就显得至关重要。
同样，在进行离散周期系统的线性二次最优状态反馈控制器的设计时，需要用到离散周期Riccati方程的解。
基于这样的研究背景，本文针对离散周期系统下的Lyapunov方程和Riccati方程，给出了其求解的迭代算法。
针对离散周期Lyapunov方程，推导出了相应的迭代算法，分别对零初始条件和任意初始条件的情况给出了严谨的收敛性证明，并通过数值仿真验证了算法的有效性。
并且将最新估计信息的思想引入了迭代算法，得到了新的基于最新估计信息的迭代算法，同样对给出了算法在零初始条件下和非零初始条件下，迭代算法的严谨的收敛性证明，利用数值仿真例子证明了算法是有效并且收敛的。
并且通过对两种算法的数值仿真对比发现，基于最新估计信息的迭代算法的收敛速度要快于原始的迭代算法，从而验证了加入最新估计信息的迭代算法的优越性。
针对推导出的离散周期Riccati方程的迭代算法，给出了其在零初始条件下的收敛性证明，并通过数值仿真验证了算法的有效性，同样，为了改进算法，加入了最新估计信息，得到了新的基于最新估计信息的迭代算法。
同样对该算法的收敛性进行了严谨的证明与数值仿真验证，说明了该算法是有效可用的。
针对两种方程的迭代算法，为了研究最新估计信息对迭代算法的影响程度，引入了加权的思想，得到了带权重因子的新的迭代算法，并进行了收敛性证明。
通过数值仿真，给出了不同权重因子下的收敛性曲线，通过对比可以看出当全部使用最新估计信息时，算法的收敛速度最快，由此可见，加入最新估计信息能有效提高迭代算法的收敛速度。
关键词：离散周期系统；Lyapunov方程；Riccati方程；迭代算法AbstractAs the simplest form of time-varying linear systems, periodic linear systems have been attracting much attention during the past several decades. This is partially because this type of systems has very wide application. To investigate the stabilization problem of the periodic linear systems, it is important to achieve the solution of the periodic Lyapunov matrix equation. Similarly, the design of linear quadratic optimal state feedback controller based on the robust control is related to the stabilizing positive definite solution of Riccati equation. Based on this research background, we propose iterative algorithms for solving discrete-time periodic Lyapunov matrix equation and discrete-time periodic Riccati matrix equation.Iterative algorithms for discrete periodic Lyapunov equations are derived, respectively to the zero initial conditions and arbitrary initial conditions. And the proof of convergence is given. The effectiveness of the algorithm is verified by numerical simulation. And the latest information estimation theory is into the iterative algorithm, the proof of the convergence is also given. The validity of the algorithm is verified by numerical simulations. Finally, the simulation analysis of the two algorithms find that the convergence rate of the iterative algorithm based on the estimation of the latest information is faster than the original algorithm. It proves the superiority of the iterative algorithm adding the latest information of the estimation.Iterative algorithm for discrete periodic Riccati equations is derived, given the zero initial condition of convergence, and the effectiveness of the algorithm is verified through numerical simulation. In order to improve the algorithm with the latest estimate information, a new iterative algorithm based on the information of the latest estimation is given. The convergence of the new algorithm is proved and the validity of the algorithm is verified by numerical simulation. Through numerical simulation, the convergence curves of different weighting factors are given. It found that using the latest estimate information, the convergence speed is the fastest. Therefore, adding the latest estimation information can effectively improve the convergence speed of iterative algorithm.Key words：discrete-time linear periodic system，periodic Lyapunov equations，periodic Riccati equations，iterative algorithms目录摘要 (I)ABSTRACT ..................................................................................................................... I I 第1章绪论 . (1)1.1课题的来源及研究的背景意义 (1)1.2国内外在该方向上的研究现状及分析 (2)1.3本文的主要研究内容 (6)第2章离散周期系统Lyapunov方程快速迭代算法 (8)2.1相关的概念与性质 (8)2.2原始迭代算法 (9)2.2.1显式迭代算法 (9)2.2.2数值仿真 (12)2.3基于最新估计信息的迭代算法 (16)2.3.1显示迭代算法 (16)2.3.2数值仿真 (19)2.4本章小结 (24)第3章离散周期Riccati方程的迭代算法 (25)3.1相关的概念与性质 (25)3.2问题的描述 (25)3.3原始迭代算法 (25)3.3.1显示迭代算法 (26)3.3.2数值仿真 (28)3.4基于最新估计信息的迭代算法 (29)3.4.1显示迭代算法 (30)3.4.2数值仿真 (32)3.5本章小结 (34)第4章离散周期Riccati方程的加权最新估计迭代算法 (35)4.1 带加权因子的快速迭代算法 (35)4.2数值仿真 (37)4.3本章小结 (39)结论 (40)参考文献 (41) (45)致谢 (46)第1章绪论1.1课题的来源及研究的背景意义随着对控制系统的研究越来越深入，人们发现，许多生活中的系统是线性周期系统。
lyapunov-schmidt方法
Lyapunov-Schmidt方法是一种用于非线性方程组的求解的数值方法。
它是由俄罗斯数学家Aleksandr Lyapunov和德国数学家Ernst Schmidt分别在19世纪和20世纪提出的。
这种方法在处理非线性问题时非常有效，并且在应用数学和工程领域得到了广泛的应用。
Lyapunov-Schmidt方法的核心思想是将原始的非线性方程组转化成一系列线性方程组，从而简化求解过程。
这种方法的优势在于可以通过有限步骤来逼近非线性方程组的解，从而大大提高了求解效率。
下面我们将详细介绍Lyapunov-Schmidt方法的原理和应用。
1. Lyapunov-Schmidt方法的原理Lyapunov-Schmidt方法的原理是通过引入一组正交归一的特征函数，将原始的非线性方程组转化为一系列正交归一的线性方程组。
这样一来，原始的非线性方程组就被分解成了一系列互相独立的线性方程组，从而使得求解过程变得更加简单和高效。
2. Lyapunov-Schmidt方法的应用Lyapunov-Schmidt方法在科学和工程领域有着广泛的应用。
比如在物理学中，通过Lyapunov-Schmidt方法可以求解复杂的非线性波动方程，从而对物质的运动和变形进行研究。
在工程领域，Lyapunov-Schmidt方法可以用于求解具有非线性特性的结构力学问题，如弹性体的变形和弹性波的传播等。
3. 使用案例我们以一个简单的非线性方程组为例来说明Lyapunov-Schmidt方法的求解过程。
假设我们有一个非线性方程组：f(x, y) = 0g(x, y) = 0我们希望求解这个方程组的解。
我们可以通过Lyapunov-Schmidt方法将原始的非线性方程组转化为一系列正交归一的线性方程组：Φ1(x, y) = 0λ1(Φ1x + Φ1y) = 0Φ2(x, y) = 0λ2(Φ2x + Φ2y) = 0...我们可以通过求解这一系列线性方程组来逼近原始的非线性方程组的解。
lyapunov指数计算
lyapunov指数计算
Lyapunov指数（Lyapunov exponent）是一种用于描述动态系统混沌性质的指标。
在数学上，Lyapunov指数是描述线性化系统的稳定性的指标，它可以判断非线性系统是否具有混沌性质。
计算Lyapunov指数的基本过程如下：
1.首先，选择一个合适的初始状态，并计算该状态在系统中的轨迹。
2.然后，选取一个邻域，计算在该邻域内的状态与初始状态的差异随时间的变化情况。
3.对于每个时间步长，计算邻域中的点向初始状态点移动的距离与时间的比值。
4. 重复上述步骤，直到获得足够的数据，然后计算Lyapunov指数。
5. Lyapunov指数表示在该系统中相邻轨迹的指数级别分离速度。
具体计算Lyapunov指数的过程比较复杂，一般需要借助计算机进行模拟和计算。
-Lyapunov指数的计算方法
【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚，看了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线，把目前已有的LE计算方法做一个汇总！1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法，或者通过求解系统的微分方程，得到微分方程解的时间序列，然后利用时间序列（即离散系统）的LE求解方法来计算得到。
关于连续系统LE的计算，主要以定义方法、Jacobian方法做主要介绍内容。
（1）定义法定义法求解Lyapunov指数.JPG关于定义法求解的程序，和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:iteratetimestspan = tstart:tstep:(tstart + 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.8924Wolf文章中计算得到的Rossler系统的LE为————0.09 0 -9.77需要注意的是——定义法求解的精度有限，对有些系统的计算往往出现计果和理论值有偏差的现象。
高维传染病模型的Lyapunov函数构造
栾培译，李 静 临沂大学数学与统计学院，山东 临沂
收稿日期：2020年11月21日；录用日期：2020年12月20日；发布日期：2020年12月28日
摘要
本文主要研究几类传染病模型如SIR，SIRS，SIS和SEIR模型的Lypunov函数构造方法，从而获得传染病 模型全局稳定性的结论。
栾培译，李静
Copyright © 2020 by author(s) and Hans Publishers Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). /licenses/by/4.0/
Received: Nov. 21st, 2020; accepted: Dec. 20th, 2020; published: Dec. 28th, 2020
Abstract
Lypunov functions for classical epidemiological models are introduced such as SIR, SIRS, SIS and SEIR. Global stability of some epidemiological models is also established.
以及
∂2V ∂S 2 ∂2V ∂S ∂I
∂2V ∂S ∂I= ∂2V ∂I 2
∂2V ∂2V ∂2V
∂S
2
∂I 2
−
∂S ∂I
=
1 Nγ S2I2 β > 0
很容易看出 E∗ 为唯一极值状态且为函数V (S, I ) 在 R+2 内的最小值。
Lyapunov指数的计算方法
【总结】Lyapunov指数的计算方法非线性理论近期为了把计算LE的一些问题弄清楚,瞧了有7～9本书！下面以吕金虎《混沌时间序列分析及其应用》、马军海《复杂非线性系统的重构技术》为主线,把目前已有的LE计算方法做一个汇总！1、关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法,或者通过求解系统的微分方程,得到微分方程解的时间序列,然后利用时间序列(即离散系统)的LE求解方法来计算得到。
关于连续系统LE的计算,主要以定义方法、Jacobian方法做主要介绍内容。
(1)定义法定义法求解Lyapunov指数、JPG关于定义法求解的程序,与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:iteratetimestspan = tstart:tstep:(tstart + 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、8924Wolf文章中计算得到的Rossler系统的LE为————0、09 0 -9、77需要注意的就是——定义法求解的精度有限,对有些系统的计算往往出现计果与理论值有偏差的现象。
Lyapunov稳定性理论概述
一， 稳定性的概念
初始值的微分变化对不同系统的影响不同，例如初始值问题
dx = ax ， x(0)=x0 , t≥0，x0≥0
(1)
dt
x e 的解为 x(t) = 0 at ，而x=0 是(1)式的一个解。当a f 0时，无论|x0|多小，只要
|x0| ≠ 0 ,在t→+∞时，总有x（t）→ ∞，即初始值的微小变化会导致解的误
x e 差任意大，而当a ≺0时， x(t) = 0 at 。与零解的误差不会超过初始误差x0，且随
着t 值的增加很快就会消失，所以，当|x0|很小时，x(t)与零解的误差也很小。
这个例子表明a f 0时的零解是“稳定”的。下面，我们就给出微分方程零解稳
定的严格定义。
设微分方程
R dx
dt
=
f
(t, x) ，
的解) 正定(>0) 半正定(≥0)且不恒为0 (对任意非零的初始状态
的解)
结论 该平衡态渐近稳定
该平衡态渐近稳定
ห้องสมุดไป่ตู้该平衡态稳定 但非渐近稳定
该平衡态不稳定
该平衡态不稳定
经过艰苦的研究证明，学者们发现，在上述三种定理中，只有Lyapunov的 渐近稳定性定理不可逆，其他定理，包括推广的一致稳定、一致渐近稳定、指数 稳定、全局指数稳定及不稳定定理等所有定理，都是可逆的。
数稳定，则可以任意给定负定矩阵-C，作 V = xT B x，其中B为线性矩阵不等式
BA+ATB=-C的解。这是根据上述方法2的思想所做出的构造过程。
四， Lyapunov方法的发展
世界著名数学大师Hirsch和Smale在他们的专著《常微分方程·动力系统·线
性代数》的序言中谈到：“有人说常微分方程这一学科是求解技巧和提示的汇集，
lyap函数
lyap函数LYAP函数是指Lyapunov函数，其名称来源于俄罗斯数学家亚历山大·米哈伊洛维奇·李亚普诺夫。
这种函数常常用于控制系统中，用来确定系统的稳定性。
Lyapunov函数类似于能量函数，它描述了系统的动态性能，能够以一种数学方式证明系统是否稳定。
在控制系统的设计中，LYAP函数是一种重要的思想工具，可以帮助工程师分析、评估和优化系统的性能。
Lyapunov函数在数学上的定义是：一个实数函数V(x)被称为系统状态x的Lyapunov 函数，当且仅当V(x)>0，且dV(x)/dt<0。
其中，dV(x)/dt表示V(x)的时间导数，即V(x)对时间的变化率。
这个定义的意义是，Lyapunov函数V(x)必须是正定义的，也就是说，V(x)只有当x=0时取到最小值；而且它的时间导数必须是负的，也就是说，V(x)随时间的演化将不断减小。
这个定义可以用来证明系统的稳定性。
如果在系统状态x处满足这个定义，那么这个状态就是一个稳态点。
系统中的任何轨迹都会朝着这个稳态点收敛。
Lyapunov函数的构造并不总是很容易，因为它必须同时满足正定和时间导数为负这两个限制条件。
如果有多个稳态点，还必须为每个稳态点构造一个相应的Lyapunov函数。
一般情况下，Lyapunov函数的形式是非常复杂的，但掌握这个方法的关键是应用恰当的技巧和方法进行构造。
可以用一些特殊的技巧来寻找Lyapunov函数，比如拉格朗日乘数法和优化理论。
Lyapunov函数例子吸引-排斥Lyapunov函数为它衡量一个非线性动力系统的稳定性。
它的形式通常如下：V(x, y) = -x^2-y^2+a*ln(x^2+y^2)其中x、y分别表示系统状态的两个维度，a是一个正数。
Lyapunov函数中的第一项-x^2-y^2表示系统状态的坐标离原点的距离的平方，它越小说明距离原点越近，系统越稳定。
第二项a*ln(x^2+y^2)描述了系统状态的坐标距离原点的距离变化率，它可以看做是一个动能项。
Lyapunov指数的计算方法
总结Lyapunov 指数的计算方法近期为了把计算LE的一些问题弄清楚,看了有7～9本书下面以吕金虎混沌时间序列分析及其应用、马军海复杂非线性系统的重构技术为主线,把目前已有的LE计算方法做一个汇总1. 关于连续系统Lyapunov指数的计算方法连续系统LE的计算方法主要有定义方法、Jacobian方法、QR分解方法、奇异值分解方法,或者通过求解系统的微分方程,得到微分方程解的时间序列,然后利用时间序列即离散系统的LE求解方法来计算得到;关于连续系统LE的计算,主要以定义方法、Jacobian方法做主要介绍内容;1定义法关于定义法求解的程序,和matlab板块的“连续系统LE求解程序”差不多;以Rossler 系统为例Rossler系统微分方程定义程序function dX = Rossler_lyt,X% Rossler吸引子,用来计算Lyapunov指数% a=,b=,c=% dx/dt = -y-z,% dy/dt = x+ay,% dz/dt = b+zx-c,a = ;b = ;c = ;x=X1; y=X2; z=X3;% Y的三个列向量为相互正交的单位向量Y = X4, X7, X10;X5, X8, X11;X6, X9, X12;% 输出向量的初始化,必不可少dX = zeros12,1;% Rossler吸引子dX1 = -y-z;dX2 = x+ay;dX3 = b+zx-c;% Rossler吸引子的Jacobi矩阵Jaco = 0 -1 -1;1 a 0;z 0 x-c;dX4:12 = JacoY;求解LE代码：% 计算Rossler吸引子的Lyapunov指数clear;yinit = 1,1,1;orthmatrix = 1 0 0;0 1 0;0 0 1;a = ;b = ;c = ;y = zeros12,1;% 初始化输入y1:3 = yinit;y4:12 = orthmatrix;tstart = 0; % 时间初始值tstep = 1e-3; % 时间步长wholetimes = 1e5; % 总的循环次数steps = 10; % 每次演化的步数iteratetimes = wholetimes/steps; % 演化的次数mod = zeros3,1;lp = zeros3,1;% 初始化三个Lyapunov指数Lyapunov1 = zerositeratetimes,1;Lyapunov2 = zerositeratetimes,1;Lyapunov3 = zerositeratetimes,1;for i=1:iteratetimestspan = tstart:tstep:tstart + tstepsteps;T,Y = ode45'Rossler_ly', tspan, y;% 取积分得到的最后一个时刻的值y = YsizeY,1,:;% 重新定义起始时刻tstart = tstart + tstepsteps;y0 = y4 y7 y10;y5 y8 y11;y6 y9 y12;%正交化y0 = ThreeGSy0;% 取三个向量的模mod1 = sqrty0:,1'y0:,1;mod2 = sqrty0:,2'y0:,2;mod3 = sqrty0:,3'y0:,3;y0:,1 = y0:,1/mod1;y0:,2 = y0:,2/mod2;y0:,3 = y0:,3/mod3;lp = lp+logabsmod;%三个Lyapunov指数Lyapunov1i = lp1/tstart;Lyapunov2i = lp2/tstart;Lyapunov3i = lp3/tstart;y4:12 = y0';end% 作Lyapunov指数谱图i = 1:iteratetimes;ploti,Lyapunov1,i,Lyapunov2,i,Lyapunov3程序中用到的ThreeGS程序如下：%G-S正交化function A = ThreeGSV % V 为33向量v1 = V:,1;v2 = V:,2;v3 = V:,3;a1 = zeros3,1;a2 = zeros3,1;a3 = zeros3,1;a1 = v1;a2 = v2-a1'v2/a1'a1a1;a3 = v3-a1'v3/a1'a1a1-a2'v3/a2'a2a2;A = a1,a2,a3;计算得到的Rossler系统的LE为————Wolf文章中计算得到的Rossler系统的LE为———— 0需要注意的是——定义法求解的精度有限,对有些系统的计算往往出现计果和理论值有偏差的现象;正交化程序可以根据上面的扩展到NN向量,这里就不加以说明了,对matlab用户来说应该还是比较简单的2Jacobian方法通过资料检索,发现论坛中用的较多的LET工具箱的算法原理就是Jacobian方法;基本原理就是首先求解出连续系统微分方程的近似解,然后对系统的Jacobian矩阵进行QR分解,计算Jacobian矩阵特征值的乘积,最后计算出LE和分数维;经过计算也证明了这种方法精度较高,对目前常见的混沌系统,如Lorenz、Henon、Duffing等的Lyapunov指数的计算精度都很高,而且程序编写有一定的规范,个人很推荐使用;虽然我自己要做的系统并不适用LET工具箱可以在网络上找到,这里就不列出了关于LET工具箱如果有问题,欢迎加入本帖讨论Jacobian法求解Lyapunov指数.JPG对离散动力系统,或者说是非线性时间序列,往往不需要计算出所有的Lyapunov指数,通常只需计算出其最大的Lyapunov指数即可;“1983年,格里波基证明了只要最大Lyapunov指数大于零,就可以肯定混沌的存在”;目前常用的计算混沌序列最大Lyapunov指数的方法主要有一下几种：1由定义法延伸的Nicolis方法2Jacobian方法3Wolf方法4P－范数方法5小数据量方法其中以Wolf方法和小数据量方法应用最为广泛,也最为普遍;下面对Nicolis方法、Wolf方法以及小数据量方法作一一介绍;1Nicolis方法这种方法和连续系统的定义方法类似,而且目前应用很有限制,因此只对其理论进行介绍,编程应用方面就省略了Nicolis方法求最大2Wolf方法Wolf方法求最大Wolf方法的Matlab程序如下：function lambda_1=lyapunov_wolfdata,N,m,tau,P% 该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法% m: 嵌入维数一般选大于等于10% tau:时间延迟一般选与周期相当,如我选2000% data:时间序列可以选1000；% N:时间序列长度满足公式：M=N-m-1tau=24000-10-11000=5000% P:时间序列的平均周期,选择演化相点距当前点的位置差,即若当前相点为I,则演化相点只能在|I－J|>P的相点中搜寻P=周期=600% lambda_1: 返回最大lyapunov指数值min_point=1 ; %&&要求最少搜索到的点数MAX_CISHU=5 ; %&&最大增加搜索范围次数%FLYINGHAWK% 求最大、最小和平均相点距离max_d = 0; %最大相点距离min_d = +100; %最小相点距离avg_dd = 0;Y=reconstitutiondata,N,m,tau; %相空间重构可将此段程序加到整个程序中,在时间循环内,可以保存时间序列的地方;见完整程序;M=N-m-1tau; %重构相空间中相点的个数for i = 1 : M-1for j = i+1 : Md = 0;for k = 1 : md = d + Yk,i-Yk,jYk,i-Yk,j;endd = sqrtd;if max_d < dmax_d = d;endif min_d > dmin_d = d;endavg_dd = avg_dd + d;endendavg_d = 2avg_dd/MM-1; %平均相点距离dlt_eps = avg_d - min_d ; %若在min_eps～max_eps中找不到演化相点时,对max_eps的放宽幅度min_eps = min_d + dlt_eps / 2 ; %演化相点与当前相点距离的最小限max_eps = min_d + 2 dlt_eps ; %&&演化相点与当前相点距离的最大限% 从P+1～M-1个相点中找与第一个相点最近的相点位置Loc_DK及其最短距离DKDK = +100; %第i个相点到其最近距离点的距离Loc_DK = 2; %第i个相点对应的最近距离点的下标for i = P+1:M-1 %限制短暂分离,从点P+1开始搜索d = 0;for k = 1 : md = d + Yk,i-Yk,1Yk,i-Yk,1;endd = sqrtd;if d < DK & d > min_epsDK = d;Loc_DK = i;endend% 以下计算各相点对应的李氏数保存到lmd数组中% i 为相点序号,从1到M-1,也是i-1点的演化点；Loc_DK为相点i-1对应最短距离的相点位置,DK为其对应的最短距离% Loc_DK+1为Loc_DK的演化点,DK1为i点到Loc_DK+1点的距离,称为演化距离% 前i个log2DK1/DK的累计和用于求i点的lambda值sum_lmd = 0 ; % 存放前i个log2DK1/DK的累计和for i = 2 : M-1 % 计算演化距离DK1 = 0;for k = 1 : mDK1 = DK1 + Yk,i-Yk,Loc_DK+1Yk,i-Yk,Loc_DK+1;endDK1 = sqrtDK1;old_Loc_DK = Loc_DK ; % 保存原最近位置相点old_DK=DK;% 计算前i个log2DK1/DK的累计和以及保存i点的李氏指数if DK1 ~= 0& DK ~= 0sum_lmd = sum_lmd + logDK1/DK /log2;endlmdi-1 = sum_lmd/i-1; 此处可以保存不同相点i对应的李氏指数,见完整程序;; % 以下寻找i点的最短距离:要求距离在指定距离范围内尽量短,与DK1的角度最小 point_num = 0; % &&在指定距离范围内找到的候选相点的个数cos_sita = 0 ; %&&夹角余弦的比较初值——要求一定是锐角zjfwcs=0 ; %&&增加范围次数while point_num == 0% 搜索相点for j = 1 : M-1if absj-i <=P-1 %&&候选点距当前点太近,跳过continue;end%计算候选点与当前点的距离dnew = 0;for k = 1 : mdnew = dnew + Yk,i-Yk,jYk,i-Yk,j;enddnew = sqrtdnew;if dnew < min_eps| dnew > max_eps %&&不在距离范围,跳过continue;end%计算夹角余弦及比较DOT = 0;for k = 1 : mDOT = DOT+Yk,i-Yk,jYk,i-Yk,old_Loc_DK+1;endCTH = DOT/dnewDK1;if acosCTH > 3./4 %&&不是小于45度的角,跳过continue;endif CTH > cos_sita %&&新夹角小于过去已找到的相点的夹角,保留cos_sita = CTH;Loc_DK = j;DK = dnew;endpoint_num = point_num +1;endif point_num <= min_pointmax_eps = max_eps + dlt_eps;zjfwcs =zjfwcs +1;if zjfwcs > MAX_CISHU %&&超过最大放宽次数,改找最近的点DK = +100;for ii = 1 : M-1if absi-ii <= P-1 %&&候选点距当前点太近,跳过continue;endd = 0;for k = 1 : md = d + Yk,i-Yk,iiYk,i-Yk,ii;endd = sqrtd;if d < DK & d > min_epsDK = d;Loc_DK = ii;endendbreak;endpoint_num = 0 ; %&&扩大距离范围后重新搜索cos_sita = 0;endendend%取平均得到最大李雅普诺夫指数此处只有一个值,若为正说明体系是混沌的,若为负则说明体系是周期性的确定性运动lambda_1=sumlmd/lengthlmd;程序中用到的reconstitution函数如下：此段程序可直接放在时间循环内部,即可以保存时间序列的地方;见完整程序范例;function X=reconstitutiondata,N,m,tau%该函数用来重构相空间% m 为嵌入空间维数% tau 为时间延迟% data 为输入时间序列% N 为时间序列长度% X 为输出,是mn维矩阵M=N-m-1tau; %相空间中点的个数for j=1:M %相空间重构for i=1:mXi,j=datai-1tau+j;endend这里声明一下,这些程序并非我自己编写的,均是转载,其使用我已经验证过,绝对可以运行3小数据量方法说小数据量方法是目前最实用、应用最广泛的方法应该不为过吧,呵呵小数据量方法求最大Lyapunov指数.JPG上面两种方法,即Wolf方法和小数据量方法,在计算LE之前,都要求对时间序列进行重构相空间,重构相空间的优良对于最大LE的计算精度影响非常大因此重构相空间的几个参数的确定就非常重要;1时间延迟主要推荐两种方法——自相关函数法、C－C方法自相关函数法——对一个混沌时间序列,可以先写出其自相关函数,然后作出自相关函数关于时间t的函数图像;根据数值试验结果,当自相关函数下降到初始值的1－1/e时,所得的时间t即为重构相空间的时间延迟;C－C方法——可以同时计算出时间延迟和时间窗口,个人推荐使用这种方法2平均周期平均周期的计算可以采用FFT方法;在matlab帮助中有一个太阳黑子的例子,现摘录如下：load %装载数据文件year = sunspot:,1; %读取年份信息wolfer = sunspot:,2; %读取黑子活动数据plotyear,wolfer %绘制原始数据图title'Sunspot Data'Y = fftwolfer; %快速FFT变换N = lengthY; %FFT变换后数据长度Y1 = ; %去掉Y的第一个数据,它是所有数据的和power = absY1:N/2.^2; %求功率谱nyquist = 1/2;freq = 1:N/2/N/2nyquist; %求频率plotfreq,power, grid on %绘制功率谱图xlabel'cycles/year'title'Periodogram'period = 1./freq; %年份周期plotperiod,power, axis0 40 0 2e7, grid on ％绘制年份－功率谱曲线ylabel'Power'xlabel'PeriodYears/Cycle'mp,index = maxpower; %求最高谱线所对应的年份下标periodindex %由下标求出平均周期3嵌入维数目前嵌入维数的主要计算方法是采用Grassberger和Procaccia提出的G－P算法计算出序列的关联维数d,然后利用嵌入维数m>=2d+1,选取合适的嵌入维数;G—P算法程序如下：function ln_r,ln_C=G_Pdata,N,tau,min_m,max_m,ss% the function is used to calculate correlation dimention with G-P algorithm% 计算关联维数的G－P算法% data:the time series 时间序列% N: the length of the time series 时间序列长度% tau: the time delay 时间延迟% min_m:the least embedded dimention m 最小的嵌入维数% max_m:the largest embedded dimention m 最大的嵌入维数% ss:the stepsize of r r的步长%skyhawkfor m=min_m:max_mY=reconstitutiondata,N,m,tau;%reconstitute state spaceM=N-m-1tau;%the number of points in state spacefor i=1:M-1for j=i+1:Mdi,j=maxabsY:,i-Y:,j;%calculate the distance of each twoend %points in state space 计算状态空间中每两点之间的距离endmax_d=maxmaxd;%the max distance of all points 得到所有点之间的最大距离d1,1=max_d;min_d=minmind;%the min distance of all points 得到所有点间的最短距离 delt=max_d-min_d/ss;%the stepsize of r 得到r的步长for k=1:ssr=min_d+kdelt;Ck=correlation_integralY,M,r;%calculate the correlation integralln_Cm,k=logCk;%lnCrln_rm,k=logr;%lnrfprintf'%d/%d/%d/%d\n',k,ss,m,max_m;endplotln_rm,:,ln_Cm,:;hold on;endfid=fopen'','w';fprintffid,'% %\n',ln_r;fclosefid;fid = fopen'','w';fprintffid,'% %\n',ln_C;fclosefid;程序中的correlation_integral函数如下：function C_I=correlation_integralX,M,r%the function is used to calculate correlation integral%C_I:the value of the correlation integral%X:the reconstituted state space,M is a mM matrix%m:the embedding demention%M:M is the number of embedded points in m-dimensional sapce%r:the radius of the Heaviside function,sigma/2<r<2sigma%calculate the sum of all the values of Heaviside%skyhawksum_H=0;for i=1:M% fprintf'%d/%d\n',i,M;for j=i+1:Md=normX:,i-X:,j,inf;%calculat the distances of each two points in matris M with sup-normsita=heavisider,d;%calculate the value of the heaviside functionsum_H=sum_H+sita;endendC_I=2sum_H/MM-1;%the value of correlation integral发布于2007-08-03 20:35:01感谢octopussheng的总结这是现有的方法的一个总结,不知道你对这些方法有些什么体会,或者说他们的局限,现在还有作这方面的改进的吗先总结总结,具体的比较会逐步贴出来的,而且这也需要大家的一起努力啊光是我在这里说,也只是我自己的一点体会,希望用过这些方法的一起来参与这个总结工作以上的各种方法在实际应用的时候要根据具体情况来选择;一般地,如果已知系统方程当然系统不能太过复杂时,则计算Lyapunov指数采用定义法、Jacobian方法要精确、简单些而如果系统方程比较复杂如超维系统、或者为一时间序列时,则推荐采样Wolf 方法、小数据量方法;Wolf方法的特点是时间序列无噪声,空间中小向量的演变高度非线性,而Jacobian 方法则是噪声大,空间中小向量的演变接近线性;小数据量方法的优点在于：1对小数据组的计算可靠；2计算量较小,比wolf方法快很多；3编程、操作较为容易;而关于时间延迟、嵌入维数、平均周期的确定,还是推荐C－C方法和G－P算法,结果更为可靠一些你这已经做得很好了;这个具体的应用我也是做得很少的,还有请注意一下matlab版中的一个帖子, 早就注意了,我这里面一些代码就是从那边套用的试过了吗都没有问题吧计算起来,原来论坛里面一些程序不能运行,套用,那你就要注意,在这里重写就要有我们的特色;不能够是简单的重请详细读帖子,上面有说明,当然不是简单的重复啦,我这里是理论结合实践,呵呵程序都测试能算的,不过对不同的系统,精度不敢保证的太高,呵呵,只要能运行就可以参照自己来编写精度可以根据自己需要来选择呵呵,确实,保证能够使用希望这样总结一下对大家有所帮助吧非常有价值;谢谢。
6.2 Lyapunov第二方法
V (x ) x 2
显 然 ， V ( x ) 0 x 0 ,且 V ( x ) 0 x 0 。
现 在 ， 我 们 考 虑 V 沿 上 述 微 分 方 程 的 解 对 时 间 t的 导 数 ， 有
V 2 x x 1 0 x2 0 x 0
由 于 V ( x ) 正 定 ， V 负 定 ， 从 而 渐 近 稳 定 。
定则说明：在任一点x处，v(x) 的值都是减小的，从
而在任一点x 处，运动的轨线都从v(x)=C的外部穿
越v(x)=C 走向内部。这表明，limt0x(t)=0，即原
点（零解）是渐近稳定精的选可。编辑ppt
8
练 习 ： P 277,2(1 ,2,3)
精选可编辑ppt
9
定理6 (不稳定性定理)
如果对微分方程组(6.20)存在函数V(x)和非负常数,满足
x2
x1
x2
x1
精选可编辑ppt
7
定理5(Lyapunov渐近稳定性定理) 如果对微分方程组(6.20)可以找到一个函数V(x),满足
（1）V(x)0
(2) V(x)0 则方程组(6.20)的零解是渐近稳定的.
几何解释：
由于v(x)正定， v(x)=C是一个闭的曲面族，层
层相套、随C 趋向于零而向原点退缩。而dv/dt 负
§6.2 李雅普诺夫第二方法
为了分析运动的稳定性,李雅普诺夫提出了两 种方法:
第一方法包含许多步骤，包括最终用微分方程 的显式解来对稳定性近行分析，是一个间接的方法。
第二方法不是求解微分方程组，而是通过构造 所谓李雅普诺夫函数（标量函数）来直接判断运动 的稳定性，因此又称为直接法。
李雅普诺夫第二方法目前仍是研究非线性、时 变系统最有效的方法，是许多系统控制律设计的 基本工具。
漫谈Lyapuno..
⎧ dx = f ( x, y ) ⎪ ⎪ dt ⎨ ⎪ dy = g ( x, y ) ⎪ ⎩ dt f (0, 0) = g (0, 0) = 0
(0, 0) , 不远离原点 (0, 0) , 远离原点 (0, 0) 来决定的,
而此信息等价于 V ( x(t ), y (t )) 是 t 的“下降”、“不增”、 “上升”函数, 而后者又等价于
1 华中科技大学控制科学与工程系, 武汉, 430070
1892 年 , 伟大的俄国数学力学专家亚历山大 · 米哈伊维 奇·Lyapunov 院士(A.M.Ляпунов), 完成了他的天才博士论文 《运 [1] 动稳定性的一般问题》 , 前四五十年之内, 它的影响范围还基 本上局限于前苏联的数学、 力学界, 甚至是曲高和寡, 束之高阁, 仅仅流连辗转在少数理论修养很高的数学力学专家手中, 后来 才得到自动控制界的青睐, 反响逐渐且越来越强烈. 1907 年论文 被译成法文, 而传至美国, 由 Princeton 大学 1947 再译成英文. 1954 年, 获美国哈佛大学数学博士学位, 后来回到国内数学研 究所工作的秦元勋教授将 Lyapunov 博士论文编译成中文讲义, 在数学所举办讲习班 , 系统地讲解研讨这篇博士论文 , 于是 Lyapunov 稳定性理论从前苏联传入中国, 开始在中国开花结果. 讲义以纪念 Lyapunov 诞生 100 周年的名义, 于 1958 年在(北京) 科学出版社出版[12]. 1952 年, 苏联著名的数学家马尔金(Малки)的专著[13] 《运动
漫谈 Lyapunov 稳定性的理论、方法和应用
lyapunov-krasovskii泛函方法
lyapunov-krasovskii泛函方法
Lyapunov-Krasovskii泛函方法是一种控制系统分析和设计中常用的数学方法。
该方法可以使用Lyapunov函数来评估系统的稳定性，并提供了计算误差界的工具。
此外，该方法还可以用于设计满足一定性能要求的控制器。
该方法的基本思想是通过构造Lyapunov-Krasovskii泛函来分析系统的稳定性。
该泛函是系统状态和控制输入的函数，并且具有一定的正定性和次级正定性条件。
使用Lyapunov-Krasovskii泛函，可以得到一些关于系统稳定性的定理。
这些定理
可以用于评估系统的稳定性，并采取控制策略来保证系统的稳定性。
在控制器设计方面，Lyapunov-Krasovskii泛函方法可以用于设计满足一定性能
要求的控制器。
这种控制器通常被称为LMIs（线性矩阵不等式）控制器，可以通
过使用LMIs工具箱来计算。
总之，Lyapunov-Krasovskii泛函方法是控制系统分析和设计中非常重要的数学
方法。
它可以用于评估系统的稳定性，提供了计算误差界的工具，以及设计满足
一定性能要求的控制器。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
推广 形式：
u w* w* u du,
w
w w*
u w* u w du, * du, w f u
*
w
w
w*
f (u ) f ( w* ) du, f u
w
w*
f (u ) f ( w* ) du. u
Lyapunov 函数的构造
四、一类SI模型的全局稳定性
x*
u
y*
u
ky*
v*
u
c
w*
?
u
dL 则 1 dt
* p 0 x x * * x 2 w w 0 * c ww* x x M
2
x* xy*v y v* x v 3 x x* yv* y *v 0,
Lyapunov 函数的一些构造方法
李建全
空军工程大学(西安)
E-mail: jianq_li@
Lyapunov 函数的构造
一、引言 二、带有免疫的病毒动力学模型 三、Lyapunov函数形式的推广 四、一类SI模型的全局稳定性
Lyapunov 函数的构造
一、引言
全局渐近稳定性 (GAS)是动力系统研究的重要 问题。 对于高维自治微分系统 ，证明GAS最常用的方法是 Lyapunov - LaSalle定理， 即若自治系统 * x' f x, x Rn
Lyapunov 函数的构造
二、带有免疫的病毒动力学模型
x 健康细胞； y 被感染细胞； v 游离病毒； w 免疫细胞.
其中f w可分别取w,1, w /1 w.
M有惟一正平衡点 假设系统 E* x* , y* , v* , w* ,
下面通过构造 Lyapunov 函数来证明其GAS.
* *
dL 并且 1 dt
y v 0 x x* , * * , w w* . y v M
OK
Lyapunov 函数的构造
函数 3、 情形f w w / 1 w. 定义Lyapunov * * x ux y u y p w 1 u u w* x*v* v u v* L du du du du
结论：
Lyapunov 函数的构造
由于
则(IN)变为
Lyapunov 函数的构造
(IN) 定义Lyapunov函数 则有
OK
Lyapunov 函数的构造
补充：
取
条件：
Lyapunov 函数的构造
上面介绍的构造Lyapunov函数方法也可 用于如下模型：
和
A class of Lyapunov functions and the global stability of some epidemic models with nonlinear incidence, Journal of Applied Analysis and Computation , 6(2016)38-46
Lyapunov 函数的构造
谢 谢 ！
函数 2、 情形f w 1. 定义Lyapunov
*
dL 则 1 dt
* x x p 0 * * x 2 w w * x x M c
2ຫໍສະໝຸດ x* xy*v y v* x v 3 x x* yv* y *v 0,
x * y * * * v *
* * *
函数 1、 情形f w w. 定义Lyapunov
Nonlinear Analysis—Real World Applications: 12(2011)2163-2173; 13(2012)2006-2016.
Lyapunov 函数的构造
* y u y p w u w* 2 p u x* x*v* v u v* w w du ,, L* du * du du * * * x y v u u ky u 2 ccw u x
Lyapunov 函数的构造
* w p u w ux u y x v u v du, L du du du w* * x y v c u u u ky u dL1 x x* p0 w w* * 2 * 则 x 2 * c w w dt M x x * * * x xy v y v * * x v 3 x x* yv* y *v 0, v dL1 * y * 并且 0 x x , * * , w w . OK y v dt M
在区域内有惟一平衡点 E *， 并且能找到Lyapunov 函数L， 使得dL/dt * 在内关于E *负定或半负定， 同时系统(*)在集 合 dL/dt * 0 上的最大不变集为E * ， 则 E * 在内GAS。
证明的关键在于构造恰当的Lyapunov函数， 并证明其沿系统解的全导数是负定或半负定的。
* *
dL 并且 1 dt
y v 0 x x* , * * , w w* . y v M
OK
Lyapunov 函数的构造
三、Lyapunov函数形式的推广
* y u y u x* x*v* v u v* 其中 F * du * du du. * * x y v u u ky u x