matlab練習程序（解西爾維斯特、李雅普諾夫方程）

西爾維斯特方程的形式：AX+XB=C

李雅普諾夫方程的形式：AX+XA'=-C

這兩種方程都是已知矩陣A,B,C,求解X的方程。

對於這種方程有兩種方法來求解，一種是朴素法，一種是Bartels-Stewart法。

以西爾維斯特方程為例，朴素法是將方程寫為下列形式進行直接求解：

其中圓圈中帶個X的符號是kron積，vec()是將X或C轉換為了n*1的列向量。

該方法將原來O(n^3)的問題變為了O(n^6)，如果矩陣比較大，估計速度會比較慢。

第二種方法為Bartels-Stewart法，下面以西爾維斯特方程為例介紹一下：

首先我們對A和B‘進行shur分解：

原方程可改寫為：

此時令：

 

得到：

此時R和S都是一個上三角矩陣，我們需要S作為下三角矩陣才能方便求解。

這里的S正好是我們是對B'的shur分解，由於shur分解的特性，shur(B)=VSV'，shur(B')=VS'V'，所以這里再對S求個轉置即可。

得到類似下面的矩陣方程：

展開得到：

再依次求出y4,y3,y2,y1即可。

matlab代碼如下：

解西爾維斯特方程：

clear all;close all;clc;

A = rand(2,2);
X = rand(2,2);
B = rand(2,2);
C = A*X+X*B;

X

%%系統函數
X = sylvester(A,B,C)

%%朴素法，自寫kron積
IA = [A zeros(2);zeros(2) A];
BI = [B(1,1)*eye(2) B(2,1)*eye(2);
      B(1,2)*eye(2) B(2,2)*eye(2)];
X = reshape(inv(IA+BI)*C(:),[2,2])

%%朴素法，系統kron積
X = reshape(inv(kron(eye(2),A)+kron(B',eye(2)))*C(:),[2,2])

%%Bartels–Stewart法
[U,R] = schur(A);   %schur分解，R是上三角
[V,S] = schur(B');
S = S';             %S是下三角
F = U'*C*V;

%解R*Y + Y*S = F方程
Y = zeros(2,2);
Y(2,2) = F(2,2)/(R(2,2) + S(2,2));
Y(2,1) = (F(2,1) - S(2,1)*Y(2,2)) / (S(1,1) + R(2,2));
Y(1,2) = (F(1,2) - R(1,2)*Y(2,2)) / (R(1,1) + S(2,2));
Y(1,1) = (F(1,1) - R(1,2)*Y(2,1) - S(2,1)*Y(1,2)) / (R(1,1) + S(1,1));

X = U*Y*V'

解李雅普諾夫方程：

clear all;close all;clc;

A = rand(2);
X = rand(2);
C = -(A*X+X*A');

X

%%系統函數
X = lyap(A,C)

%%朴素法，自寫kron積
IA = [A zeros(2);zeros(2) A];
BI = [A(1,1)*eye(2) A(1,2)*eye(2);
      A(2,1)*eye(2) A(2,2)*eye(2)];
X = reshape(inv(IA+BI)*(-C(:)),[2,2])

%%朴素法，系統kron積
X = reshape(inv(kron(eye(2),A)+kron(A,eye(2)))*(-C(:)),[2,2])

%%Bartels–Stewart法
[U,R] = schur(A);   %schur分解，R是上三角
S = R';             %S是下三角
F = U'*(-C)*U;

%解R*Y + Y*S = F方程
Y = zeros(2,2);
Y(2,2) = F(2,2)/(R(2,2) + S(2,2));
Y(2,1) = (F(2,1) - S(2,1)*Y(2,2)) / (S(1,1) + R(2,2));
Y(1,2) = (F(1,2) - R(1,2)*Y(2,2)) / (R(1,1) + S(2,2));
Y(1,1) = (F(1,1) - R(1,2)*Y(2,1) - S(2,1)*Y(1,2)) / (R(1,1) + S(1,1));

X = U*Y*U'