To execute an example as defined, one needs to define a proper Maple environmentrestart;with(LinearAlgebra): with(MatrixPolynomialAlgebra):Next, one has to activate all procedures defined in the next section.

For the sake of convenience, we have included also all procedures defined in section 11.6 which are needed in this section, especially to execute the example below.LCMDenomMatrixPolynom:=proc(W,x) local k,m,nc,D0,D1,cofmax,mat; mat := convert(W(x),Matrix): nc := ColumnDimension(mat): D0:=[]: for k from 1 to nc do for m from 1 to nc do D0:=[op(D0),denom(mat[k,m])]; end do end do; D1:=lcm(op(D0)): cofmax:=coeff(D1,x,degree(D1,x)): return(expand(D1/cofmax));end;GetPolesandZeros:=proc(pf,zf) local poles, zeros, mu; poles:=[solve(pf(x),x)]: print('poles'=poles); zeros:=[solve(zf(x),x)]: print('zeros'=zeros); mu:=convert([op({op(convert(poles,list)), op(convert(zeros,list))})],Vector[row]): print(`Set of different poles and zeros`=mu); return([convert(poles,Vector[row]), convert(zeros,Vector[row]),mu]); end proc;GetMultiplicity:=proc(p::Vector,mu::Vector) local dimp,dimgv,mP,k,m; dimp:=Dimension(p): dimgv:=Dimension(mu): mP:=Vector[row](dimgv,0): for k from 1 to dimgv do for m from 1 to dimp do if evalb(p[m]=mu[k]) then mP[k]:=mP[k]+1; end if: end do:end do: return(mP);end;GetAllMOrderings:=proc(mA::Vector,mZ::Vector,mu::Vector) local h,kordering,orderings, newperm; newperm:=true:h:=0:kordering:=0:while (kordering=0) do h:=h+1:orderings:=GetMOrderingsH(mA,mZ,mu,h,newperm): newperm:=false: kordering:=orderings[1]: end do: return(h,orderings);end;GetMOrderingsH:=proc(muA,muZ,mu,h,newperm) local nmu,nperm,kk,k,perm,orderingP,orderingZ; global allperm; nmu:= Dimension(mu): if newperm=true then allperm:= combinat[permute](nmu): end if: nperm:= combinat[numbperm](nmu): kk:=0: for k from 1 to nperm do perm:=allperm[k]: if (TestOrderingMAZ(muA,muZ,perm,h)) then kk:=kk+1: orderingP[kk]:=GetOrderedVector(muA,mu,convert(perm,list)): orderingZ[kk]:=GetOrderedVector(muZ,mu,convert(perm,list)): end if: end do: if kk>0 then return(kk,orderingP,orderingZ) else return(kk); end if: end proc;TestOrderingMAZ:= proc(ma,mz,perms,h) local s, bol, t; s:=Dimension(mz): bol:=true: t:=0: while (t<s) and bol do t:=t+1: bol:=not((add(mz[perms[m1]],m1=1..t)- add(ma[perms[m2]],m2=1..(t-1)))>(h)): end do;return(bol);end;GetOrderedVector:=proc(ma,mu,ordering) local no,n1,tk,nc, k, m, mup, mpp, orderedV; mup:=mu[ordering]: mpp:=ma[ordering]: no:=Dimension(mu): n1:=add(ma[k],k=1..no): nc:=0: orderedV:=Vector[row](n1,0): for k from 1 to no do tk:=mpp[k]: for m from 1 to tk do nc:=nc+1: orderedV[nc]:=mup[k]: end do: end do: return(orderedV);end proc;MakeFactorization:=proc(A,B,C,x) local WW, nv, m, Im, k, R, oa; m:=RowDimension(C);nv:=Dimension(A):Im:=IdentityMatrix(m): WW:=Vector(nv): for k from 1 to nv[1] do R:=MakeRmatrix(B,C,k): WW[k]:=unapply(map(factor,Im+ScalarMultiply(R,1/(x-A[k,k]))),x): end do: return(WW);end;MakeRmatrix:=proc(B,C,k) local nc, mat, i, j; nc:=RowDimension(C); mat:=Matrix(nc,nc,0); for i from 1 to nc do for j from 1 to nc do mat[i,j]:=C[i,k]*B[k,j]; end do end do: return(mat);end;Factors2Transfer := proc(AllFactors,x) local Wtest,DimS,k,n,Wdum,ResultW; DimS:=ColumnDimension(AllFactors[1](x)): n:=Dimension(AllFactors):Wtest:=IdentityMatrix(DimS): for k from 1 to n do Wtest:=map(simplify,Wtest.AllFactors[k](x)): end do: Wdum:=convert(map(factor,map(simplify,evalm(Wtest))),Matrix): ResultW:=unapply(Wdum,x): return(ResultW);end proc;End of procedures which are defined already in section 11.6Johnson's ruleord:= proc(x) local i,s; s:=sort([seq([x[i],i],i=1..nops(x))],(a,b)->evalb(a[1]<b[1])); [map(x->x[1],s),map(x->x[2],s)] end;JohnsonRule:=proc(mA,mZ,mu) local V1,V2,VZ1,VA2,S,lv1,nm,k,m1,m2,h; nm:=Dimension(mA): VZ1:=[]: VA2:=[]: V1:=[]: V2:=[]: for k from 1 to nm do if (mZ[k]<mA[k]) then V1:=[op(V1),k]: VZ1:=[op(VZ1),mZ[k]] : else V2:=[op(V2),k]: VA2:=[op(VA2),-1*mA[k]]: end if: end do: lv1:=ord(VZ1)[2]: S:=[]: for k from 1 to nops(V1) do S:=[op(S),V1[lv1[k]]]: end do: lv1:=ord(VA2)[2]: for k from 1 to nops(V2) do S:=[op(S),V2[lv1[k]]]: end do: h:=1: while (h<nm) and not (TestOrderingMAZ(mA,mZ,S,h)) do h:=h+1: end do: return(h,1,[GetOrderedVector(mA,mu,S),GetOrderedVector(mZ,mu,S)]); end proc;Procedures specific for section 12.5QCmatrices:=proc(Am,Bm,Cm,Dm,dqw,gammav) local v, k, Kmat, Xmat, Fmat, Gmat, Ahat, Bhat, Chat, Dhat, LastRowA, na; if not (IsCompanionForm(Am)) then error "First input matrix should be in companion form but get %1",Am; end if; na:=Dimension(Am): if (dqw>na[2]) then LastRowA:=Row(Am,na[1]): v:=Vector([0,1]); Kmat:=MakeKmatrix(LastRowA,v,gammav): Xmat:=MakeXmatrix(gammav,dqw,na[2]): Gmat:=DiagonalMatrix(gammav[1..(dqw-na[2])],(dqw-na[2]), (dqw-na[2])): Fmat:=Kmat.Xmat: Ahat:=<<Am,ZeroMatrix(dqw-na[2],na[2])>|<Bm.Fmat,Gmat>>: Bhat:=<Bm,ZeroMatrix(dqw-na[2],ColumnDimension(Bm))>: Chat:=<<Cm>|<Fmat>>: else Ahat:=Am: Bhat:=Bm: Chat:=Cm: end if: Dhat:=Dm: return(Ahat,Bhat,Chat,Dhat); end proc;MakeGamma:=proc(dw,mu,symbols) local nmu,mm,k,numval,gammav; if symbols then gammav:=Vector(dw,symbol=gamma): else gammav:=Vector(dw): nmu:=Dimension(mu):k:=0: numval:=[0]: for k from 1 to nmu do if (type(mu[k],complex)) then numval:=[op(numval),Re(mu[k])]:end if:end do: mm:=max(op(numval)): for k from 1 to dw do gammav[k]:=mm+k:end do:end if: return(gammav);end proc;IsCompanionForm:=proc(A) local na, zm, tm, k, bol; na:=Dimension(A): zm:=Matrix(na[1]-1,na[2],0): for k from 1 to (na[1]-1) do zm[k,k+1]:=1: end do: tm:=SubMatrix(A,[1..na[1]-1],[1..na[2]]): bol:=Equal(map(simplify,(tm-zm)),ZeroMatrix(na[1]-1,na[2])): return(bol);end proc;GenVandermondeMatrix:=proc(r,av::Vector) local GVM,m,v,j,c; c:=Dimension(av):v:=0:GVM:=GenVandermondeVector(r,av[1],v): for j from 2 to c do v:=0:for m from 1 to (j-1) do if (av[m]=av[j]) then v:=v+1: end if: end do: GVM:=<GVM|GenVandermondeVector(r,av[j],v)>:end do: return(GVM);end proc;GenVandermondeVector:=proc(r,a,k) local colgvm,j; colgvm:=Vector(r,0): for j from 1 to r do if (j>k) then colgvm[j]:=binomial(j-1,k)*a^(j-1-k):end if: end do: return(colgvm); end proc;MakeKmatrix:= proc(a,v,gammav) local nc,k,Kmat,pol; nc:=Dimension(a):pol:=1: for k from 1 to (nc) do pol:=pol*(lambda-gammav[k]):end do: pol:=expand(pol): Kmat:=ScalarMultiply(v,-a[1]-coeff(expand(pol),lambda,0)): for k from 2 to (nc) do Kmat:=<Kmat|ScalarMultiply(v,-a[k]-coeff(expand(pol),lambda,k-1))>: end do:return(Kmat);end proc;MakeXmatrix:=proc(gammav,d,n) local gd,k,m,X; X:=GenVandermondeVector(n,gammav[1],0): for k from 2 to (d-n) do X:=<X|GenVandermondeVector(n,gammav[k],0)>: end do: return(convert(X,Matrix)); end proc;MakeBasisA:=proc(orderedA,dqw,gammav) local BasisMat,id,zerod,n,j,jj,GVM,npoles,k; npoles:=Dimension(orderedA): GVM:=GenVandermondeMatrix(npoles,orderedA): n:=Dimension(GVM)[1]:if dqw>n then id:=IdentityMatrix(dqw-n):zerod:=Vector(dqw-n,0): BasisMat:=<GenVandermondeVector(n,gammav[1],0),Column(id,1)>: for j from 2 to (dqw-n) do BasisMat:= <BasisMat|<GenVandermondeVector(n,gammav[j],0),Column(id,j)>>: end do: for j from (dqw-n+1) to dqw do jj:=j+n-dqw: BasisMat:=<BasisMat|<Column(GVM,jj),zerod>>:end do: else BasisMat:=GVM:end if: return(BasisMat);end proc;MakeBasisZ:=proc(orderedZ,dqw,gammav) local BasisMat,id,zerod,n,j,jj,nzeros,ReorderedZ,GVM,k; nzeros:=Dimension(orderedZ): ReorderedZ:=ReverseOrder(orderedZ): GVM:=GenVandermondeMatrix(nzeros,ReorderedZ): n:=Dimension(GVM)[1]:if dqw>n then id:=IdentityMatrix(dqw-n):zerod:=Vector(n,0): BasisMat:=<zerod,Column(id,1)>: for j from 2 to (dqw-n) do BasisMat:=<BasisMat|<zerod,Column(id,j)>>: end do: zerod:=Vector(dqw-n,0):for j from (dqw-n+1) to dqw do jj:=dqw+1-j:jj:=j-dqw+n: BasisMat:=<BasisMat|<Column(GVM,jj),zerod>>: end do:else BasisMat:=GVM:end if; return(BasisMat);end proc;ReverseOrder:=proc(v) local rv,n,k; n:=Dimension(v): rv:=Vector(n): for k from 1 to n do rv[k]:=v[n-k+1]: end do: return(rv);end proc;UpperLowerTransformation:=proc(TA,TZ,dqw) local S1,S,k; S:=op(IntersectionBasis([[Column(TA,1)],[Column(TZ,[1..dqw])]])): for k from 2 to dqw do S1:=op(IntersectionBasis([[Column(TA,[1..k])],[Column(TZ,[1..dqw-k+1])]])): S:=<S|S1>:end do: return(map(simplify,S));end proc;

Construction of WExample of a transfer function W, defined in symbolic variablesn:=3; #size of Wp:='p':r:='r':q:='q': sw1:=<<1,0>|<'s'(lambda)/ 'q(lambda)', '''q^(x)''(lambda)'/ 'q(lambda)'>>: 'W(lambda)'=sw1;# ppol is a routine to create polynomialsppol:=proc(x,p,nn) local r,k; r:=1: for k from 1 to (nn) do r:= r*(x-p[k-1]): end do: return(r); end proc;p:='p': p:=array(0..(n-1)): pp:=x->ppol(x,p,n): 'q(lambda)'=pp(lambda); z:='z': z:=array(0..(n-1)): zz:=x->ppol(x,z,n): 'q^(x)'(lambda) = zz(lambda); r:='r': r:=array(0..(n-1)): rs:=x->ppol(x,r,n-2): 's(lambda)'=rs(lambda);p:='p': z:='z': r:='r': p[1]:=p[0]; z[0]:=p[0]; z[1]:=p[1]; z[2]:=p[2];W:=x-><<1,0>|<rs(x)/pp(x),zz(x)/pp(x)>>:'W(lambda)'=W(lambda);End of construction of WStart calculation of factorization for example defined above:DimS:=ColumnDimension(W(x));q:=unapply(LCMDenomMatrixPolynom(W,x),x):ppoles:=q:pzeros:=unapply(simplify(W(x)[2,2]*q(x)),x):'p(lambda)'=sort(collect(ppoles(lambda),lambda),lambda);'(p^(x))(lambda)'=sort(collect(pzeros(lambda),lambda),lambda);res1:=GetPolesandZeros(ppoles,pzeros):poles := res1[1]: zeros:=res1[2]: mu:=res1[3]:npoles:=Dimension(poles);nzeros:=Dimension(zeros); nmu:=Dimension(mu);Construct companion based realization:r:=unapply(simplify(W(x)[1,2]*ppoles(x)),x):r(lambda):Amin:= Transpose(CompanionMatrix(ppoles(x),x)):Bmin:= Matrix(npoles,DimS,0):Bmin[npoles,DimS]:=1:Cmin:= Matrix(DimS,npoles): for k from 1 to npoles do Cmin[1,k]:= coeff(r(lambda),lambda,k-1): Cmin[2,k] := coeff(pzeros(lambda)-ppoles(lambda),lambda,k-1): end do:Dmin:=IdentityMatrix(DimS,DimS):Amincross:=Transpose(CompanionMatrix(pzeros(x),x)):'A'=convert(Amin,array); `transpose(B)`=convert(Transpose(Bmin),array); `C`=convert(Cmin,array);Getting orderingsmuA:=GetMultiplicity(poles,mu); muZ:=GetMultiplicity(zeros,mu);Example using Johnson's Rule (see Johnson's Rule):L1:=JohnsonRule(muA,muZ,mu);L1[3][1];Continue with main example (not using Johnson's rule):ResultOrdering:=GetAllMOrderings(muA,muZ,mu);In this example, we take the first ordering. h:=ResultOrdering[1];print(`number of orderings`=ResultOrdering[2]); orderedA:=ResultOrdering[3][1]; orderedZ:=ResultOrdering[4][1];dw := npoles: dqw:=h-1+npoles:'delta[q](W)' = dqw;print('delta[q](W)-delta(W)'=dqw-dw);Triangularization quasicomplete:gammav:=MakeGamma(dw,mu,true); #symbolic!Allhat:=QCmatrices(Amin,Bmin,Cmin,Dmin,dqw,gammav):Ahat := map(simplify,Allhat[1]): Bhat := Allhat[2]:Chat := map(simplify,Allhat[3]): Dhat := Allhat[4]:Ahatcross:=Ahat-Bhat.Chat:'Ahat'=convert(Ahat,array);TA:=MakeBasisA(orderedA,dqw,gammav):TZ:=MakeBasisZ(orderedZ,dqw,gammav):S:=UpperLowerTransformation(TA,TZ,dqw):'S'=convert(S,array);Sinv:=MatrixInverse(S):Atr:=convert(map(simplify,Sinv.Ahat.S),Matrix):Atrcross:=convert(map(simplify,Sinv.Ahatcross.S),Matrix):Btr:=convert(Sinv.Bhat,Matrix):Ctr:=convert(map(simplify,Chat.S),Matrix):Dtr:=Dhat:'Atr'=convert(Atr,array);'Atrcross'=convert(Atrcross,array);Elementary factors:Allfactors:=map(simplify,MakeFactorization(Atr,Btr,Ctr,lambda)):afactors := Vector[row](dqw,0): for k from 1 to dqw do afactors[k]:=Allfactors[k](lambda): end do: print(`Elementary factors`=convert(afactors,array));Is the factorization correct?Wtest:=Factors2Transfer(Allfactors,lambda): simplify(Wtest(lambda)-W(lambda));

It is ready to hand to make one procedure which has one single argument W, a 2 x 2 rational matrix function and which outputs all factorizations for all calculated orderings.FactorizationFactors:=proc(W) local x, w1, Wtest,res1, ppoles, poles, pzeros, zeros, npoles, nzeros, mu, DimS, nmu, dw, dqw, k, muA, muZ, AllMOrdering, orderedA, orderedZ, r, Amin, Bmin, Cmin, Dmin, Amincross, kw, Allhat, TA, TZ, S, Sinv, Atr, Acrosstr, Btr, Ctr, Dtr, polesA, Allfactors, afactors, gammav, morder,mmorder; print('W(lambda)'=W(lambda)); w1:=Column(W(x),1): if not ((w1[1]=1) and (w1[2]=0)) then print(`The transfer function has not a correct form. Procedure will exit`); return(``);end if: DimS:= ColumnDimension(W(x)): ppoles:=x->expand(LCMDenomMatrixPolynom(W,x)): pzeros:=unapply(simplify(W(x)[2,2]*ppoles(x)),x): res1:=GetPolesandZeros(ppoles,pzeros):poles:=res1[1]: zeros:=res1[2]: mu:=res1[3]: npoles:=Dimension(poles): nzeros:=Dimension(zeros): nmu:=Dimension(mu): muA:=GetMultiplicity(poles,mu): muZ:=GetMultiplicity(zeros,mu): dw:=npoles: AllMOrdering:= GetAllMOrderings(muA,muZ,mu): print('h'=AllMOrdering[1]);morder:=AllMOrdering[2]: print(`Number of orderings`=morder); dqw:=AllMOrdering[1]-1+npoles: print('delta[q](W)'=dqw); r:=unapply(simplify(W(x)[1,2]*ppoles(x)),x): Amin:=Transpose(CompanionMatrix(ppoles(x),x)): Bmin:=Matrix(npoles,DimS,0): Bmin[npoles,DimS]:=1: Cmin:=Matrix(DimS,npoles): for k from 1 to npoles do Cmin[1,k]:=coeff(r(lambda),lambda,k-1): Cmin[2,k]:=coeff(pzeros(lambda)-ppoles(lambda),lambda,k-1): end do: Dmin:=IdentityMatrix(DimS,DimS): Amincross:=Transpose(CompanionMatrix(pzeros(x),x)): print(``);print(``);print('REALIZATION'); print('A'=convert(Amin,array));print('B'=convert(Bmin,array)); print('C'=convert(Cmin,array));print('D'=convert(Dmin,array)); print('A^x'=convert(Amincross,array)); gammav:=MakeGamma(dw,mu,true): Allhat:=QCmatrices(Amin,Bmin,Cmin,Dmin,dqw,gammav):print(` `); for mmorder from 1 to morder do print(`Ordering index`=mmorder); orderedA:=AllMOrdering[3][mmorder]: orderedZ:= AllMOrdering[4][mmorder]: print(`ordering poles`=convert(orderedA,list)); print(`ordering zeros`=convert(orderedZ,list)); TA:=MakeBasisA(orderedA,dqw,gammav): TZ:=MakeBasisZ(orderedZ,dqw,gammav): S:=UpperLowerTransformation(TA,TZ,dqw): Sinv:=MatrixInverse(S): Atr:=convert(map(simplify,Sinv.Allhat[1].S),Matrix): Btr:=convert(Sinv.Allhat[2],Matrix): Ctr:=convert(map(simplify,Allhat[3].S),Matrix): Allfactors[mmorder]:=map(simplify,MakeFactorization(Atr,Btr,Ctr,lambda)): print(` `);print(`FACTORIZATION`);afactors:=Vector[row](npoles,0): for k from 1 to npoles do afactors[k]:=Allfactors[mmorder][k](lambda): end do: print(`Elementary factors`=convert(afactors,array)); Wtest:=Factors2Transfer(Allfactors[mmorder],lambda): print(`Test`=simplify(Wtest(lambda)-W(lambda))); end do: return(Allfactors); end proc;

Define a rational matrix functionWj:=x-><<1,0>|<1/((x-1)^2*(x-2)^4*(x-4)^4*(x-5)^2),((x-3)*(x-4))/((x-2)*(x-5))>>:Wj(lambda);Now apply the routine. Caveat: it will take some time to get the factorizations for all available orderings!AllfactorsWj:=FactorizationFactors(Wj):