How can i generate an Hankel matrix of the lorenz system?

hi, i'd like to know how i can create the Hankel matrix of the lorenz system to perform it over a DMD algorithm. Thank all for the help!!

[ANTONINO BIANCUZZO]

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) ;

Answer 1

[Nithin]

Hi Antonino,

I understand that you want to create the Hankel matrix of the Lorenz system, perform it over a DMD algorithm and then plot it.

To implement this, kindly refer to the following code snippet:

clc;
clear all;
close all;
% Define the Lorenz system function
lorenz = @(t,x,Beta) [Beta(1)*(x(2)-x(1)); x(1)*(Beta(2)-x(3))-x(2); x(1)*x(2)-Beta(3)*x(3)];
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;

% Plotting
figure;
surfl(real(Xdmd).');
xlabel('Time');
ylabel('DMD Mode');
zlabel('Amplitude');
title('DMD Reconstruction of Lorenz System');

For more information regarding “hankel” and “surf1”, kindly refer to the following MATLAB documentation:

https://www.mathworks.com/help/matlab/ref/hankel.html

https://www.mathworks.com/help/matlab/ref/surf.html

I hope it helps.

Regards,

Nithin Kumar.