1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | |
---|
3 | version="$Id: jordan.lib,v 1.9 1999-02-18 19:29:01 mschulze Exp $"; |
---|
4 | info=" |
---|
5 | LIBRARY: jordan.lib PROCEDURES TO COMPUTE THE JORDAN NORMAL FORM |
---|
6 | AUTHOR: Mathias Schulze, email: mschulze@mathematik.uni-kl.de |
---|
7 | |
---|
8 | jordan(M[,opt]); eigenvalues, Jordan block sizes, transformation matrix of M |
---|
9 | jordanmatrix(l); Jordan matrix with eigenvalues l[1], Jordan block sizes l[2] |
---|
10 | jordanform(M); Jordan normal form of constant square matrix M |
---|
11 | inversemat(M); inverse matrix of invertible constant matrix M |
---|
12 | "; |
---|
13 | |
---|
14 | LIB "ring.lib"; |
---|
15 | /////////////////////////////////////////////////////////////////////////////// |
---|
16 | |
---|
17 | proc jordan(matrix M,list #) |
---|
18 | "USAGE: jordan(M[,opt]); M constant square matrix, opt integer |
---|
19 | ASSUME: the eigenvalues of M are in the coefficient field |
---|
20 | RETURN: a list l with 3 entries of type ideal, list of intvecs, matrix with |
---|
21 | l[1] : eigenvalues of M |
---|
22 | l[2][i][j] : size of j-th Jordan block with eigenvalue l[1][i] |
---|
23 | l[3] : transformation matrix U with U^(-1)*M*U in Jordan normal form |
---|
24 | which entries of l are computed depends on opt: |
---|
25 | opt=-1 : compute l[1] |
---|
26 | opt= 0 : compute l[1] and l[2] (default) |
---|
27 | opt= 1 : compute l[1],l[2] and l[3] |
---|
28 | opt= 2 : compute l[1] and l[3] |
---|
29 | NOTE: a non constant polynomial matrix M is replaced by its constant part |
---|
30 | EXAMPLE: example jordan; shows an example |
---|
31 | " |
---|
32 | { |
---|
33 | // test if square matrix |
---|
34 | int n=nrows(M); |
---|
35 | if(n!=ncols(M)) |
---|
36 | { |
---|
37 | "//no square matrix"; |
---|
38 | return(); |
---|
39 | } |
---|
40 | |
---|
41 | // get constant part |
---|
42 | def br=basering; |
---|
43 | map zero=br,0; |
---|
44 | M=zero(M); |
---|
45 | kill zero; |
---|
46 | |
---|
47 | // change to polynomial ring for factorization |
---|
48 | changeord("pr","dp"); |
---|
49 | matrix M=imap(br,M); |
---|
50 | |
---|
51 | // factorize characteristic polynomial |
---|
52 | list l=factorize(det(M-var(1)*freemodule(n)),2); |
---|
53 | |
---|
54 | // get multiplicities mM |
---|
55 | def eM,mM=l[1..2]; |
---|
56 | kill l; |
---|
57 | |
---|
58 | // test if eigenvalues in the coefficient field |
---|
59 | int i; |
---|
60 | for(i=size(eM);i>=1;i--) |
---|
61 | { |
---|
62 | if(deg(eM[i])>1) |
---|
63 | { |
---|
64 | kill pr; |
---|
65 | "//eigenvalues not in the coefficient field"; |
---|
66 | return(); |
---|
67 | } |
---|
68 | } |
---|
69 | |
---|
70 | // get eigenvalues eM |
---|
71 | map inv=pr,-var(1); |
---|
72 | eM=simplify(inv(eM),1); |
---|
73 | setring br; |
---|
74 | map zero=pr,0; |
---|
75 | ideal eM=zero(eM); |
---|
76 | kill pr; |
---|
77 | |
---|
78 | // sort eigenvalues |
---|
79 | int j; |
---|
80 | poly e; |
---|
81 | int m; |
---|
82 | for(i=size(eM);i>=2;i--) |
---|
83 | { |
---|
84 | for(j=i-1;j>=1;j--) |
---|
85 | { |
---|
86 | if(eM[i]<eM[j]) |
---|
87 | { |
---|
88 | e=eM[i]; |
---|
89 | eM[i]=eM[j]; |
---|
90 | eM[j]=e; |
---|
91 | m=mM[i]; |
---|
92 | mM[i]=mM[j]; |
---|
93 | mM[j]=m; |
---|
94 | } |
---|
95 | } |
---|
96 | } |
---|
97 | kill e,m; |
---|
98 | |
---|
99 | // get option parameter |
---|
100 | int opt=0; |
---|
101 | if(size(#)>0) |
---|
102 | { |
---|
103 | if(typeof(#[1])=="int") |
---|
104 | { |
---|
105 | opt=#[1]; |
---|
106 | } |
---|
107 | } |
---|
108 | if(opt<0) |
---|
109 | { |
---|
110 | return(list(eM)); |
---|
111 | } |
---|
112 | |
---|
113 | // define needed variables |
---|
114 | int k,l; |
---|
115 | matrix E=freemodule(n); |
---|
116 | matrix Mi,Ni; |
---|
117 | module sNi; |
---|
118 | list K; |
---|
119 | if(opt>=1) |
---|
120 | { |
---|
121 | module V,K1,K2; |
---|
122 | matrix v[n][1]; |
---|
123 | } |
---|
124 | if(opt<=1) |
---|
125 | { |
---|
126 | list bM; |
---|
127 | intvec bMi; |
---|
128 | } |
---|
129 | |
---|
130 | // do the following for all eigenvalues eM[i] |
---|
131 | for(i=ncols(eM);i>=1;i--) |
---|
132 | { |
---|
133 | Mi=M-eM[i]*E; |
---|
134 | |
---|
135 | // compute kernels K of powers of Mi |
---|
136 | K=list(module()); |
---|
137 | for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi) |
---|
138 | { |
---|
139 | sNi=syz(Ni); |
---|
140 | K=K+list(sNi); |
---|
141 | } |
---|
142 | |
---|
143 | if(opt<=1) |
---|
144 | { |
---|
145 | // compute Jordan block size vector corresponding to eigenvalue eM[i] |
---|
146 | bMi=0; |
---|
147 | bMi[size(K[2])]=0; |
---|
148 | for(j=size(K);j>=2;j--) |
---|
149 | { |
---|
150 | for(k=size(bMi);k>size(bMi)+size(K[j-1])-size(K[j]);k--) |
---|
151 | { |
---|
152 | bMi[k]=bMi[k]+1; |
---|
153 | } |
---|
154 | } |
---|
155 | bM=list(bMi)+bM; |
---|
156 | } |
---|
157 | |
---|
158 | if(opt>=1) |
---|
159 | { |
---|
160 | // compute generating vectors for Jordan basis vectors corresponding to |
---|
161 | // eigenvalue eM[i] |
---|
162 | if(size(K)>1) |
---|
163 | { |
---|
164 | for(j,K1=2,0;j<=size(K)-1;j++) |
---|
165 | { |
---|
166 | K2=K[j]; |
---|
167 | K[j]=interred(reduce(K[j],std(K1+module(Mi*K[j+1])))); |
---|
168 | K1=K2; |
---|
169 | } |
---|
170 | K[j]=interred(reduce(K[j],std(K1))); |
---|
171 | } |
---|
172 | |
---|
173 | // compute Jordan basis vectors corresponding to eigenvalue eM[i] from |
---|
174 | // generating vectors |
---|
175 | for(j=size(K);j>=2;j--) |
---|
176 | { |
---|
177 | for(k=size(K[j]);k>=1;k--) |
---|
178 | { |
---|
179 | v=K[j][k]; |
---|
180 | for(l=j;l>=1;l--) |
---|
181 | { |
---|
182 | V=module(v)+V; |
---|
183 | v=Mi*v; |
---|
184 | } |
---|
185 | } |
---|
186 | } |
---|
187 | } |
---|
188 | } |
---|
189 | kill l; |
---|
190 | |
---|
191 | // create return list |
---|
192 | list l=eM; |
---|
193 | if(opt<=1) |
---|
194 | { |
---|
195 | l[2]=bM; |
---|
196 | } |
---|
197 | if(opt>=1) |
---|
198 | { |
---|
199 | l[3]=V; |
---|
200 | } |
---|
201 | |
---|
202 | return(l); |
---|
203 | } |
---|
204 | example |
---|
205 | { "EXAMPLE:"; echo=2; |
---|
206 | ring R=0,x,dp; |
---|
207 | matrix M[3][3]=3,2,1,0,2,1,0,0,3; |
---|
208 | print(M); |
---|
209 | jordan(M); |
---|
210 | } |
---|
211 | /////////////////////////////////////////////////////////////////////////////// |
---|
212 | |
---|
213 | proc jordanmatrix(list l) |
---|
214 | "USAGE: jordanmatrix(l); l list of ideal and list of intvecs |
---|
215 | RETURN: the Jordan matrix J defined by l: |
---|
216 | l[1] : eigenvalues of J |
---|
217 | l[2][i][j] : size of j-th Jordan block with eigenvalue l[1][i] |
---|
218 | EXAMPLE: example jordanmatrix; shows an example |
---|
219 | " |
---|
220 | { |
---|
221 | // check parameters |
---|
222 | if(size(l)<2) |
---|
223 | { |
---|
224 | "//not enough entries in argument list"; |
---|
225 | return(); |
---|
226 | } |
---|
227 | def eJ,bJ=l[1..2]; |
---|
228 | kill l; |
---|
229 | if(typeof(eJ)!="ideal") |
---|
230 | { |
---|
231 | "//first entry in argument list not an ideal"; |
---|
232 | return(); |
---|
233 | } |
---|
234 | if(typeof(bJ)!="list") |
---|
235 | { |
---|
236 | "//second entry in argument list not a list"; |
---|
237 | return(); |
---|
238 | } |
---|
239 | if(size(eJ)<size(bJ)) |
---|
240 | { |
---|
241 | int s=size(eJ); |
---|
242 | } |
---|
243 | else |
---|
244 | { |
---|
245 | int s=size(bJ); |
---|
246 | } |
---|
247 | |
---|
248 | // get size of Jordan matrix |
---|
249 | int i,j,k,n; |
---|
250 | for(i=s;i>=1;i--) |
---|
251 | { |
---|
252 | if(typeof(bJ[i])!="intvec") |
---|
253 | { |
---|
254 | "//second entry in argument list not a list of integer vectors"; |
---|
255 | return(); |
---|
256 | } |
---|
257 | else |
---|
258 | { |
---|
259 | for(j=size(bJ[i]);j>=1;j--) |
---|
260 | { |
---|
261 | k=bJ[i][j]; |
---|
262 | if(k>0) |
---|
263 | { |
---|
264 | n=n+k; |
---|
265 | } |
---|
266 | } |
---|
267 | } |
---|
268 | } |
---|
269 | matrix J[n][n]; |
---|
270 | |
---|
271 | // create Jordan matrix |
---|
272 | int l; |
---|
273 | for(i,l=1,1;i<=s;i++) |
---|
274 | { |
---|
275 | for(j=1;j<=size(bJ[i]);j++) |
---|
276 | { |
---|
277 | k=bJ[i][j]; |
---|
278 | if(k>0) |
---|
279 | { |
---|
280 | while(k>=2) |
---|
281 | { |
---|
282 | J[l,l]=eJ[i]; |
---|
283 | J[l,l+1]=1; |
---|
284 | k,l=k-1,l+1; |
---|
285 | } |
---|
286 | J[l,l]=eJ[i]; |
---|
287 | l++; |
---|
288 | } |
---|
289 | } |
---|
290 | } |
---|
291 | |
---|
292 | return(J); |
---|
293 | } |
---|
294 | example |
---|
295 | { "EXAMPLE:"; echo=2; |
---|
296 | ring R=0,x,dp; |
---|
297 | list l; |
---|
298 | l[1]=ideal(2,3); |
---|
299 | l[2]=list(intvec(1),intvec(2)); |
---|
300 | print(jordanmatrix(l)); |
---|
301 | } |
---|
302 | /////////////////////////////////////////////////////////////////////////////// |
---|
303 | |
---|
304 | proc jordanform(matrix M) |
---|
305 | "USAGE: jordanform(M); M constant square matrix |
---|
306 | ASSUME: the eigenvalues of M are in the coefficient field |
---|
307 | RETURN: the Jordan normal form of M |
---|
308 | NOTE: a non constant polynomial matrix M is replaced by its constant part |
---|
309 | EXAMPLE: example jordanform; shows an example |
---|
310 | " |
---|
311 | { |
---|
312 | return(jordanmatrix(jordan(M))); |
---|
313 | } |
---|
314 | example |
---|
315 | { "EXAMPLE:"; echo=2; |
---|
316 | ring R=0,x,dp; |
---|
317 | matrix M[3][3]=3,2,1,0,2,1,0,0,3; |
---|
318 | print(M); |
---|
319 | print(jordanform(M)); |
---|
320 | } |
---|
321 | /////////////////////////////////////////////////////////////////////////////// |
---|
322 | |
---|
323 | proc inversemat(matrix M) |
---|
324 | "USAGE: inversemat(M); M constant square matrix |
---|
325 | ASSUME: M is invertible |
---|
326 | RETURN: the inverse matrix of M |
---|
327 | NOTE: a non constant polynomial matrix M is replaced by its constant part |
---|
328 | EXAMPLE: example inversemat; shows an example |
---|
329 | " |
---|
330 | { |
---|
331 | // test if square matrix |
---|
332 | int n=nrows(M); |
---|
333 | if(n!=ncols(M)) |
---|
334 | { |
---|
335 | "//no square matrix"; |
---|
336 | return(); |
---|
337 | } |
---|
338 | |
---|
339 | // get constant part |
---|
340 | def br=basering; |
---|
341 | map zero=br,0; |
---|
342 | M=zero(M); |
---|
343 | |
---|
344 | // compute inverse |
---|
345 | return(lift(M,freemodule(n))); |
---|
346 | } |
---|
347 | example |
---|
348 | { "EXAMPLE:"; echo=2; |
---|
349 | ring R=0,x,dp; |
---|
350 | matrix M[3][3]=3,2,1,0,2,1,0,0,3; |
---|
351 | print(M); |
---|
352 | print(inversemat(M)); |
---|
353 | } |
---|
354 | /////////////////////////////////////////////////////////////////////////////// |
---|
355 | |
---|