maple_mode(1); restart; with(linalg): Rotate:=proc(v,n) local s,j,v0,vv; v0:=convert(v,list); vv:=v0; s:=nops(v0); if (n>=0) then for j from 1 to s-n do vv[j]:=v0[j+n]; od; for j from s-n+1 to s do vv[j]:=v0[j+n-s]; od; RETURN(convert(vv,array)); else for j from -n+1 to s do vv[j]:=v0[j+n]; od; for j from 1 to -n do vv[j]:=v0[j+n+s]; od; RETURN(convert(vv,array)); fi; end ; diago:=(i,j)->if i=j then 1; else 0; fi; # Ce programme evalue le polynome, vu comme une liste l de # coefficients ranges par ordre croissant, en a par la methode de Horner. evalpolymat:=proc(l,a) local p,b,n,i; n:=nops(l);b:=l[n]; p:=[b]; i:=n-1; while(i<>0) do b:=evalm(l[i]+a*b);p:=[b,op(p)];i:=i-1;od; RETURN(evalm(b)); end: # Calcule les coefficents matriciels du polynome matrice adjointe de (x*I-A). fadeev:=proc(A) local Aj,AAj,Id,lmat,pcara,coef,n,i; n:=rowdim(A); Id:=matrix(n,n,diago); Aj:=Id; lmat:=[]; pcara:=[1]; for i from 1 to n do lmat:=[evalm(Aj),op(lmat)]; AAj:=evalm(map(normal,Aj&*A)); coef:=-trace(AAj)/i; pcara:=[coef,op(pcara)]; Aj:=evalm(map(normal,AAj&+coef&*Id)); od; RETURN(lmat,pcara); end: # Prend une liste de matrices rangee par ordre de puissances # croissantes, calcule la derivee derive__listemat:=proc(liste) local i,n,l,j,p; l:=liste; n:=nops(l);p:=[]; for i from 1 to n do l[i]:=evalm((i-1)*l[i]); od; for j from 1 to n-1 do p:= [op(p),l[j+1]]; od; RETURN(p); end: # Prend une liste de racines ranges par "sort" et renvoie une liste de # listes de la forme [racine,multiplicite], comme "roots". arrangement__racines:=proc(l) local i,n,v,j,u,l1; n:=nops(l); l1:=[op(l),0]; v:=[]; i:=0; while(i<=n-2) do i:=i+1; j:=1; u:=[]; while(l1[i]=l1[i+1]) and (i<=n-1) do j:=j+1; i:=i+1; od; u:=[l1[i],j]; v:=[op(v),u]; od; if (l[n-1]<>l[n]) then v:=[op(v),[l[n],1]];fi; RETURN(v); end: # Prend une racine et sa multiplicite aisi que la liste des coeffs # matriciels de B(x) pour creer la matrice colonne des # B^(i)(racine)/i! collees les unes sous les autres. colonneB:=proc(v,Bliste) local g,l,C,i,A; g:=Bliste; l:=[]; for i from 0 to (v[2]-1) do l:=[op(l),simplify(evalpolymat(g,v[1]))]; g:=derive__listemat(g)/(i+1); od; A:=blockmatrix(v[2],1,l); RETURN(evalm(A)); end: # Arguments: resultat de fadeev # Renvoie une liste de listes de la forme # [[racine, multiplicite],matrice] ou la matrice est calulee par colonneB. construction__colonneB:=proc(Bliste,pliste) local u,f,v,l,g,A,i,n,j,B; n:=nops(Bliste); f:=x->sum('pliste[k]*x^(k-1)','k'=1..n+1); u:=sort([solve(f(x)=0,x)]); u:=arrangement__racines(u); B:=[]; for j from 1 to nops(u) do v:=u[j]; B:=[op(B),[v,colonneB(v,Bliste)]]; od; RETURN(B); end: # teste si la ligne i de b est nulle test__ligne__nulle:=proc(b,i) local k,n; n:=rowdim(b); k:=1; while(k<=n) do if (b[i,k]<>0) then RETURN(0);break; else k:=k+1;fi; od; if (k=n+1) then RETURN(1);fi; end: # decalage des lignes vers la gauche decalage__ligne:=proc(B,i) local n,k,l,m; n:=rowdim(B); m:=coldim(B); for l from 1 to (m/n-1) do for k from 1 to n do B[i,m-l*n+k]:=B[i,m-(l+1)*n+k]; od;od; RETURN(evalm(B)); end: # enleve les #lignes premieres colonnes de B coupe__matrice:=proc(B) local n; n:=rowdim(B); delcols(B,1..n); end: cycle:=proc(B,j) local i,k,n,m,l,Resultat,u; m:=coldim(B); n:=rowdim(B); l:=row(B,j);Resultat:=[]; for k from 0 to (m/n-1) do u:=[]; for i from 1 to n do u:=[op(u),l[n*k+i]];od; Resultat:=[op(Resultat),u] ; od; RETURN(Resultat); end: # construction permet de construire la matrice de passage a partir # de la liste des cycles. construction:=proc(l,q) local i,j,k,n,m,p,v,u,B,s; n:=nops(l); B:=matrix(q,1,[0$q]); for i from 1 to n do v:=l[i][2]; m:=nops(v); for j from 1 to m do u:=v[j]; p:=nops(u); for k from 1 to p do B:= concat(matrix(q,1,u[k]),B); od; od; od; s:=coldim(B); B:=delcols(B,s..s); RETURN(evalm(B)); end: test1:=proc(cycle1,C,i) local g,j,k,l,A,v,u,f,n; g:=nops(cycle1); n:=rowdim(C); u:=[]; if (g=0) then RETURN(1); else do v:=[]; for j from 1 to n do v:=[op(v),C[i,j]]; od; u:=[matrix(1,nops(v),v)]; for k from 1 to g do l:=cycle1[k][1]; u:=[op(u),matrix(1,nops(l),l)]; od; A:=blockmatrix(nops(u),1,u); f:=rank(A); if(f= numero__ligne) then cycle1:=[op(cycle1),cycle(C,k)]; taille__cycle:=taille__cycle + nops(cycle(C,k)); numero__ligne:=numero__ligne+1 ; fi; C:=decalage__ligne(C,k); fi; od; if(coldim(C)>rowdim(C)) then C:=coupe__matrice(C);fi; od; l:=[op(l),v[1][1],cycle1]; res:=[op(res),l]; RETURN(res); end: # Forme normale de Jordan (si polynome caracteristique scinde) TER__Jordan:=proc(A) local Bliste,pliste,B,n,i,v,C,k,j,l,taille__cycle,cycle1,Resultat,P,Q,J,numero__ligne; (Bliste,pliste):=fadeev(A); B:=construction__colonneB(Bliste,pliste); n:=nops(B);Resultat:=[]; for i from 1 to n do Resultat:=dedans(B[i],Resultat); od; P:=construction(Resultat,rowdim(A)); Q:=P&^(-1); J:=Q&*A&*P; RETURN(evalm(P),evalm(J)); end: # Exemples A:=matrix(3,3,[3,-1,1,2,0,1,1,-1,2]); TER__Jordan(A); B:=matrix(3,3,[3,2,-2,-1,0,1,1,1,0]); TER__Jordan(B); ex1:=BlockDiagonal(JordanBlock(2,2),JordanBlock(2,1),JordanBlock(1,4),JordanBlock(4,1));P:=randmatrix(8,8); ex2:=evalm(P&^(-1)&*ex1&*P): TER__Jordan(ex2); # traduction transforme un polynome de matrices vu comme une liste # de matrices coefficients en une matrice de polynome. traduction:=proc(Bliste) local C,n,i,j; n:=coldim(Bliste[1]); C:=JordanBlock(1,n); for i from 1 to n do for j from 1 to n do C[i,j]:=sum('Bliste[k][i,j]*x^(k-1)','k'=1..n); od; od; RETURN(evalm(C)); end: # nouvelle__ecriture trouve les coefficient de B dans la base des Q^k, # seulement pour k1) then v := l[i][2]; m := nops(v); for j to m do u := v[j]; p := nops(u); for k from 1 to p do w := u[k]; h := nops(w); for o from 1 to h do B := concat(B,matrix(q,1,w[o])); od; od; od; else v:=l[i][2]; m:=nops(v); for j to m do u := v[j]; p := nops(u); for k from 1 to p do B:=concat(B,matrix(q,1,u[k])); od; od; fi; od; if coldim(B)=1 then RETURN([]) fi; B := delcols(B,1 .. 1); RETURN(evalm(B)); end: # fabriq__cycles prend un cycle de Q(A) et ajoute les puissances de A # des vecteurs pour renvoyer un cycle fabriq__cycles:=proc(A,ccycle,d) local n,i,C,j,k,nouveaux__cycles; n:=nops(ccycle); nouveaux__cycles:=ccycle; if(d>1) then for j from 1 to n do nouveaux__cycles[j]:=[[op(nouveaux__cycles[j])],evalm(A&*nouveaux__cycles[j])]; if (d>3) then for k from 3 to d do nouveaux__cycles[j]:=[op(nouveaux__cycles[j]),evalm(A&*nouveaux__cycles[j][k-1])]; od;fi; od; fi; RETURN(nouveaux__cycles); end: suite__A:=proc(nouveaux__cycles) local i,n,m,M,colonne,p,j,nulle; nulle:=(i,j)->0; n:=nops(nouveaux__cycles); m:=nops(nouveaux__cycles[1]); p:=nops(nouveaux__cycles[1][1]); M:=matrix(1,m,nulle); for i from 1 to n do colonne[i]:=blockmatrix(1,m,nouveaux__cycles[i]); od; for j from 1 to n do M:=blockmatrix(2,1,[M,colonne[j]]); od; M:=delrows(M,1..1); RETURN(evalm(transpose(M))); end: test:=proc(cycle1,C,n,i,d2) local v,j,A,k,u,f,g,l; g:=nops(cycle1); if (g=0) then RETURN(1); else do v:=[]; for j from 1 to n do v:=[op(v),C[i,j]]; od; k:=1; while(k<=g) do l:=cycle1[k][1]; if (d2=1) then u:=[l,v]; else u:=[op(l),v]; fi; A:=blockmatrix(1,nops(u),u); f:=rank(A); if(f=numero__ligne) then if(test(cycle1,C,n,i,d2)=1 and taille__cyclen) then C:=coupe__matrice(C);fi; od; l:=[[f,d],cycle1]; RETURN(l); end: corps:=proc(A) local Bliste,pliste,f,i,ccycle,Resultat,u ,P,Q,J ; (Bliste,pliste):=fadeev(A); Resultat:=[]; f:=sum('pliste[k]*x^(k-1)','k'=1..nops(pliste)); u:=factors(f)[2]; for i from 1 to nops(u) do ccycle:=demarrage(A,Bliste,u[i][1],u[i][2]); Resultat:=[op(Resultat),ccycle]; od; P:=construction__special(Resultat,rowdim(A)); Q:=P&^(-1); J:=Q&*A&*P; RETURN(evalm(P),evalm(J)); end: final:=proc(A) local Bliste,pliste,f,i,ccycle,Resultat1,Resultat2,u,k ,P,P1,P2,Q,J,B ; (Bliste,pliste):=fadeev(A); Resultat1:=[]; Resultat2:=[]; f:=sum('pliste[k]*x^(k-1)','k'=1..nops(pliste)); u:=factors(f)[2]; for i from 1 to nops(u) do if (degree(u[i][1])>1) then ccycle:=demarrage(A,Bliste,u[i][1],u[i][2]); Resultat1:=[op(Resultat1),ccycle]; else k:=solve(u[i][1]=0); B:=[[k,u[i][2]],colonneB([k,u[i][2]],Bliste)]; Resultat2:=dedans(B,Resultat2); fi; od; P1:=construction__special(Resultat1,rowdim(A)); if(nops(Resultat2)>0) then P2:=construction(Resultat2,rowdim(A)); if P1=[] then P:=P2 else P:=concat(P1,P2); fi; else P:=P1; fi; Q:=P&^(-1); J:=Q&*A&*P; RETURN(Resultat1,Resultat2,evalm(P),evalm(J)); end: coefficients:=proc(f) local n,i,l; n:=degree(f); l:=tcoeff(f); for i from 1 to n do l:=[op(l),coeff(f,x^i)]; od; RETURN(l); end: # coeffq=list of coefficients of q, degq=degree of q, A=matrice # k no du bloc de vecteur a calculer (numerote a partir de 0) # v2 la matrice des vecteurs trouves, par blocs annexe:=proc(A,v2,k,degq,coeffq) local i,vk0,vki,l,m; vk0:=evalm(sum('coeffq[l+1]*sum(binomial(l,m)*v2[k-m+1][l-m+1],m=1..min(l,k))',l=1..degq)); # replace vk0 by preimage by A^degq: shift by degq coefficients vk0:=evalm(Rotate(vk0,-degq)); l:=[evalm(vk0)]; for i from 1 to degq-1 do vki:=evalm(A&^i&*vk0-sum('binomial(i,m)*v2[k-m+1,i-m+1]','m'=1..min(i,k))); l:=[op(l),evalm(vki)]; od; RETURN(l); end: ratjord:=proc(P,x,n) local C,s,M,suite,N,j; C:=companion(P,x); s:=coldim(C); M:=matrix(s,s,0); N:=copy(M); M[1,s]:=1; suite:=NULL; for j from 1 to n-1 do suite:=suite,N$(j-1),C,M,N$(n-1-j); od; suite:=suite,N$(n-1),C; RETURN(blockmatrix(n,n,[suite])); end: Jordan2:=proc(A) local P1,P2,P3,J,Q,l,liste,Resultat2,P,B,coeffq,degq,n,i,p,q,v,v1,v2,v3,v4,m,j,k,n1,n2,passage,tmp1,tmp2,tmp3,Atmp,ii,jj,kk; (l,Resultat2,P,B):=final(A); n:=nops(l); for i from 1 to n do v:=l[i]; degq:=degree(v[1][1],x); coeffq:=coefficients(v[1][1]); v1:=v[2]; p:=nops(v1); # nombre de blocs de Jordan for j from 1 to p do n1:=1; v2:=v1[j]; q:=nops(v2); # nombre de cycles du bloc Atmp:=ratjord(v[1][1],x,q); passage := [0$(q*degq)]; tmp1 := passage ; tmp2 := passage ; for k from 1 to q*degq do tmp2:=tmp1; tmp2[k]:=1; passage[k]:=tmp2; od; # passage:=convert(matrix(q*degq,q*degq,diago),list); tmp1:=NULL; # sauve en ligne la 1ere matrice de passage # remplace v2 par vect. base canonique for k from 1 to q do tmp1:=tmp1,op(v2[k]); n2:=n1+nops(v2[k]); v2[k]:=passage[n1..(n2-1)]; n1:=n2; od; tmp2:=op(v2[1]); # contient en ligne la 2eme matrice de passage for k from 2 to q do liste:=annexe(Atmp,v2,k-1,degq,coeffq); v2[k]:=liste; tmp2:=tmp2,op(v2[k]); od; # produit de matrices tmp3:=[tmp2]*[tmp1] tmp2:=[tmp2]; tmp1:=[tmp1]; tmp3:=tmp2; k:=nops(tmp1); for ii from 1 to nops(tmp2) do tmp3[ii]:=[0$nops(tmp1[1])]; for kk from 1 to nops(tmp1[1]) do tmp4:=0; for jj from 1 to k do # print(ii,jj,kk); tmp4:=tmp4+tmp2[ii][jj]*tmp1[jj][kk]; od; tmp3[ii,kk]:=tmp4; od; od; # remplace dans v1 n1:=1; for k from 1 to q do n2:=n1+nops(v2[k]); v2[k]:=tmp3[n1..(n2-1)]; n1:=n2; od; v1[j]:=v2; od; v[2]:=v1; l[i]:=v; od; P1:=construction__special(l,rowdim(A)); if (nops(Resultat2)>0) then P2:=construction(Resultat2,rowdim(A)); if P1=[] then P3:=P2; else P3:=concat(P1,P2); fi; else P3:=P1; fi; Q:=P3&^(-1); J:=Q&*A&*P3; RETURN(evalm(P3),evalm(J)); end: # Exemples A:=ratjord(x^2-2,x,2); Jordan2(A); A:=matrix(6,6,[1,-2,4,-2,5,-4,0,1,5/2,-7/2,2,-5/2,1,-5/2,2,-1/2,5/2,-3,0,-1,9/2,-7/2,3,-7/2,0,0,2,-2,3,-1,1,-3/2,-1/2,1,3/2,1/2]); Jordan2(A); ex1:=BlockDiagonal(ratjord(x^2-2,x,2),ratjord(x^2+2,x,2));P:=randmatrix(rowdim(ex1),rowdim(ex1)); ex2:=evalm(P&^(-1)&*ex1&*P): Jordan2(ex2);