Add files via upload

Brusselator_Section.m

@@ -0,0 +1,197 @@
+% ------------------------------------------------------------------
+% SINDy method for discovering mappings in Poincaré sections 
+% ------------------------------------------------------------------
+% Application to the driven Brusselator
+%
+%           x' = a + alpha*sin(2*pi*t/T) - (b+1)*x + x^2*y
+%           y' = b*x - x^2*y     
+%
+% Here a,b,alpha are real-valued parameters and T > 0 controls the 
+% period of forcing.
+%
+% This code is associated with the paper 
+% "Poincaré maps for multiscale physics discovery and nonlinear Floquet
+% theory" by Jason J. Bramburger and J. Nathan Kutz (Physica D, 2020). 
+% This script is used to obtain the results in Section 3.4.
+% ------------------------------------------------------------------
+
+% Clean workspace
+clear all
+close all 
+clc
+
+format long
+
+%Model parameters
+a = 0.4;
+b = 1.2;
+T = 1; 
+alpha = 0.1;
+
+%ODE generation parameters
+m = 3; %Dimension of ODE
+n = m-1; %Dimension of Poincaré section
+dt = 0.005;
+tspan = (0:100000-1)*dt;
+options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,m));
+
+%Generate Trajectories
+x0(1,:) = [0; 0; 0]; 
+[~,xdat(1,:,:)]=ode45(@(t,x) Brusselator(x,a,b,T,alpha),tspan,x0(1,:),options);
+kfinal = 2;
+if kfinal >= 2
+    for k = 2:kfinal
+        x0(k,:) = [4*rand; 4*rand; 0]; %Initial conditions start in section 
+        [~,xdat(k,:,:)]=ode45(@(t,x) Brusselator(x,a,b,T,alpha),tspan,x0(k,:),options);
+    end
+end
+
+%% Poincaré section data
+
+%Counting parameter
+count = 1;
+
+%Initialize
+Psec = [];
+PsecNext = [];
+
+%Create Poincaré section data
+for i = 1:kfinal
+    for j = 1:length(xdat(i,:,1))-1 
+        if  (mod(xdat(i,j,3),T) >= T-dt && mod(xdat(i,j+1,3),T) <= dt) %&& j >= length(xdat(i,:,1))/50) 
+            temp(count,:) = xdat(i,j+1,1:2); %nth iterate
+            count = count + 1;
+        end
+    end
+    Psec = [Psec; temp(1:length(temp)-1,:)];
+    PsecNext = [PsecNext; temp(2:length(temp),:)];
+   	count = 1;
+    temp = [];
+end
+
+%% SINDy for Poincaré Sections
+
+% Access SINDy directory
+addpath Util
+
+% Create the recurrence data
+xt = Psec;
+xtnext = PsecNext;
+
+% pool Data  (i.e., build library of nonlinear time series)
+polyorder = 5; %polynomial order 
+usesine = 0; %use sine on (1) or off (0)
+
+Theta = poolData(xt,n,polyorder,usesine);
+
+% compute Sparse regression: sequential least squares
+lambda = 0.01;      % lambda is our sparsification knob.
+
+% apply iterative least squares/sparse regression
+Xi = sparsifyDynamicsAlt(Theta,xtnext,lambda,n);
+if n == 4
+[yout, newout] = poolDataLIST({'x','y','z','w'},Xi,n,polyorder,usesine);
+elseif n == 3
+[yout, newout] = poolDataLIST({'x','y','z'},Xi,n,polyorder,usesine);
+elseif n == 2
+ [yout, newout] = poolDataLIST({'x','y'},Xi,n,polyorder,usesine);
+elseif n == 1 
+  [yout, newout] = poolDataLIST({'x'},Xi,n,polyorder,usesine);
+end 
+
+fprintf('SINDy model: \n ')
+for k = 2:size(newout,2) 
+    SINDy_eq = newout{1,k}; 
+    SINDy_eq = [SINDy_eq  ' = '];
+    new = 1;
+   for j = 2:size(newout, 1)
+       if newout{j,k} ~= 0 
+           if new == 1 
+             SINDy_eq = [SINDy_eq  num2str(newout{j,k}) newout{j,1} ];  
+             new = 0;
+           else 
+             SINDy_eq = [SINDy_eq  ' + ' num2str(newout{j,k}) newout{j,1} ' '];
+           end 
+       end
+   end
+  fprintf(SINDy_eq)
+  fprintf('\n ')
+end 
+
+%% Simulate SINDy Map
+
+a = zeros(1000,1); %SINDy map solution
+b = zeros(1000,1);
+a(1) = Psec(1,1);
+b(1) = Psec(1,2);
+
+for k = 1:999
+
+    % Constant terms
+    a(k+1) = Xi(1,1);
+    b(k+1) = Xi(1,2);
+
+    %Polynomial terms
+   for p = 1:polyorder
+       for j = 0:p
+           a(k+1) = a(k+1) + Xi(1 + j + p*(p+1)/2,1)*(a(k)^(p-j))*(b(k)^j);
+           b(k+1) = b(k+1) + Xi(1 + j + p*(p+1)/2,2)*(a(k)^(p-j))*(b(k)^j);
+       end
+   end
+
+   if usesine == 1
+        a(k+1) = a(k+1) + Xi((p+1)*p/2+p+2,1)*sin(a(k)) + Xi((p+1)*p/2+p+3,1)*sin(b(k))+ Xi((p+1)*p/2+p+4,1)*cos(a(k)) + Xi((p+1)*p/2+p+5,1)*cos(b(k));
+        b(k+1) = b(k+1) + Xi((p+1)*p/2+p+2,2)*sin(a(k)) + Xi((p+1)*p/2+p+3,2)*sin(b(k))+ Xi((p+1)*p/2+p+4,2)*cos(a(k)) + Xi((p+1)*p/2+p+5,2)*cos(b(k));
+   end
+end
+
+%% Plot Results
+
+% Figure 1: Continuous-time solution x(t)
+figure(1)
+plot(tspan,xdat(1,:,1),'b','LineWidth',2)
+set(gca,'FontSize',16)
+xlabel('$t$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+ylabel('$x(t)$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+title('Solution of the ODE','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+
+% Figure 2: Continuous-time solution y(t)
+figure(2)
+plot(tspan,xdat(1,:,2),'r','LineWidth',2)
+set(gca,'FontSize',16)
+xlabel('$t$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+ylabel('$y(t)$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+title('Solution of the ODE','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+
+% Figure 3: Simulations of the discovered Poincaré map
+figure(3)
+plot(1:100,a(1:100),'b.','MarkerSize',10)
+set(gca,'FontSize',16)
+xlabel('$n$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+ylabel('$x_n$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+title('Iterates of the Discovered Poincaré Mapping','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+
+% Figure 4: Simulations of the discovered Poincaré map
+figure(4)
+plot(1:100,b(1:100),'r.','MarkerSize',10)
+set(gca,'FontSize',16)
+xlabel('$n$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+ylabel('$y_n$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+title('Iterates of the Discovered Poincaré Mapping','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+
+%% Driven Brusselator right-hand-side
+function dx = Brusselator(x,a,b,T,alpha)
+
+    dx = [a + alpha*sin(2*pi*x(3)/T) - b*x(1) + x(1)*x(1)*x(2) - x(1); b*x(1) - x(1)*x(1)*x(2); 1];
+
+end
+
+
+
+
+
+
+
+
+
+

Hopf_Section.m

@@ -0,0 +1,151 @@
+% ------------------------------------------------------------------
+% SINDy method for discovering mappings in Poincaré sections 
+% ------------------------------------------------------------------
+% Application to the Hopf normal form
+%
+%           x' = x - omega*y - x*(x^2 + y^2)
+%           y' = omega*x + y - y*(x^2 + y^2)     
+%
+% Here omega > 0 is a real-valued parameter.
+%
+% This code is associated with the paper 
+% "Poincaré maps for multiscale physics discovery and nonlinear Floquet
+% theory" by Jason J. Bramburger and J. Nathan Kutz (Physica D, 2020). 
+% This script is used to obtain the results in Section 3.2.
+% ------------------------------------------------------------------
+
+% Clean workspace
+clear all
+close all 
+clc
+
+format long
+
+%Model parameters
+T = 10*pi; % Period 
+
+%Generate Trajectories Hopf normal form
+m = 2; %Dimension of ODE
+n = m-1; %Dimension of Poincaré section
+dt = 0.005;
+tspan = (0:10000-1)*dt;
+options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,m));
+
+%Generate Trajectories
+x0(1,:)=[0.0001; 0];  % Initial condition (Guarantee one starts near origin)
+x0(2,:) = [0; 0]; %Trajectory at the ubstable equilibrium
+[~,xdat(1,:,:)]=ode45(@(t,x) Hopf(x,T),tspan,x0(1,:),options);
+[~,xdat(2,:,:)]=ode45(@(t,x) Hopf(x,T),tspan,x0(2,:),options);
+kfinal = 5; % Number of trajectories
+if kfinal >= 3
+    for k = 3:kfinal
+        x0(k,:) = [20*rand; 0]; %Initial conditions start in section 
+        [~,xdat(k,:,:)]=ode45(@(t,x) Hopf(x,T),tspan,x0(k,:),options);
+    end
+end
+
+%% Poincaré section data
+
+%Counting parameter
+count = 1;
+
+%Initialize
+Psec(1) = xdat(1,1,1);
+count = count + 1;
+
+%Create Poincare section data
+for i = 1:kfinal
+    for j = 1:length(xdat(1,:,1))-1 
+        if  (xdat(i,j,2) < 0) && (xdat(i,j+1,2) >= 0) 
+            Psec(count) = xdat(i,j+1,1); %nth iterate
+            PsecNext(count - 1) = xdat(i,j+1,1); %(n+1)st iterate
+            count = count + 1;
+        end
+    end
+    Psec = Psec(1:length(Psec)-1);
+    count = count - 1;
+end
+
+%% SINDy for Poincaré Sections
+
+% Access SINDy directory
+addpath Util
+
+% Create the recurrence data
+xt = Psec';
+xtnext = PsecNext';
+
+% pool Data  (i.e., build library of nonlinear time series)
+polyorder = 5; %polynomial order 
+usesine = 0; %use sine on (1) or off (0)
+
+Theta = poolData(xt,n,polyorder,usesine);
+
+% compute Sparse regression: sequential least squares
+lambda = 0.01;      % lambda is our sparsification knob.
+
+% apply iterative least squares/sparse regression
+Xi = sparsifyDynamicsAlt(Theta,xtnext,lambda,n);
+if n == 4
+[yout, newout] = poolDataLIST({'x','y','z','w'},Xi,n,polyorder,usesine);
+elseif n == 3
+[yout, newout] = poolDataLIST({'x','y','z'},Xi,n,polyorder,usesine);
+elseif n == 2
+ [yout, newout] = poolDataLIST({'x','y'},Xi,n,polyorder,usesine);
+elseif n == 1 
+  [yout, newout] = poolDataLIST({'x'},Xi,n,polyorder,usesine);
+end 
+
+fprintf('SINDy model: \n ')
+for k = 2:size(newout,2) 
+    SINDy_eq = newout{1,k}; 
+    SINDy_eq = [SINDy_eq  ' = '];
+    new = 1;
+   for j = 2:size(newout, 1)
+       if newout{j,k} ~= 0 
+           if new == 1 
+             SINDy_eq = [SINDy_eq  num2str(newout{j,k}) newout{j,1} ];  
+             new = 0;
+           else 
+             SINDy_eq = [SINDy_eq  ' + ' num2str(newout{j,k}) newout{j,1} ' '];
+           end 
+       end
+   end
+  fprintf(SINDy_eq)
+  fprintf('\n ')
+end 
+
+%% Simulate SINDy Map
+
+a = zeros(length(Psec),1); %SINDy map solution
+a(1) = 2*rand;
+
+for k = 1:length(Psec)-1
+   for j = 1:length(Xi)
+      a(k+1) = a(k+1) + Xi(j)*a(k).^(j-1);
+   end
+end
+
+%% Plot Results
+
+% Figure 1: Simulations of the discovered Poincaré map
+figure(1)
+plot(1:100,a(1:100),'b.','MarkerSize',10)
+set(gca,'FontSize',16)
+xlabel('$n$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+ylabel('$x_n$','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+title('Iterates of the Discovered Poincaré Mapping','Interpreter','latex','FontSize',20,'FontWeight','Bold')
+
+
+%% Hopf Right-hand-side
+function dx = Hopf(x,T)
+
+    dx = [x(1) - T*x(2) - x(1)*(x(1)^2 + x(2)^2); T*x(1) + x(2) - x(2)*(x(1)^2 + x(2)^2)];
+
+end
+
+
+
+
+
+