www.gusucode.com > matlab矩阵函数工具箱 > matlab矩阵函数工具箱/matlab矩阵函数工具箱/rootpm_real.m
function X = rootpm_real(A,p) %ROOTPM_REAL Pth root of real matrix via real Schur form. % X = ROOTPM_REAL(A,P) is a Pth root of the nonsingular matrix A % (A = X^P). When A has no eigenvalues on the negative real axis % then X is the principal pth root, which is real when A is real. % Implementation is complicated by use of zero (and even "-1") % subscripts in algorithm. Dealt with by increasing "k" by 1 % each time and moving "-1" cases in front of loops. if norm(imag(A),1), error('A must be real.'), end if p == 1, X = A; return; end n = length(A); [Q,R] = schur(A,'real'); % m blocks: i'th has order s(i), starts at t(i). [m,s,t] = quasitriang_struct(R); U = zeros(n); B = cell(p); V = cell(p-1); V{1} = eye(n); % U^k == V{k+1}. for j=1:m rj = t(j):t(j)+s(j)-1; if s(j) == 1 % If A is real, X will be real unless R(p,p)<0 on the next line. U(rj,rj) = R(rj,rj)^(1/p); else U(rj,rj) = root_block(R(rj,rj),p); end for k = 0:p-2, kk = k+1; V{kk}(rj,rj) = U(rj,rj)^(k+1); end for i=j-1:-1:1 ri = t(i):t(i)+s(i)-1; for k = 0:p-2 kk = k+1; B{kk} = zeros(s(i),s(j)); for ell = i+1:j-1 rell = t(ell):t(ell)+s(ell)-1; B{kk} = B{kk} + U(ri,rell)*V{kk}(rell,rj); end end rhs = R(ri,rj) - B{p-2 +1}; for k = 0:p-3 rhs = rhs - V{p-3-k +1}(ri,ri)*B{k+1}; end coeff = kron( eye(s(j)), V{p-2 +1}(ri,ri) ) ... + kron( V{p-2 +1}(rj,rj).', eye(s(i)) ); for k = 1:p-2 coeff = coeff + kron( V{k-1 +1}(rj,rj).', V{p-2-k +1}(ri,ri) ); end y = coeff\rhs(:); rhs(:) = y; % `Un-vec' the solution. U(ri,rj) = rhs; for k = 0:p-2 if k == 0 S = U(ri,rj); else S = V{k-1 +1}(ri,ri)*U(ri,rj) + U(ri,rj)*V{k-1 +1}(rj,rj) ... + B{k-1 +1}; for ell = 1:k-1 S = S + V{k-ell-1 +1}(ri,ri)*U(ri,rj)*V{ell-1 +1}(rj,rj); end for ell = 0:k-2 S = S + V{k-2-ell +1}(ri,ri)*B{ell +1}; end end V{k+1}(ri,rj) = S; end end end X = Q*U*Q'; %%%%%%%%%%%%%%%%%%%%%%%%%%%% function X = root_block(R,p) %ROOT_BLOCK Pth root of a real 2x2 matrix with complex conjugate eigenvalues. if norm(imag(R)) ~= 0 || any(size(R) - [2,2]) error('Matrix must be real, of dimension 2.') end r11 = R(1,1); r12 = R(1,2); r21 = R(2,1); r22 = R(2,2); theta = (r11 + r22) / 2; musq = (-(r11 - r22)^2 - 4*r21*r12) / 4; mu = sqrt(musq); if musq <= 0 error('Matrix must have non-real complex conjugate eigenvalues.') end r = sqrt(theta^2+musq); phi = angle(complex(theta,mu)); rootp = r^(1/p)*exp(i*phi/p); alpha = real(rootp); beta = imag(rootp); X = alpha*eye(2) + (beta/mu)*(R - theta*eye(2));