/
eigGet.m
89 lines (80 loc) · 2.44 KB
/
eigGet.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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
function [sn,un,cn,Ws,Wu,Wc] = eigGet(A,discrete)
% [sn,un,cn,Ws,Wu,Wc] = eigGet(A,discrete) ;
%
% Compute the eigenvalues and eigenvectors of the matrix A spanning the
% three local subspaces <Es,Eu,Ec> where A is MxM and s+u+c=M, which
% locally approximate the invariant manifolds <Ws,Wu,Wc>
%
% if discrete time system, discrete = 1
% if continuous time system, discrete = 0
%
% Shane Ross (revised 2.19.04)
sn=[];
un=[];
cn=[];
Ws=[];
Wu=[];
Wc=[];
% arbitrary small displacement for use in discrete case
delta = 1.e-4; % <== maybe needs tweeking?
[M,dum]=size(A);
[V,D] =eig(A); % obtain eigenvectors (V) and eigenvalues (D) of matrix A
V = cleanUpMatrix(V) ;
D = cleanUpMatrix(D) ;
s=0;u=0;c=0;
for k = 1:M
if discrete == 1 % discrete time system
if abs(D(k,k))-1 < -delta
s=s+1;
sn(s) = D(k,k);
jj=1; while abs(V(jj,k)) == 0, jj=jj+1; end
Ws(:,s) = V(:,k)./V(jj,k);% stable (s dimensional)
Ws(1:M,s) = RemoveInfinitesimals(Ws(:,s));
elseif abs(D(k,k))-1 > delta
u=u+1;
un(u) = D(k,k);
jj=1; while abs(V(jj,k)) == 0, jj=jj+1; end
Wu(:,u) = V(:,k)./V(jj,k);% unstable (u dimensional)
Wu(1:M,u) = RemoveInfinitesimals(Wu(:,u));
else
c=c+1;
cn(c) = D(k,k);
jj=1; while abs(V(jj,k)) == 0, jj=jj+1; end
Wc(:,c) = V(:,k)./V(jj,k);% center (c dimensional)
Wc(1:M,c) = RemoveInfinitesimals(Wc(:,c));
end
elseif discrete == 0 % continuous time system
if real(D(k,k)) < 0
s=s+1;
sn(s) = D(k,k);
jj=1; while abs(V(jj,k)) == 0, jj=jj+1; end
Ws(:,s) = V(:,k)./V(jj,k);% stable (s dimensional)
Ws(1:M,s) = RemoveInfinitesimals(Ws(:,s));
elseif real(D(k,k)) > 0
u=u+1;
un(u) = D(k,k);
jj=1; while abs(V(jj,k)) == 0, jj=jj+1; end
Wu(:,u) = V(:,k)./V(jj,k);% unstable (u dimensional)
Wu(1:M,u) = RemoveInfinitesimals(Wu(:,u));
else
c=c+1;
cn(c) = D(k,k);
jj=1; while abs(V(jj,k)) == 0, jj=jj+1; end
Wc(:,c) = V(:,k)./V(jj,k);% center (c dimensional)
Wc(1:M,c) = RemoveInfinitesimals(Wc(:,c));
end
end
end
end
function A=RemoveInfinitesimals(A)
% Remove all entries with absolute value smaller than TOL
TOL=1.e-14;
for k=1:length(A),
if abs(real(A(k))) < TOL,
A(k)=A(k) - real(A(k));
end
if abs(imag(A(k))) < TOL,
A(k)=A(k) - i*imag(A(k));
end
end
end