how to accelate the 'for loop' command' It's not end processing !!
4 次查看(过去 30 天)
显示 更早的评论
clear all, close all, clc
%import data from excel sheet I have attached
x = x.'; %just x variable
n =5000; % number of rows
m =2726; % number of columns
index1 = 1:n;
index2 = n:n+m-1;
X = []; Xprime=[];
for ir = 1:size(x,1)
% Hankel blocks ()
c = x(ir,index1).'; r = x(ir,index2);
H = hankel(c,r).';
c = x(ir,index1+1).'; r = x(ir,index2+1);
UH= hankel(c,r).';
X=[X,H]; Xprime=[Xprime,UH];
end
%% DMD of x
dt = 1;
t=[0:1:7725];
r = 16 ;
[U, S, V] = svd(X, 'econ');
r = min(r, size(U,2));
U_r = U(:, 1:r); % truncate to rank-r
S_r = S(1:r, 1:r);
V_r = V(:, 1:r);
Atilde = U_r' * Xprime * V_r / S_r; % low-rank dynamics
[W_r, D] = eig(Atilde);
Phi = Xprime * V_r / S_r * W_r; % DMD modes
%PhiP = U_r * V_r; % Projected DMD modes
lambda = diag(D); % discrete-time eigenvalues
omega = log(lambda)/dt; % continuous-time eigenvalues
% Compute DMD mode amplitudes b
x1 = X(:, 2);
b = Phi\x1; % equal to b = pinv(Phi)*x1
% DMD reconstruction
time_dynamics = zeros(r, length(t));
for iter = 1:length(t)
time_dynamics(:,iter) = (b.*exp(omega*t(iter)));
end
Xdmd = Phi * time_dynamics;
% DMD reconstruction
rr = length(lambda) ;
T = size(X,2) ;
time_dmd = zeros(T-1,rr);
for iter = 1:T-1
for p = 1:rr
time_dmd(iter,p) = b(p)*(exp(omega(p)*t(iter)));
Xdm(:,iter,p) = real(Phi(:,p)*(b(p).*exp(omega(p)*t(iter))));
end
end
X_dmd = sum(Xdm,3) ;
采纳的回答
DGM
2023-6-30
编辑:DGM
2023-6-30
Again, you're trying to construct large (~1.7 GB) arrays by growing them (in the reconstruction part). That will take an impractically long time (at least an hour or so). Like the other answer I gave, the loops are unnecessary anyway, so the whole issue can be avoided. Even on my dumpster computer with minimal ram, the loopless replacement section takes about 5 seconds.
As to whether the rest of the code works, I don't know. That's all over my head. I just tried to make sure my edits didn't change the results.
n = 5000; % number of rows
m = 2726; % number of columns
x = rand(1,n+m);
index1 = 1:n;
index2 = n:n+m-1;
% preallocate
szx = size(x,1);
X = zeros(m,n*szx);
Xprime = zeros(m,n*szx);
tic
for ir = 1:szx
% Hankel blocks ()
c = x(ir,index1).';
r = x(ir,index2);
H = hankel(c,r).';
c = x(ir,index1+1).';
r = x(ir,index2+1);
UH = hankel(c,r).';
cols = n*(ir-1)+1:n*ir;
X(:,cols) = H;
Xprime(:,cols) = UH;
end
toc
% preallocating would have been more important if x had more than one row
% if x never has more than one row, then the entire loop is unnecessary
%% DMD of x
dt = 1;
t = 0:1:7725;
r = 16 ;
[U, S, V] = svd(X, 'econ');
%%
r = min(r, size(U,2));
U_r = U(:, 1:r); % truncate to rank-r
S_r = S(1:r, 1:r);
V_r = V(:, 1:r);
Atilde = U_r' * Xprime * V_r / S_r; % low-rank dynamics
[W_r, D] = eig(Atilde);
Phi = Xprime * V_r / S_r * W_r; % DMD modes
%PhiP = U_r * V_r; % Projected DMD modes
lambda = diag(D); % discrete-time eigenvalues
omega = log(lambda)/dt; % continuous-time eigenvalues
% Compute DMD mode amplitudes b
x1 = X(:, 2);
b = Phi\x1; % equal to b = pinv(Phi)*x1
% DMD reconstruction
time_dmd = b.'.*exp(omega.'.*reshape(t(1:size(X,2)-1),[],1));
Xdm = permute(Phi,[3 2 1]).*time_dmd;
Xdm = real(permute(Xdm,[3 1 2]));
X_dmd = sum(Xdm,3);
0 个评论
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!