-
Notifications
You must be signed in to change notification settings - Fork 0
/
jordanmatris.m
47 lines (36 loc) · 1.07 KB
/
jordanmatris.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
function J = jordanmatris(A,tol)
if nargin < 2, tol = eps*1e7; end
[ev,mult] = heltalsev(A,tol)
jordan = zeros(size(A));
progress = 1; %How far the construction has come
for i=1:length(ev) %Loop over eigenvalues
j=1;
%Initialize
ri = []; pi = []; bi = []; ni = [];
while length(ri) < 2 | pi(end) ~= pi(end-1) %Compute table
ri(j) = rank((A - ev(i)*eye(length(A)))^j,tol);
pi(j) = length(A) - ri(j);
if j == 1
bi(j) = pi(j);
else
bi(j) = pi(j) - pi(j-1);
end
if j > 1
ni(j-1) = bi(j-1) - bi(j); %Fix, for last pos
end
j = j+1;
end
ni(j) = 0;
nValues = unique(ni);
for l=1:length(nValues)
positions = find(ni == nValues(l))
for m=1:nValues(l)
for k=1:length(positions)
n = positions(k);
jordan(progress:progress+n-1,progress:progress+n-1) = jordanblock(ev(i), n)
progress = progress + n;
end
end
end
J = jordan;
end