lyap.m

数字通信第四版原书的例程

function X = lyap(A, B, C)
%LYAP	Lyapunov equation.
%	X = LYAP(A,C) solves the special form of the Lyapunov matrix 
%	equation:
%
%		A*X + X*A' = -C
%
%	X = LYAP(A,B,C) solves the general form of the Lyapunov matrix
%	equation:
%
%		A*X + X*B = -C
%
%	See also DLYAP.

%	S.N. Bangert 1-10-86
%	Copyright (c) 1986-93 by The MathWorks, Inc.
%	Last revised JNL 3-24-88

	if nargin == 2
		C = B;
		B = A';
	end
	[ma,na] = size(A);
	[mb,nb] = size(B);
	[mc,nc] = size(C);

	if (ma ~= na) | (mb ~= nb) | (mc ~= ma) | (nc ~= mb)
		error('Dimensions do not agree.')
	elseif ma==0,
		X = []; return
	end

% check if problem has any complex inputs
	real_flg = 1;
	if any(any(imag(A))) | any(any(imag(B))) | any(any(imag(C)))
		real_flg = 0;
	end

% perform schur decomposition on A and B (note: complex schur form forced by
% adding small complex part so ua and ub are complex upper triangular)
	[ua,ta] = schur(A);
	[ua,ta] = rsf2csf(ua,ta);
	[ub,tb] = schur(B);
	[ub,tb] = rsf2csf(ub,tb);

% check all combinations of ua(i,i)+ub(j,j) for zero
	p1 = ones(mb,1)*diag(ta)';
	p2 = diag(tb)*ones(1,ma);
	sum = abs(p1) + abs(p2);
	if any(any(sum == 0)) | any(any(abs(p1 + p2) < 1000*eps*sum))
		error('Solution is not unique.');
	end

% transform C
	ucu = -ua'*C*ub;

% solve for first column of transformed solution
	y(ma,mb) = 0;
	ema = eye(ma);
	y(:,1) = (ta+ema*tb(1,1))\ucu(:,1);

% solve for remaining columns of transformed solution
	for k=2:mb
		km1 = 1:(k-1);
		y(:,k) = (ta+ema*tb(k,k))\(ucu(:,k)-y(:,km1)*tb(km1,k));
	end

% find untransformed solution 
	X = ua*y*ub';

% ignore complex part if real inputs (better be small)
	if real_flg
		X = real(X);
	end