一、连续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');