一、连续Lyapunov方程

连续Lyapunov方程可以表示为



Lyapunov方程来源与微分方程稳定性理论，其中要求C为对称正定的n×n方阵，从而可以证明解X亦为n×n对称矩阵，这类方程直接求解比较困难，不过有了Matlab中lyap()函数，就简单多了。

>> A=[1 2 3;4 5 6;7 8 0]

A =

    1     2     3
    4     5     6
    7     8     0

>> C=-[10 5 4;5 6 7;4 7 9]

C =

  -10    -5    -4
   -5    -6    -7
   -4    -7    -9

>> X=lyap(A,C)

X =

  -3.9444    3.8889    0.3889
   3.8889   -2.7778    0.2222
0.3889    0.2222   -0.1111

二、Lyapunov方程的解析解

利用Kroncecker乘积的表示方法，可以将Lyapunov方程写为

function x=lyap2(A,C)
%Lyapunov方程的符号解法
n=size(C,1);
A0=kron(A,eye(n))+kron(eye(n),A);
c=C(:);
x0=-inv(A0)*c;
x=reshape(x0,n,n)

例子

>>A=[1 2 3;4 5 6;7 8 0];
>>C=-[10 5 4;5 6 7;4 7 9];
>>x=lyap2(sym(A),sym(C))

x =

[ -71/18,   35/9,   7/18]
[   35/9,  -25/9,    2/9]
[   7/18,    2/9,   -1/9]

三、离散Lyapunov方程

离散Lyapunov方程的一般形式为



Matlab中直接提供了dlyap()函数求解该方程：X=dlyap(A,Q)

其实，如果A矩阵非奇异，则等式两边同时右乘得到



就可以将其变换成连续的Sylvester方程



而Sylvester方程是广义Lyapunov方程，故离散的Lyapunov方程还可以使用下面的方法求解

B=-inv(A’)
C=Q*inv(A’)
X=lyap(A,B,C)

下面总结下我们上面的讲到的知识点：

X=lyap(A,C)                               连续Lyapunov方程数值解法

X=lyap2(A,C)                             连续Lyapunov方程符号解法

X=lyap(A,B,C)                              广义Lyapunov方程，即Sylvester方程

X=dlyap(A,Q)或者X=lyap(A,-inv(A’),Q*inv(A’))    离散Lyapunov方程



Sylvester方程Matlab求解Sylvester方程的一般形式为



该方程又称为广义的Lyapunov方程，式中A是n×n方阵，B是m×m方阵，X和C是n×m矩阵。Matlab控制工具箱提供了直接的求解该方程的lyap()函数




A=[8 1 6;3 5 7;4 9 2]
B=[2 3;4 5]
C=[1 2;3 4;5 6]
X=lyap(A,B,C)


A =

    8     1     6
    3     5     7
    4     9     2


B =

    2     3
    4     5


C =

    1     2
    3     4
    5     6


X =

   0.2011    0.2016
   0.0393    0.1554
  -0.6428   -0.8966


同理，我们使用Kronecker乘机的形式将原方程进行如下变化



故可以编写Sylvester方程的解析解函数


function X=lyap3(A,B,C)
%Sylvester方程的解析解法
%rewrited by dynamic
%more information http://www.ilovematlab.cn

If nargin==2,C=B;B=A';end
[nr,nc]=size(C);
A0=kron(A,eye(nc))+kron(eye(nr),B');
try 
   C1=C';
   X0=-inv(A0)*C1(:);
   X=reshape(X0,nc,nr);
catch
   error('Matlabsky提醒您：矩阵奇异！');
end



用上面的数据，我们试验下该解析解法的


>>X=lyap3(sym(A),B,C)

X =

[   2853/14186,    557/14186,  -9119/14186]
[  11441/56744,   8817/56744, -50879/56744]


Riccati方程的Matlab求解Riccati方程是一类很著名的二次型矩阵形式，其一般形式为



由于含有矩阵X的二次项，所有Riccati方程求解要Lyapunov方程更难，Matlab控制工具箱提供了are()函数，可以直接求解该函数


A=[-2 1 -3;-1 0 -2;0 -1 -2]
B=[2 2 -2;-1 5 -2;-1 1 2]
C=[5 -4 4;1 0 4;1 -1 5]
X=are(A,B,C)


A =

   -2     1    -3
   -1     0    -2
    0    -1    -2


B =

    2     2    -2
   -1     5    -2
   -1     1     2


C =

    5    -4     4
    1     0     4
    1    -1     5


X =

   0.9874   -0.7983    0.4189
   0.5774   -0.1308    0.5775
  -0.2840   -0.0730    0.6924


如何用matlab求解lyapunov指数我是需要分析计算LOGISTIC数据  ,都是用来说明对初值的敏感.以下是LOGISTIC求解的程序 ,希望得到LYAPUNOV的程序.
clc 
clear 
close all 

lambda = 3:5e-4:4; 
x = 0.4*ones(1,length(lambda)); 

N1 = 400;                  % 前面的迭代点数
N2 = 100;                  % 后面的迭代点数

f = zeros(N1+N2,length(lambda)); 
for i = 1:N1+N2 
   x = lambda .* x .* (1 - x); 
   f(i,:) = x; 
end 
f = f(N1+1:end,:); 
plot(lambda,f,'k.','MarkerSize',1) 
xlabel('lambda') 
ylabel('x');

转载自 http://blog.sciencenet.cn/blog-1373583-798129.html 感谢作者！