-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathCompanionMatrix_DMD.m
61 lines (41 loc) · 1.93 KB
/
CompanionMatrix_DMD.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
function [KEv, KModes,Norms] = CompanionMatrix_DMD( Data,eps_norm )
% Code is provided by Hassan Arbabi
% and modified by Keisuke Fujii
% Dynamic Mode Decomposition as presented by
% "Spectral analysis of nonlinear flows" by Rowley et al., Journal of FLuid
% Mechanics, 2009
% inputs :
% Data - the data matrix: each column is a a set of measurements done at
% each instant - the sampling frequency is assumed to be constant
% outputs:
% 1 - KModes- Koopman or Dynamic Modes: structures that scale exponentially
% with time - the exponent is given by Koopman or Dynamic Eigenvalues and
% could be imaginary as well
% 2- KEv - Koopman or Dynamic Eigenvalues: the exponents and frequencies
% that make up the time-varaiation of data in time
% 3- Norms - Euclidean (vector) norm of each mode, used to sort the data
X=Data(:,1:end-1);
% pseudo inverse approach
c = pinv(X)*Data(:,end);
m = size(Data,2)-1;
C = spdiags(ones(m,1),-1,m,m);
C(:,end)=c; %% companion matrix
[P,D,Z]=eig(full(C)); %% Ritz evalue and evectors
KEunsrtd = diag(D); %% Koopman Evalues
KMunsrtd = X*P; %% Koopman Modes
if nargout>2 % norm
Cnorm = sum(abs(KMunsrtd).^2,1); %% Euclidean norm of Koopman Modes
[Norms, ind]=sort(Cnorm,'descend'); %% sorting based on the norm
% ind = ind(abs(Norms)>eps_norm); % absolute tolerance
ind = ind(abs(Norms)/abs(Norms(1))>eps_norm); % relative tolerance
KEv= KEunsrtd(ind); %% Koopman Eigenvalues
KModes = KMunsrtd(:,ind); %% Koopman Modes
Norms = Norms(1:length(ind))' ;
end
end
%=========================================================================%
% Hassan Arbabi - 08-17-2015
% Mezic research group
% UC Santa Barbara
%=========================================================================%