西尔维斯特方程的形式：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\'
请发表评论