From 32e5eabfc3a2be01fce2badefbb521c8d0d8c3e3 Mon Sep 17 00:00:00 2001 From: Jason Bramburger <68243199+jbramburger@users.noreply.github.com> Date: Mon, 13 Jul 2020 15:08:47 -0400 Subject: [PATCH] Add files via upload --- Brusselator_Section.m | 197 ++++++++++++++++++++++++++++++++++++++++++ Hopf_Section.m | 151 ++++++++++++++++++++++++++++++++ Logistic_Section.m | 172 ++++++++++++++++++++++++++++++++++++ RC_Section.m | 161 ++++++++++++++++++++++++++++++++++ Rossler_Section.m | 172 ++++++++++++++++++++++++++++++++++++ Spiral_Section.m | 174 +++++++++++++++++++++++++++++++++++++ sparsifyDynamicsAlt.m | 36 ++++++++ spiral_data.mat | Bin 0 -> 12528 bytes 8 files changed, 1063 insertions(+) create mode 100644 Brusselator_Section.m create mode 100644 Hopf_Section.m create mode 100644 Logistic_Section.m create mode 100644 RC_Section.m create mode 100644 Rossler_Section.m create mode 100644 Spiral_Section.m create mode 100644 sparsifyDynamicsAlt.m create mode 100644 spiral_data.mat diff --git a/Brusselator_Section.m b/Brusselator_Section.m new file mode 100644 index 0000000..7cb9277 --- /dev/null +++ b/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 + + + + + + + + + + diff --git a/Hopf_Section.m b/Hopf_Section.m new file mode 100644 index 0000000..721ed17 --- /dev/null +++ b/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 + + + + + + diff --git a/Logistic_Section.m b/Logistic_Section.m new file mode 100644 index 0000000..08a7737 --- /dev/null +++ b/Logistic_Section.m @@ -0,0 +1,172 @@ +% ------------------------------------------------------------------ +% SINDy method for discovering mappings in Poincaré sections +% ------------------------------------------------------------------ +% Application to the singularly perturbed non-autonomous Logistic +% equation +% +% x' = eps*x*(1 + sin(T*t) - x) +% +% Here eps > 0 is taken to be small and T > is real-valued. +% +% 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.3. +% ------------------------------------------------------------------ + +% Clean workspace +clear all +close all +clc + +format long + +%Model parameters +eps = 0.1; % epsilon value +T = 2*pi; % period of forcing terms + +%Inializations for generating trajectories +m = 2; %Dimension of ODE +n = m-1; %Dimension of Poincaré section +dt = 0.005; +tspan = (0:20000-1)*dt; +options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,m)); + +%Generate Trajectories +kfinal = 5; +if kfinal >= 1 + for k = 1:kfinal + x0(k,:) = [4*rand; 0]; %Initial conditions start in section + [~,xdat(k,:,:)]=ode45(@(t,x) Logistic(x,eps,T),tspan,x0(k,:),options); + end +end + +%% Poincaré section data + +%Counting parameter +count = 1; + +%Known equilibrium entry +Psec(1) = 0; +PsecNext(1) = 0; + +%Create Poincare section data +for i = 1:kfinal + for j = 1:length(xdat(i,:,1))-1 + if (mod(xdat(i,j,2),1) >= 0.95 && mod(xdat(i,j+1,2),1) <= 0.05) %&& j >= length(xdat(i,:,1))/2) + temp(count) = xdat(i,j+1,1); %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(10000,1); %SINDy map solution +b = zeros(10000,1); +a(1) = 0.1; +b(1) = x0(1,1); + +for k = 1:9999 + for j = 1:length(Xi) + a(k+1) = a(k+1) + Xi(j)*a(k)^(j-1); + b(k+1) = b(k+1) + Xi(j)*b(k)^(j-1); + end +end + +%% Plot Results + +% Figure 1: Continuous-time solution +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: Simulations of the discovered Poincaré map +figure(2) +plot(1:100,a(1:100),'b.','MarkerSize',10) +hold on +plot(1:100,b(1:100),'r.','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') + +%% Logistic right-hand-side +function dx = Logistic(x,eps,T) + + dx = [eps*(x(1)*(1- x(1) + sin(T*x(2)))); 1]; + +end + + + + + + + + + + + + + + + + + diff --git a/RC_Section.m b/RC_Section.m new file mode 100644 index 0000000..63bf611 --- /dev/null +++ b/RC_Section.m @@ -0,0 +1,161 @@ +% ------------------------------------------------------------------ +% SINDy method for discovering mappings in Poincaré sections +% ------------------------------------------------------------------ +% Application to the non-autonomous RC circuit equation +% +% x' = A*sin(omega*t)-x, A,omega > 0 +% +% The Poincaré section can be computed explicitly and is linear. +% +% 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.1. +% ------------------------------------------------------------------ + +% Clean workspace +clear all +close all +clc + +format long + +%Model parameters +A = 1; +omega = 2*pi; + +%Generate Trajectories Hopf normal form +m = 2; %Dimension of ODE +n = m-1; %Dimension of Poincaré section +dt = 0.01; +tspan = (0:20000-1)*dt; +options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,m)); + +%Generate More Trajectories +x0(1,:) = [0; 0]; %At least one solution +[~,xdat(1,:,:)]=ode45(@(t,x) RC(x,A,omega),tspan,x0(1,:),options); +kfinal = 5; %Number of trajectories +if kfinal >= 2 + for k = 2:kfinal + x0(k,:) = [10*rand-5; 0]; %Initial conditions start in section + [~,xdat(k,:,:)]=ode45(@(t,x) RC(x,A,omega),tspan,x0(k,:),options); + end +end + +%% Poincaré section data + +%Counting parameter +count = 1; + +%Initialize +Psec(1) = xdat(1,1,1); +count = count + 1; + +%Create Poincaré section data +for i = 1:kfinal + for j = 1:length(xdat(i,:,1))-1 + if (j == 1) && (i > 1) %Trajectories start in the section + Psec(count) = xdat(i,j,1); + count = count + 1; + elseif (mod(xdat(i,j,2),2*pi/omega) >= 2*pi/omega-dt && mod(xdat(i,j+1,2),2*pi/omega) <= dt) + 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(1:length(Psec)-1)'; +xtnext = Psec(2:length(Psec))'; + +% 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.005; % 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 +b = zeros(length(Psec),1); +a(1) = x0(1,1,1); +b(1) = 3*rand; + +% Simulate section +for k = 1:length(Psec)-1 + for j = 1:length(Xi) + a(k+1) = a(k+1) + Xi(j)*a(k)^(j-1); + b(k+1) = b(k+1) + Xi(j)*b(k)^(j-1); + end +end + +%% Plot Solutions + +% Figure 1: Continuous-time solution +figure(1) +plot(tspan,xdat(1,:,1),'b') +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 with $x(0) = 0$','Interpreter','latex','FontSize',20,'FontWeight','Bold') + +% Figure 2: Simulations of the discovered Poincaré map +figure(2) +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') + +%% Rc Circuit right-hand-side + +function dx = RC(x,A,omega) + +dx = [A*sin(omega*x(2)) - x(1); 1]; + +end + + diff --git a/Rossler_Section.m b/Rossler_Section.m new file mode 100644 index 0000000..dc15b86 --- /dev/null +++ b/Rossler_Section.m @@ -0,0 +1,172 @@ +% ------------------------------------------------------------------ +% SINDy method for discovering mappings in Poincaré sections +% ------------------------------------------------------------------ +% Application to the Rossler system +% +% x' = -y - z +% y' = x + a*y +% z' = b + z*(x-c) +% +% Here a,b,c are real-valued parameters. Fixing a = b = 0.1 and increasing +% c from 0 causes a sequence of period-doubling bifurcations leading to +% chaos. Attractor at certain parameter values: +% c | Attractor +% ------------------- +% 6 | period 2 orbit +% 9 | 2-banded chaos +% 12.6 | period 6 orbit +% 13 | 3-banded chaos +% +% 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.5. +% ------------------------------------------------------------------ + +% Clean workspace +clear all +close all +clc + +format long + +%Model parameters +a = 0.1; +b = 0.1; +c = 6; + +%ODE generation parameters +m = 3; %Dimension of ODE +n = m-2; %Dimension of Poincaré section +dt = 0.005; +tspan = (0:1000000-1)*dt; +options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,m)); + +%Generate Trajectories +x0(1,:) = [0; -8; 0]; +[~,xdat(1,:,:)]=ode45(@(t,x) Rossler(x,a,b,c),tspan,x0(1,:),options); +kfinal = 1; %Number of trajectories +if kfinal >= 2 + for k = 2:kfinal + x0(k,:) = [0; -4*rand-8; 0]; + [~,xdat(k,:,:)]=ode45(@(t,x) Rossler(x,a,b,c),tspan,x0(k,:),options); + end +end + +%% Poincare section data + +%Counting parameter +count = 1; + +%Initialize +Psec = []; +PsecNext = []; + +%Create Poincare section data +for i = 1:kfinal + temp = []; + for j = 1:length(xdat(i,:,1))-1 + if (xdat(i,j,1) < 0 && xdat(i,j+1,1) >= 0) + temp(count,:) = xdat(i,j+1,2); %nth iterate + count = count + 1; + end + end + Psec = [Psec; temp(1:length(temp)-1,:)]; + PsecNext = [PsecNext; temp(2:length(temp),:)]; + count = 1; +end + +%% SINDy for Poincaré Sections + +% Access SINDy directory +addpath Util + +% Create the recurrence data +xt = Psec(:,1); +xtnext = PsecNext(:,1); + +% 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(500,1); %SINDy map solution +a(1) = Psec(1,1); + +for k = 1:499 + for j = 1:length(Xi) + a(k+1) = a(k+1) + Xi(j)*a(k)^(j-1); + end + + if usesine == 1 + a(k+1) = a(k+1) + Xi(j+1)*sin(a(k)) + Xi(j+2)*cos(a(k)); + end +end + +%% Plot Results + +% Figure 1: Simulations of the discovered Poincaré map +figure(1) +plot(1:100,a(1:100),'b.','MarkerSize',10) +hold on +plot(1:20,Psec(1:20,1),'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') +legend({'SINDy Mapping','Training Data'}, 'Interpreter','latex','FontSize',16,'Location','best') + +%% Rossler right-hand-side + +function dx = Rossler(x,a,b,c) + + dx = [-x(2) - x(3); x(1) + a*x(2); b + x(3)*(x(1) - c)]; + +end + + + + + + + + + + diff --git a/Spiral_Section.m b/Spiral_Section.m new file mode 100644 index 0000000..398a9aa --- /dev/null +++ b/Spiral_Section.m @@ -0,0 +1,174 @@ +% ------------------------------------------------------------------ +% SINDy method for discovering mappings in Poincaré sections +% ------------------------------------------------------------------ +% Application to coarse-grained dynamical evolution of spiral +% wave solutions to the lambda-omega PDE +% +% u_t = D?u + u*(1 - u^2 + v^2) - beta*v*(u^2 + v^2) +% v_t = D?v + v*(1 - u^2 + v^2) + beta*u*(u^2 + v^2) +% +% Here D,beta > 0 are real-valued parameters. Data is gathered using +% D = 0.1 and beta = 1 on a spatial domain with x,y in [-10,10] with +% periodic boundary conditions. Simulations are performed using spectral +% methods. +% +% Solutions are simulated and projected onto their two dominant PCA modes. +% Temporal data for the PCA modes is contained in spiral_data.mat. +% +% 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.6. +% ------------------------------------------------------------------ + +% Clean workspace +clear all +close all +clc + +format long + +%Initializations +n = 2; %Using only the two principal components from u + +%% Aggregate Data + +load spiral_pca_series.mat + +%Create section data +for j = 1:80 + xt(j,:) = [pcaSeries_u(5*(j-1)+1,:)]; %pcaSeries_v(10*(j-1)+1,:)]; + xtnext(j,:) = [pcaSeries_u(5*j+1,:)]; %pcaSeries_v(10*j+1,:)]; +end + +%% SINDy for Coarse-Grained Forecasting + +% Access SINDy directory +addpath Util + +% 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 = 10; % 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 Poincare Map + +a = zeros(1000,1); %SINDy map solution +b = zeros(1000,1); +a(1) = xt(1,1); +b(1) = xt(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 + +%% Uncertainty Quantification + +sample = 5000; + +A(:,1) = normrnd(a(1),0.1,sample,1); +B(:,1) = normrnd(b(1),0.1,sample,1); + +for k = 1:999 + + % Constant terms + A(:,k+1) = Xi(1,1)*ones(sample,1); + B(:,k+1) = Xi(1,2)*ones(sample,1); + + %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 +end + +%Eliminate blow-up terms +A(any(isnan(A), 2), :) = []; +B(any(isnan(B), 2), :) = []; + +% Figure 1: Simulations of the discovered mapping +figure(1) +plot(-A(:,1:50)','k.') +hold on +plot(-a(1:50),'b.','MarkerSize',40) +set(gca,'FontSize',16) +xlabel('$n$','Interpreter','latex','FontSize',20,'FontWeight','Bold') +ylabel('$\tilde{x}(n)$','Interpreter','latex','FontSize',20,'FontWeight','Bold') +title('Iterates of the Discovered Mapping: Component 1','Interpreter','latex','FontSize',20,'FontWeight','Bold') + +% Figure 1: Simulations of the discovered mapping +figure(2) +plot(B(:,1:50)','k.') +hold on +plot(b(1:50),'r.','MarkerSize',40) +set(gca,'FontSize',16) +xlabel('$n$','Interpreter','latex','FontSize',20,'FontWeight','Bold') +ylabel('$\tilde{y}(n)$','Interpreter','latex','FontSize',20,'FontWeight','Bold') +title('Iterates of the Discovered Mapping: Component 2','Interpreter','latex','FontSize',20,'FontWeight','Bold') + + + + + + + + + + + + diff --git a/sparsifyDynamicsAlt.m b/sparsifyDynamicsAlt.m new file mode 100644 index 0000000..2d83821 --- /dev/null +++ b/sparsifyDynamicsAlt.m @@ -0,0 +1,36 @@ +function Xi = sparsifyDynamics2(Theta,dXdt,lambda,n) +% Copyright 2015, All Rights Reserved +% Code by Steven L. Brunton +% For Paper, "Discovering Governing Equations from Data: +% Sparse Identification of Nonlinear Dynamical Systems" +% by S. L. Brunton, J. L. Proctor, and J. N. Kutz + +%SVD +[U,S,V]=svd(Theta,0); +Xi = V*((U'*dXdt)./diag(S)); + +%Xi = pinv(Theta)*dXdt; +k = 1; +Xi_new = Xi; + +% lambda is our sparsification knob. +while sum(sum(abs(Xi - Xi_new))) > 0 || k == 1 + Xi = Xi_new; + smallinds = (abs(Xi)m@O=Y15H}}>aOqV zv)*IHE5CbY{jApQ>HYi)`{3Z9C28+KZ%9qF@3Z0aN$0KK=MR)nU%0>#LLrM2Bu~KP z;u3s%c|{2MTtT>rfN>Jm5U>IFGQwKl2zx}|<9^iZ#`&ju5e9@!T#ppJ7lu6@K?-XF zwYwcD3xw5D5AwuBP2G<=UE#0}AZ!@IpH)MGBs03d=$4I1J+(Tn>Lfd|!cdh2@o!;? zzi)kUUpmjda#$Y{$q(8bEf34qPTS%{=NORa9(lM3DPxE$gk)Zs4*pr-4RXmRvG>2? znD-|y2%QBWe$#zMIDf@tcB~aIdP<&e5^#cowhZD=2V+=UGhMnv_<%@+#2)4D7}}r; zc4I-YUC&9?^gIM8am+Ep?4BQuuHvVDAWZ1Ve|lV@>;j;?H@ex!Em)g>o#LSEMqurs z%aNa`YY0Wu%6Vuw%nWYPC;YO~qWWlh$Z;$AwTBapz1p<}qu=>5{(PK-`!;teOIs0E zG}vyqZw=9505tRZ-bQ0b7~o*^S)HN&5jv;$k;M_U_0c~_jZ|q&ZFTGB_rmj+(s%~Z zR}78KlhC;!W)}%JB$Zp%I4s}@SZCSsmR$Y<$sD&Ckj%l9R}0c$h@;+xnbHyGYIw%zCFHE~q< z?ZS_DmDn!{hGbz&x8 z%X@pZkq@;bmLId`;Hor-P9Ji~!k!8qcdQCOBcT8WU&B9y zX%W^6dA9w_R)+C5l<(}eiDiCcU3Yr-&8Ezf5Wd$AgunzM(XQ8~;YSgCb3_){!U`|< zWIbIN&1?$0=$$-SVM$;j8Lfv+yaIi!@SWniZSksx8?w#bS@nrhky|I>@-mhb`d!1!$){D_{ZF-Nnl_-eMHUmub zNSyi&;#!71b_D2NGgQ7sqs;!Wg4LJy?`PksHVhwe(%Gm3|3o+k+ahir0y>%8>_l>y zfZZTSJGjl|T5+ye0l@h+)ih$o0exOhp_Z_XM5|nF8CoHa51-ofkV=tM1={@9DyW*?VT;qU<)o_Rq$?I> zXwGmf^=cZFxFyQjL$rS--c$%}vp z8aPAdZ7X}wByf?d4mRfx3}%1tz#IdWF$KM|+7dsMDsrTlkRwI0L0#?RUl~Ac39L9> zv$S1XpwW591@HH~R&vk?2z7`6OpVf0M{mHwpE)HvZ zs}$0dZ{8BUvj-^;^CKIGfT@o)us3Ydz!Xyv(+U|Gska388S__S@{Cn4PIS`G`Ar4CPjGt?iu{=$r{AHrPG;wa$O2K}JGIJwlK5a;0p0rNg&iJ}l(M+Ow8B z#nNGonmzxa+;6v-mR5s6U6ysmXOQ?HewOd=-BWEjtR*7<524LWL{SSrKh+Gi-Gr@* zH6BG+pw|Pi{vME}OsKeTjtO-%bWFh6(q;YcT;Q0A|Qv29PiUneYm=iC!Zq>zvL86jrV81#$NRYFUBX@px(|SXuq_s+} z9fYmrT&Ji+}G%C!AVT?ea5S~ zW^kvk^)*q=?#6QWZ}2)N{s#hYuW=9>k)URGWE2YfwQcosNL{GBeG)HIb4h4;4X;-j z8tWR5cX|VeDJ;~FQvuCy0)SgJy52KiI8f?}`00L$UfD-==duIb%h2hw(K3K|<_D0^OU4KqiF#P7 zuZ;T1l?Ra+pfDSpAD|;!;vv1=cK-z+d`VPqR3*3n@*P8wTsWHO;l)Etx=5!b@YaI_ zN(&TM-9s;4^x>Gk8utYiai`6mv(mJ5=TKU9G5~vLR`6Hm>_K_kvl~jN(6{QSgoWwN z)030l@`d1{Vt2i@GQCVY6)M(yV7LZKz1>rv2Bk`T_pI;4Z%|o@rie~y0uC$~v@)Rt zuepC59#(EI!&~55XA82>iNOarY-EhZB(Zx+O0TeyfCEM>+i(7f310Ni#ouN_oZZ`R zMj|o8R!~}tM<&cgdgb**R9hY-5fIMBV9ik?oNL*5#aNokUMNjcnKL11Zwdzi* z5)&1acD;vJ2D(cT>sKC{Y~({AAv@dyGJK$n9=TI0%KFVb-{j2Up!xD7f?)zMcT^;& z(+fen9%;6>vyG4|GVk@Zwh?OM3bKOySqCFn*7X2AR%uMae??}X2dkK{N-ZF1SG zWs>E}7*#vH?8R0JOfj|tH{uX+Sn5pnis4;vk1fq}$07U<(kl1Pw-Q~T?Sf=t+rvBH zH*%{U55%1+m+$DIp6C0B+hzfZ-hb?lt&l^NM+q;bmcqO_Fu+ro{d%G~yff9=5z}AEg4Jff&l(Q3MSv*2iA+bf zoQV0%sf@bNzd9mt!5F?nAw>T`3Ptc)QPEodo1YTi9;sL1j^o%kY9(azNvYXIj!J=Y zvoL%!K&c|m$b4V%CL+_xAD0Ki{1wuNLF5LoaoX^($qH;L7IpGqWSZ`xz+qR351p5U z#h9Y&A?>riErY{)KWin}@D6;93rxO@=*BvsJFNAduOs|1X$ykPW*Vnx9!(gA2J5>jQt>^j0f+KpfrJGPT&b78M|dD?`i>sP9#zbxD#{%a6+-v26L7HuBp0g55r z`HZMwKJngdV`ja1=5G(SOp%Usc*Xmxf#0&=kAsX@caY#}%Zk2yLgNT7_yW~D+zz+6 zNLpo=dzOfkZG2*QOY6r1cgxQvwr_L{E$*TIr?M1Y9t=K9NG@x9H{CC5dz#n-WqzH0 zJ)_?Pb&CD7=72R{$YQ<{DiyPsjUb9|UM^Dk`!A!l0l0gT-<6wcH@JS*P|BY9t+t)a zDH@gUJnrMV{&(WnE(MDcyk7Z6m89xbJ2N2Iux#t|#D0!jK=*?_4jSBHTY23;xw{e* zH@i0o6ZO{%l|z91ixUx#e>AwI-+?WAS39*x-p50ZEE*<_uF=ehTGT~1K;9}V>}gYj zOTX~P9+)&26968X|7a78VbeDrS;Jr#2;e2Q`raRh1o{z<$p4uLJ^E*di}o zumtqzQ&|RQ&Zg#WgM?fu_Wp>JE~@YXMliIKkqR-CQ<{{X-H^YU&CDx3IguMbyaDgZ zPYIGk%UGINpL(AhYjsLpNHgF$unlfyL)p%7|@7n6bo3@RJ z)pb(^1HsYoI#t9J18yf@#?Ri0@SixB>@M5G8_NnJ#?OBF{_z+4?$(?JQjF!N$V*({ z9T&k{Efu1MtMS+pMdH4V6qE|KL3wnIm3P{49!%wYaSCmLQrVOwqnn{bJZsZ&^D$yN zh^r}DviMZ@_ifgoW+BDMZCg;R*0{S>E8$)4WG7L|sG)Aiq8UmqD&uPrmP;-_DeLuV zCnP|XeHS=4h*ls{!9aw@>@j!>xNK9za(oDEJCJnbTnOAwW5cTv9e*rKp`GC1utbVe zBq}yL<$Lc-X89(Y46|g#eq%HNwNeAX=7K-VYCCG0KaUx`Jf%VXH)5yzz9nl6-XND% zqyjdHb2w~RL0#;4w!Q5w(<<>gbf}F}Fgfl#z5b+*_MrMSf5CQZGBTb2!-m4dRsC3) z+kRk5J}mxWLW(2lMCpT|zeg@#FH@@s-iNgaJv*||k8_x3&5JoQmHOo&TZe{&Meh-~ z-L_wHa9pgjIjtXuu?V<+RWRLL>cHeQ%t zK^_P0W@OLPVy2!066jqzObqTD**!22twTWMGt{{U4Ei2>)q6(lR6?3+nJ1vZfmK$v zVNtmg$mpA?ma-LNJ-?Ti-id`_VTnVh*+IRnejQ&{drM}Fpwqo4p^%cOZSK;K_%6h) z_he$RM4B=uAm8cFcsnMNYI)&ZqGA z?W*;))mF~!tU_?)g(+L2R=IeRhc=WgzHD~aUY?l zjD$m{Wa=}$x`>`-bX@-8#le)akUC52O6n++zqq0Xz?kCH5nH4x<8l+?d&G>@4Th>B zW!}Afw!8|eIiV}NuzX}(s6f39Wj@v(J|IwvCBq69uFowba zF8hoJP-D(trC^9sEES1-@PCGK9Yt5{z1o^s$Nt0US6aR{bvj$e(9v+NNmqaANJqRi zD(FPqst@bEuzv`6!D%${FtjNhm(Wj)XS*ELP2UT)OVrZTYp`X+|_`(FuFtU={8(d2V#fird%8n z&)1I`YOBVx#yPrlK1X_OGuD`;NCphw%$a@5qf%a35q@q*h!%P=QhYv4I(e0CFu8|~ zwkcLz5FY9D4oYu#uI0kHs_*Jy+DeZldonF233qcD$4_*a5lHyO_?AV$NJ$++1lOscU*&ac(WRYbYgUD z&B(B(L3MAl76h>Af&F?`}c$2d@CQkN?gUfJ~7kkOGthL}kDa2oV*{OBL{rY#7zgxEi7svLMeJ0CH z=k-{AS(IH+RQ zU0Fhp)c&U~?YhY%WHF=9sW%5Xs_s|AJGodSwBJxmUE%Uj(dnqNuPw?55_+kH?s8JO z`3Z}VM7T-N;Qlrajj}ykiE&>gK)>Un>+VM%VPEBES)aZ-;xPG31tm2S@UjeIE+t3Y zv$AEzlaYeb?sdiWm}3!qIQ5-&{?wDJ1P_NY%-!AI@C0-o4jTMc z`i-B)7qg!NxbLusS|l+cSn0n;ll$|w9^mY+}N#neF| zE{#e|$72!9hI;zJsA$;ZPJJqI4?pN($~Pyj3A@ymugqclgjwsGN()!D95I_9#Y&@i z9igTD)z!+aLZ&zAUAA#AOK-=efsB2DJ(dd6m?#g>LNS!+0`)w~s1fe$&K}(DtBI^l z^X;?AtTSRQe>!&x0=af%(Qhqw_KSWN1ct?a2clHKm37;}zX^Aa9PR(#-2m&8;~1}3eXY`2)~hH})o1@&v+_L(crC7YWX?QJoI3LS!s#G_|A<$l-f z-5WTM_lh0tZ^i3oX?V42Rh*QRRF{`Dr!N|V^#1wQpS$*%kQx>qCILeWUkn?vJv=(T zvwW6!mYyo~>69O4Y(&8$EO0jtve4C0e(y(47OXL?%0OBcC2``G#uT|%{Ac~l0z32m z&Mj!63Qvo2R4;SG(dbRHadO=Q|JQ)XZ)FN|c9+;2Jc9m<+%07j2P>_mGB%1QF%7HqM=pvLFc1JDuXp8cI4zqQ1!WngTXpNJmxdZ zpKgM1LbEZ%%nTD@h19XNi6;1vkIrt2J4rTWhb1+iB^c$im`YZ~*odM1QUtYC#{t5g z6I$6DxZ^^zb@m+`kxab&K;~uigxv!?O;j9Q11!xOMm$XT%=~Tbw7@ ziGQzT{bS%DcRq~o9qaLiIsn6SQ?;V=FwhgB_k>AkDE!#o%%CV>;M}1*mD~Rnm0v1_ z8m%B75E{`^Hb^wnp9Kmc{a@voO_!3l+jh(%H)iz34-A=WVCP+NZ^6Wo`)YCFCO{GbD)RYq=tf zn+Y#HgguyGX(E+x$s+QtJB}u^`^4%(THnk6%iF~Wmrr`>ZeeQO%q<_J?Zj1#gGeq^ z^6}?6xk%|&fka^>v@B&@b)9qExy#a6=P`jpyWpnEw!wj_Hann+?#LL`|BaGxzK@%P z_+qwlD_;Sq*WRVqJc4Ml?omt>Ze?TK<=6imDONGSdz1<_^P7421vFc|NDlfByX`*_qAkp7+e*AeQojfbOqVB?};uVKLe#P|e4i%HwqU4*6 zer;Cbe6K8s_Aj1^`Q?fP$M#C~z=+mu&_2s2XKMl|{=~rlh}&%Pp&LykyHVEY_x7dN zUjt#u#gR4JOy1Z@BQ2)dB3RN{=USVG2)mAu<1ZIYL@?fFw=6UVDS^Ji!wEyZ>m7I5Ntx9ZMo#EOjzZj zCK?)R$=+iV4ED=zZI2G{T z)Hr_sr0{{l(Wrs1sT?y{y77?+7A3Ejj$5bpr-H>$Emdp{D?3O#bKSd4V^9yg-gWo> z$@G#9IA?T$T{G#dHF6`sGOFJ!wiRI%S znJQ+QzLNJi0Y{Orp3hIPq`2R?u!7cWB%SshbhIp}Xlf$-8o%9iA$pM?L=k})S>NO(1o*iW7p?5>|hr}4j^o15r{sc-N< zya96nJLEJCR}DmVM}IaN{Vt5u!LxN;0qS$)dFbP&kBcVG?JsV|uF)(%PBs*(#_U{9 zE}8>molS=R`JYWPq;>eO)H+N^2!`4ertTxZNzkG`_^8JvX1;S8o#}*4WvIE*XU;nh zcPp47`KbHzpLQVvLo96=UfaV2D0t7HCjc!swu7z6XFtE?DE{PF>#Qg-G_b2Vs;9V* z#x}rgq(oJ6PV@7g9dZFHGSfaUkj9x&ti;_2sxYl6O23^3V%n{x?uqcj7UysuzpVR) z7DW3zBsB-q*U10qP5DvGuymh3t1Zy8>s5Yfkcs?{JFAWGIO#OiK6xak|GEoppb3+KwjVG-Lv^7|Z=cZJq0C z@B3FRh&Tx+0H#lM!LRLv80j87MWkrY1pK61l_0)g3a5`vVCsAUqEJhy^%42b03()A zJ=JM(_yR|!^XKQC9l@h~x2ad+qPy^k znFB=_D3-O8|NJW!mD~3!V#FFuU`&HwM{ExY0IQ&VEJ5G=eD*dzMTVeO&?g07c z0{wnCOIhYW^*lKs-pV*C`=R+O z`;9`Jx;c6~!-pliOvC|TG4x-(iKaLNN*?}!OcDeqV zb;#*^)KSNzPp6va@2ga7SW5N@Nv}iYrhK((3n9HI`l)t#Tq0Vr{IKm6!zoh+W340G)ay+6#*BNR) zMm_c4FaoEijlv@*+4MaGUVV6;l;zsOv4CAjQ|)4?K;Fb|9KQ#EG-F3fy}rV}{fES# zd(p3Ons^?JXg_XTv_S*f-=jBT^cH^9NS{C3di6I)kEFazt9I<9tXTrfF`$f8ZlVD% zTGZD#q16Rqs=Mr5I)*_hu&c4^I?SF0BN=FpilHwQfv=rQ`1*I%~fy>my_c z3F8N)cd;6!nG_OIBqo^gWzcM2Ke>T8j`W+Nd0WjiXYoOyj`&gY4Y|PFuYhNY_z14R zJ|lWh>=NBv11D{2SU3{Cv!{+@v%2ra^q%p-iJy}tt5QP84CJJ4=3Wz1j>IFxgb?ei zq401OT_Y`nz2NG-(M88T`#BPf757^j*g}<@tYkBU&Hx5?kjYR;P^yF3FG39B$~TmE zK501iPY)Nab5q=a;eeI{`|sY13Rawl){P-0_BhW`h}`82O0AvBq@Ob}MJQ@M%GkAO z<2w1!BmaIRzdo>=Jpa0Nnl!MPYlgh`Q^lIE5xzRFyc{!1e|k~bt*+B_x*S4M>wcP` zuk8CYYaZjztxhTxF7hIg8^vEF^e&A28E&(~+WCF3F=xc&)(x=UiA2)1%YEoWfy%t}Ay9o`I*eo3zt=a1_oqmlEl~d* zY+~O#m3B}2z^I9&fv0%=GekoW=^SNK z;A_qUOyHWA-`Sf4@s9al|JH`x{mx-CMR$5QG-@8nzfa97TJ6&n(C2kG8lM_e2N^40mNTK(m*pQBi z!$?8zNM=wT>7Zbb*$2|C10F6LwBM>|qlXOSBCBV8j$5+}V{xO^;?p844dJy|YAOa2 zT6WL-UD8VJ);&9Lg1Eo^kP*ZqqT37ar+KVp9y?Di2LJl5=7!iz*hhFK{}hp3?6I>k zQ{1TKtI;t&D}BACR#|x(fYiK05Pxy67=hBdHM@Av;oQ9epj}Ku2*NzHn`(EM*gL^hyEMM9$jN!Ui*h$w(7 zx|f1rRVOz5E==YxlCkLGV?!BY!JXB1BMys2r^$0ebG+%7vnD_Oj;9 zVbvnhgm;4W-uQ>$&J^xOa(kW|J4U^A$z#w!8~Zh<(f*JgPe0v2C=w>Qb1s`%y!y7FAZADgCw`@pm&@fLy;Opb7pMZNY z@C3gXy7P!<|HqUE)SBkae2W}{GusT|QrkYz;O0*)VzJAGA3I?Guu_3b`C1ZRcdJ*kuMMK4FFYTA#Os&o}uS z8rW95^}TFJtzOVc2_1H;a73n-U>f=&0>8}2VV&hfCp~>^GPtCQxcJ+r`dK! zkV#^2V>L6jeH|A!{}nP^F72dHa=!-{}M|mijqcCX2gIx2$AG#2p!o0>^ zzt^ImYT6JKQJm7_kfmzL(KpWBT2?1L@_(`f48Ohy7N9BbO?`Jb**swJd{zR=WaQ zh2@#;C)FN=&`==>Gxqw~i%TJJIceYRt@o?uDV>~hF&yujl~i%~!QxbW^A9M4Z9LpmK%|u@8L113h`bw98ltr(@k0K``j|Fo-EU-K!7c{h)~dO8xHqAakKT z`)!%{jif=NGpk>CKnf#tr6WSf-#wKe^@V75EQdcPHxPWCEYu+^o;f77eD0><)%+}N zz+tmXFgX^NF#cPaS!iFbO`(^|KaIDj7H&*LbnCojS-C3D-OFvT zO!!c@&4d_6GMSN#de}?%IkV56sKk*P?L041U%_syU<(K5+j@pp9QyvIKpQ^GeC^|= zCQdxR3s*|nLkyO)`R22!U>d@N`+&t;Hq1yPepw!ySL58+jSinaMVeBQx63U;zd1dL zS>}yP*!ZQ;$$zpTdxIyYm$%pO(zR~hHDfsX!&rs{jMPLtfczKpePXz6OY%3>-6fph_^0+|8r_#XLYJ;aV_ zeJ0`bKQcDpTPs&b&*Pj#94F^MRHgB zi557bzEW>jI*jQOfGnJzDa;wDWaRGEh7ltnQ^<8oPzD^A)=*OM> zhMr~j{Nlbuvc+i2=t<9z+U8)gE;0=`r4>DK6&eP*V2RnF2qoAzO6>I(RPr)DNB zx%npD_tBdQMPzrEsm}ZXQ`0?ov#t8Xa8?l7NcVW!en=l^;dky@XXUN*jjb6fYH#y) zrrp%&r-uqmd3jU=JKp3F4lEt--VJ2-}{i4A0)uOT5=L2zAataA2d$G-T8pYrbcBey5I zOn*NNj^8Yj`a%-Uu7&+G|G@0;-!p2X$J_azQiVgBM`m614z}Pd<9gv7viy%E2}=vt zfllG*@1Q0}WY`fcanc!HPGaKOFn&0Y2kU$r^HvHU*;#15AEw|47ypdPY59*N`c$)tpZ zMNoOk?--W-DX|$nCww76qZI*i)hHnBB*OtY4Ncx?6;5tyVh~rHoFl=^&*tJuk06C_ zD>Dt@>z(R^-hq6NuYD98NeKraL-fdfD$jsYoT2%Dcl&JHQ6lbhNtpsw`@nxst_u6a zxuwlpbf?+u_1ESoqhl(UIYqNvN^;w1&`$!B4X9%!CeNeeJ9= literal 0 HcmV?d00001