longinv.m

来自「小波变换所有常见例程」· M 代码 · 共 38 行

function Pinv = longinv(P)

% LONGINV -- brute force inverse of matrix polynomial
%
%        Pinv = longinv(P)
%
% Caution: this is done by the adjoint formula for matrix inverses
% (Cramer's rule), so it is extremely slow for large matrices.
% Use only as a last resort, for cases where INV does not work.

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

D = det(P);
if (size(D.coef,3) > 1)
    error('det(P) is not a monomial');
end
if (iszero(D))
    error('det(P) = 0');
end

n = size(P,1);
Pinv = P;
for i = 1:n
    for j = 1:n
%       Pinv(i,j) = (-1)^(i+j) * det(P([1:j-1,j+1:n],[1:i-1,i+1:n]));
	Pinv = subsasgn(Pinv,substruct('()',{i,j}),...
	    (-1)^(i+j) * det(subsref(P,substruct('()',{[1:j-1,j+1:n],[1:i-1,i+1:n]}))));
    end
end

Pinv = trim(Pinv / D);