Update: I built the matrix and I performed the DMD algorithm, but, as my final step in this code, i could not to plot the system yet. I attach the code afterwards: I tried two way without result
clear all, close all, clc
Beta = [10; 28; 8/3]; % chaotic values
x0 = [0; 1; 20]; % initial condition
dt = 0.001;
tspan = dt:dt:9; % 9.000 time span
%% ode options
options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,3));
[t,x] = ode45(@(t,x)lorenz(t,x,Beta),tspan,x0,options);
%% hankel matrix
x = x.';
n = 6000; % number of rows
m = 3000; % 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
r = 50;
[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;
surfl(real(Xdmd).');
% DMD reconstruction from "Data-driven spectral analysis for coordinative structures in ..."
%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)));% time dynamics
% Xdm(:,iter,p) = real(Phi(:,p)*(b(p).*exp(omega(p)*t(iter))));% reconstruct
% end
%end
%X_dmd = sum(Xdm,3) ;
