CondNumInvPairs := proc(coeffs,X,S,C,R) # Input: coeffs - cell array containing coefficients of matrix polynomial: # P(x) = A_l x^l + ...+ A_2 x^2 + A_1 x + A_0... coeffs := [A_0,...,A_l] # X - nxk matrix invariant pair: P(X,S)=0 # S - kxk matrix invariant pair: P(X,S)=0 # C - center of the circle # R - radius of the circle # Output: cond - condition number of P(X,S) local n,len,P,i,invP,ga,alp,Fp,Bk,N,M,M1,B,cond; n := Size(coeffs[1],2): len := numelems(coeffs)-1: P:=0: Bk:=[]; ########### Compute the matrix polynomial P(z) ############### for i from 1 by 1 to len+1 do P:= P+z^(i-1)*coeffs[i]; end do: ############################################################## ####### Compute gamma ####### ga := C+R*exp(Complex(1)*t); ############################# ####### Compute inverse of (ga*I-S) ######### invP:= MatrixInverse(ga*IdentityMatrix(n)-S); ############################################# ###################### Compute alpha_i’s ######################## alp:= [seq(Norm(coeffs[i], Frobenius), i = 1..len+1)]; ################################################################# ############################## Compute matrices B_k ############################### for i from 0 to len do Fp := KroneckerProduct(Transpose(ga^(i)*X.invP),IdentityMatrix(n))*exp(Complex(1)*t); Bk := [Bk[], (R/(2*Pi))*map(int, evalc(Fp), t=0..(2*Pi))]; end do: ################################################################################### ######################### Compute matrices M and N ############################ Fp:= KroneckerProduct(Transpose(invP), subs(z = ga,P).X.invP)*exp(Complex(1)*t); M:= (R/(2*Pi))*map(int, evalc(Fp), t=0..(2*Pi)); Fp:= KroneckerProduct(Transpose(invP), subs(z = ga,P))*exp(Complex(1)*t); N:= (R/(2*Pi))*map(int, evalc(Fp), t=0..(2*Pi)); ############################################################################### ################################ Compute matrix B ############################# B:= convert([seq(alp[len+2-i]*Bk[len+2-i], i = 1..len+1)], Matrix); M1:= convert([N,M], Matrix); ############################################################################### ###################################### Compute condition number ################################################ cond:= evalf( Norm( MatrixInverse(M1).B ,2)/Norm( convert([Norm(X, Frobenius), Norm(S, Frobenius)],Vector) ,2) ); ################################################################################################################ return cond; end proc: