% Manifolds for Duffing's oscillator. 
% Comparison with FTLE and FSLE.

% Francisco J. Beron-Vera, 2006/15/05

clear

%% Parameters, etc.
% Equations.
vep=.1;  
T=2*pi/3;
% Manifold computation.
v=linspace(0,1e-3,1e2);
ds=1e-3;
Nt=4;
% Color scales.
BwdOrFwd=[1.00 1.00 1.00
   0.80 0.89 0.90
   0.49 0.75 0.75
   0.65 0.84 0.52
   1.00 1.00 0.50
   1.00 0.50 0.25
   1.00 0.00 0.00
   0.78 0.00 0.00];
BwdAndFwd=[0.00  0.00  1.00
   0.00  0.70  1.00
   0.70  1.00  1.00
   1.00  1.00  1.00
   1.00  1.00  0.70
   1.00  0.70  0.00
   1.00  0.00  0.00];
% ODE tolerances.
options=odeset('RelTol',1e-3,'AbsTol',1e-3);

%% UNPERTURBED MANIFOLDS.
figure(1), clf

% Levels.
H=@(x,y) y.^2/2-x.^2/2+x.^4/4;
[x y]=meshgrid(-1.5:.05:1.5,-1:.05:1);
contour(x(1,:),y(:,1),H(x,y),[-.2 -.1 .1 .2],'k'), hold on
contour(x(1,:),y(:,1),H(x,y),[0 0],'r')
arrow([0 0;0 0;1 -1;-1 1]/4,[1 1;-1 -1;0 0;0 0]/4,10)
text([1 1]/3.5,[1 -1]/3.5,{'$W_u$';'$W_s$'},'Interpreter','latex','FontSize',14)

% Attacting/repelling nature of Wu/Ws.
m=5;
theta0=linspace(0,2.1*pi,10);
r0=.05;
[xc yc]=ginput(2);

x0=xc(1)+r0*cos(theta0);
y0=yc(1)+r0*sin(theta0);
xy01=[x0;y0];
[x y]=manifold(x0',y0',ds,'spline');
h=fill(x,y,'g'); set(h,'EdgeColor','none')
x0=xc(2)+r0*cos(theta0);
y0=yc(2)+r0*sin(theta0);
xy02=[x0;y0];
[x y]=manifold(x0',y0',ds,'spline');
h=fill(x,y,'b'); set(h,'EdgeColor','none')
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('$\dot{x} = y$, $\dot{y} = x-x^3$','Interpreter','latex','FontSize',14)
drawnow, pause

for n=2:m
   [t xy]=ode23(@duffeq,(n-2)*T/m:T/m:n*T/m,xy01(:),...
      options,...
      'forward',0,T);
   x=xy(end,1:2:end-1);
   y=xy(end,2:2:end);
   [x y]=manifold(x(:),y(:),ds/(n-1),'spline');
   h=fill(x,y,'g'); set(h,'EdgeColor','none')
   xy01=[x';y'];
   [t xy]=ode23(@duffeq,(n-2)*T/m:T/m:n*T/m,xy02(:),...
      options,...
      'forward',0,T);
   x=xy(end,1:2:end-1);
   y=xy(end,2:2:end);
   [x y]=manifold(x(:),y(:),ds/(n-1),'spline');
   h=fill(x,y,'b'); set(h,'EdgeColor','none'), drawnow
   xy02=[x';y'];
end, hold off

%% UNPERTURBED MANIFOLDS (FTLE).
N=250;
x0=linspace(-1.5,1.5,N); dx0=min(diff(x0));
y0=linspace(-1,1,N); dy0=min(diff(y0));
[x0 y0]=meshgrid(x0,y0);
xy0=[x0(:)';y0(:)'];

% Attracting LCS's.
[t xy]=ode23(@duffeq,[0 T Nt*T],xy0(:),...
   options,...
   'backward',0,T);
x=xy(end,1:2:end-1); x=reshape(x,[N N]);
y=xy(end,2:2:end); y=reshape(y,[N N]);
dxdx0bwd=prmy(difx(x)/dx0); 
dxdy0bwd=prmx(dify(x)/dy0);
dydx0bwd=prmy(difx(y)/dx0); 
dydy0bwd=prmx(dify(y)/dy0); 

% Repelling LCS's.
[t xy]=ode23(@duffeq,[0 T Nt*T],xy0(:),...
   options,...
   'forward',0,T);
x=xy(end,1:2:end-1); x=reshape(x,[N N]);
y=xy(end,2:2:end); y=reshape(y,[N N]);
dxdx0fwd=prmy(difx(x)/dx0); 
dxdy0fwd=prmx(dify(x)/dy0);
dydx0fwd=prmy(difx(y)/dx0); 
dydy0fwd=prmx(dify(y)/dy0); 

FTLEbwd=zeros(N-1);
FTLEfwd=zeros(N-1);
for I=1:(N-1)^2
   dxydxy0=[dxdx0bwd(I) dxdy0bwd(I)
            dydx0bwd(I) dydy0bwd(I)];
   FTLEbwd(I)=norm(dxydxy0,2);
   dxydxy0=[dxdx0fwd(I) dxdy0fwd(I)
            dydx0fwd(I) dydy0fwd(I)];
   FTLEfwd(I)=norm(dxydxy0,2);
end
FTLEbwd=log(abs(FTLEbwd))/Nt/T;
FTLEfwd=log(abs(FTLEfwd))/Nt/T; 

A=FTLEbwd; A(A<2.5/3*max(max(A)))=NaN;
R=FTLEfwd; R(R<2.5/3*max(max(R)))=NaN;

x0=min(min(x0))+dx0/2:dx0:max(max(x0))-dx0/2;
y0=min(min(y0))+dy0/2:dy0:max(max(y0))-dy0/2;

figure(2), clf
pcolor(x0,y0,A)
colormap(BwdOrFwd), shading flat
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (Backward-Time FTLEs)','Interpreter','latex','FontSize',14)

figure(3), clf
pcolor(x0,y0,R)
colormap(BwdOrFwd), shading flat
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('aLCSs (Forward-Time FTLEs)','Interpreter','latex','FontSize',14)

figure(4), clf
pcolor(x0,y0,A), hold on
pcolor(x0,y0,-R), hold off
colormap(BwdAndFwd), shading flat
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (bleu) \& aLCSs (red)','Interpreter','latex','FontSize',14)

%% UNPERTURBED MANIFOLDS (FSLE).
r=100;
tspan=linspace(0,75*T,250);
N=200;
x0=linspace(-1.5,1.5,N); d0=min(diff(x0));
y0=linspace(-1,1,N); 
[x0 y0]=meshgrid(x0,y0);
xy0=[x0(:)';y0(:)'];

% Attracting LCS's.
[tbwd xy]=ode45(@duffeq,tspan,xy0(:),...
   options,...
   'backward',0,T);
xbwd=xy(:,1:2:end-1); xbwd=reshape(xbwd,[length(tbwd) N N]);
ybwd=xy(:,2:2:end); ybwd=reshape(ybwd,[length(tbwd) N N]);

% Repelling LCS's.
[tfwd xy]=ode45(@duffeq,tspan,xy0(:),...
   options,...
   'forward',0,T);
xfwd=xy(:,1:2:end-1); xfwd=reshape(xfwd,[length(tfwd) N N]);
yfwd=xy(:,2:2:end); yfwd=reshape(yfwd,[length(tfwd) N N]);

taubwd=repmat(0,[N-2 N-2]);
taufwd=repmat(0,[N-2 N-2]);
for I=1:N-2
   for J=1:N-2
      d1=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I,J+1)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I,J+1)).^2);
      d2=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I+1,J)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I+1,J)).^2);
      d3=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I+2,J+1)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I+2,J+1)).^2);
      d4=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I+1,J+2)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I+1,J+2)).^2);
      K1=min(find(d1>r*d0));
      K2=min(find(d2>r*d0));
      K3=min(find(d3>r*d0));
      K4=min(find(d4>r*d0));
      K=[K1 K2 K3 K4];
      if ~isempty(K)
         taubwd(I,J)=sum(tbwd(K))/length(K);
      end
      d1=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I,J+1)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I,J+1)).^2);
      d2=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I+1,J)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I+1,J)).^2);
      d3=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I+2,J+1)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I+2,J+1)).^2);
      d4=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I+1,J+2)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I+1,J+2)).^2);
      K1=min(find(d1>r*d0));
      K2=min(find(d2>r*d0));
      K3=min(find(d3>r*d0));
      K4=min(find(d4>r*d0));
      K=[K1 K2 K3 K4];
      if ~isempty(K)
         taufwd(I,J)=sum(tfwd(K))/length(K);
      end
   end
end
FSLEbwd=log(r)./taubwd; FSLEbwd(isinf(FSLEbwd))=NaN;
FSLEfwd=log(r)./taufwd; FSLEfwd(isinf(FSLEfwd))=NaN;

A=FSLEbwd; A(A<1/3*max(max(A)))=NaN;
R=FSLEfwd; R(R<1/3*max(max(R)))=NaN;

figure(5), clf
pcolor(x0(1,2:end-1),y0(2:end-1,1),A)
colormap(BwdOrFwd), shading flat
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (Backward-Time FSLEs)','Interpreter','latex','FontSize',14)

figure(6), clf
pcolor(x0(1,2:end-1),y0(2:end-1,1),R)
colormap(BwdOrFwd), shading flat
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('aLCSs (Forward-Time FSLEs)','Interpreter','latex','FontSize',14)

figure(7), clf
pcolor(x0(1,2:end-1),y0(2:end-1,1),A), hold on
pcolor(x0(1,2:end-1),y0(2:end-1,1),-R), hold off
colormap(BwdAndFwd), shading flat
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (bleu) \& aLCSs (red)','Interpreter','latex','FontSize',14)


%% PERTURBED MANIFOLDS.
figure(8), clf

% Poincare section.
x0=[linspace(-.9,-.15,5) ...
   linspace(-.125,.125,50) ...
   linspace(.15,.9,5)];
y0=zeros(size(x0));
xy0=[x0;y0];
[t xy]=ode45(@duffeq,0:T:100*T,xy0(:),...
   options,...
   'forward',vep,T);
x=xy(2:end,1:2:end-1);
y=xy(2:end,2:2:end);
plot(x,y,'.k','MarkerS',1), hold on

% Hyperbolic point.
x0=-vep/(1+(2*pi/T)^2); 
y0=0;
plot(x0,y0,'k+')

% Compute manifolds.
xs=[]; ys=[];
xu=[]; yu=[];
% IC's on the stable direction (STABLE MANIFOLD).
xy0=[x0+v;y0+v];
for n=2:Nt
   [t xy]=ode23(@duffeq,(n-2)*T:T:n*T,xy0(:),...
      options,...
      'forward',vep,T);
   x=xy(end,1:2:end-1);
   y=xy(end,2:2:end);
   [x y]=manifold(x(:),y(:),ds/(n-1),'spline'); xs=[xs;x]; ys=[ys;y];
   plot(x,y,'r','LineW',1)
   xy0=[x';y'];
end
xy0=[x0-v;y0-v];
for n=2:Nt
   [t xy]=ode23(@duffeq,(n-2)*T:T:n*T,xy0(:),...
      options,...
      'forward',vep,T);
   x=xy(end,1:2:end-1);
   y=xy(end,2:2:end);
   [x y]=manifold(x(:),y(:),ds/(n-1),'spline'); xs=[xs;x]; ys=[ys;y];
   hs=plot(x,y,'r','LineW',1);
   xy0=[x';y'];
end
% IC's on the unstable direction (UNSTABLE MANIFOLD).
xy0=[x0+v;y0-v];
for n=2:Nt
   [t xy]=ode23(@duffeq,(n-2)*T:T:n*T,xy0(:),...
      options,...
      'backward',vep,T);
   x=xy(end,1:2:end-1);
   y=xy(end,2:2:end);
   [x y]=manifold(x(:),y(:),ds/(n-1),'spline'); xu=[xu;x]; yu=[yu;y];
   plot(x,y,'b','LineW',1)
   xy0=[x';y'];
end
xy0=[x0-v;y0+v];
for n=2:Nt
   [t xy]=ode23(@duffeq,(n-2)*T:T:n*T,xy0(:),...
      options,...
      'backward',vep,T);
   x=xy(end,1:2:end-1);
   y=xy(end,2:2:end);
   [x y]=manifold(x(:),y(:),ds/(n-1),'spline'); xu=[xu;x]; yu=[yu;y];
   hu=plot(x,y,'b','LineW',1);
   xy0=[x';y'];
end, hold off

axis([-1.5 1.5 -1 1])
legend([hu hs],{'$W_s$';'$W_u$'},'Interpreter','latex','FontSize',14)
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('$\dot{x} = y$, $\dot{y} = x-x^3+\varepsilon\cos\omega t$','Interpreter','latex','FontSize',14)

%% PERTURBED MANIFOLDS (FTLE).
N=250;
x0=linspace(-1.5,1.5,N); dx0=min(diff(x0));
y0=linspace(-1,1,N); dy0=min(diff(y0));
[x0 y0]=meshgrid(x0,y0);
xy0=[x0(:)';y0(:)'];

% Attracting LCS's.
[t xy]=ode23(@duffeq,[0 T Nt*T],xy0(:),...
   options,...
   'backward',vep,T);
x=xy(end,1:2:end-1); x=reshape(x,[N N]);
y=xy(end,2:2:end); y=reshape(y,[N N]);
dxdx0bwd=prmy(difx(x)/dx0); 
dxdy0bwd=prmx(dify(x)/dy0);
dydx0bwd=prmy(difx(y)/dx0); 
dydy0bwd=prmx(dify(y)/dy0); 

% Repelling LCS's.
[t xy]=ode23(@duffeq,[0 T Nt*T],xy0(:),...
   options,...
   'forward',vep,T);
x=xy(end,1:2:end-1); x=reshape(x,[N N]);
y=xy(end,2:2:end); y=reshape(y,[N N]);
dxdx0fwd=prmy(difx(x)/dx0); 
dxdy0fwd=prmx(dify(x)/dy0);
dydx0fwd=prmy(difx(y)/dx0); 
dydy0fwd=prmx(dify(y)/dy0); 

FTLEbwd=zeros(N-1);
FTLEfwd=zeros(N-1);
for I=1:(N-1)^2
   dxydxy0=[dxdx0bwd(I) dxdy0bwd(I)
            dydx0bwd(I) dydy0bwd(I)];
   FTLEbwd(I)=norm(dxydxy0,2);
   dxydxy0=[dxdx0fwd(I) dxdy0fwd(I)
            dydx0fwd(I) dydy0fwd(I)];
   FTLEfwd(I)=norm(dxydxy0,2);
end
FTLEbwd=log(abs(FTLEbwd))/Nt/T;
FTLEfwd=log(abs(FTLEfwd))/Nt/T;

A=FTLEbwd; A(A<2/3*max(max(A)))=NaN;
R=FTLEfwd; R(R<2/3*max(max(R)))=NaN;

x0=min(min(x0))+dx0/2:dx0:max(max(x0))-dx0/2;
y0=min(min(y0))+dy0/2:dy0:max(max(y0))-dy0/2;

figure(9), clf
imagesc(x0,y0,A), hold on
colormap(BwdOrFwd), shading interp
plot(xs,ys,'k.','MarkerS',1), hold off
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (Backward-Time FTLEs)','Interpreter','latex','FontSize',14)

figure(10), clf
imagesc(x0,y0,R), hold on
colormap(BwdOrFwd), shading interp
plot(xu,yu,'k.','MarkerS',1), hold off
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('aLCSs (Forward-Time FTLEs)','Interpreter','latex','FontSize',14)

figure(11), clf
pcolor(x0,y0,A), hold on
pcolor(x0,y0,-R), hold off
colormap(BwdAndFwd), shading interp
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (bleu) \& aLCSs (red)','Interpreter','latex','FontSize',14)

%% PERTURBED MANIFOLDS (FSLE).
r=125;
tspan=linspace(0,100*T,50);
N=350;
x0=linspace(-1.5,1.5,N); d0=min(diff(x0));
y0=linspace(-1,1,N); 
[x0 y0]=meshgrid(x0,y0);
xy0=[x0(:)';y0(:)'];

% Attracting LCS's.
[tbwd xy]=ode23(@duffeq,tspan,xy0(:),...
   options,...
   'backward',vep,T);
xbwd=xy(:,1:2:end-1); xbwd=reshape(xbwd,[length(tbwd) N N]);
ybwd=xy(:,2:2:end); ybwd=reshape(ybwd,[length(tbwd) N N]);

% Repelling LCS's.
[tfwd xy]=ode23(@duffeq,tspan,xy0(:),...
   options,...
   'forward',vep,T);
xfwd=xy(:,1:2:end-1); xfwd=reshape(xfwd,[length(tfwd) N N]);
yfwd=xy(:,2:2:end); yfwd=reshape(yfwd,[length(tfwd) N N]);

taubwd=repmat(0,[N-2 N-2]);
taufwd=repmat(0,[N-2 N-2]);
for I=1:N-2
   for J=1:N-2
      d1=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I,J+1)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I,J+1)).^2);
      d2=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I+1,J)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I+1,J)).^2);
      d3=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I+2,J+1)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I+2,J+1)).^2);
      d4=sqrt((xbwd(:,I+1,J+1)-xbwd(:,I+1,J+2)).^2+(ybwd(:,I+1,J+1)-ybwd(:,I+1,J+2)).^2);
      K1=min(find(d1>r*d0));
      K2=min(find(d2>r*d0));
      K3=min(find(d3>r*d0));
      K4=min(find(d4>r*d0));
      K=[K1 K2 K3 K4];
      if ~isempty(K)
         taubwd(I,J)=sum(tbwd(K))/length(K);
      end
      d1=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I,J+1)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I,J+1)).^2);
      d2=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I+1,J)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I+1,J)).^2);
      d3=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I+2,J+1)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I+2,J+1)).^2);
      d4=sqrt((xfwd(:,I+1,J+1)-xfwd(:,I+1,J+2)).^2+(yfwd(:,I+1,J+1)-yfwd(:,I+1,J+2)).^2);
      K1=min(find(d1>r*d0));
      K2=min(find(d2>r*d0));
      K3=min(find(d3>r*d0));
      K4=min(find(d4>r*d0));
      K=[K1 K2 K3 K4];
      if ~isempty(K)
         taufwd(I,J)=sum(tfwd(K))/length(K);
      end
   end
end
FSLEbwd=log(r)./taubwd; FSLEbwd(isinf(FSLEbwd))=NaN;
FSLEfwd=log(r)./taufwd; FSLEfwd(isinf(FSLEfwd))=NaN;

A=FSLEbwd; A(A<1/3*max(max(A)))=NaN;
R=FSLEfwd; R(R<1/3*max(max(R)))=NaN;

figure(12), clf
pcolor(x0(1,2:end-1),y0(2:end-1,1),A), hold on
colormap(BwdOrFwd), shading flat
%plot(xs,ys,'k.','MarkerS',1), hold off
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (Backward-Time FSLEs)','Interpreter','latex','FontSize',14)

figure(13), clf
pcolor(x0(1,2:end-1),y0(2:end-1,1),R), hold on
colormap(BwdOrFwd), shading flat
%plot(xu,yu,'k.','MarkerS',1), hold off
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('aLCSs (Forward-Time FSLEs)','Interpreter','latex','FontSize',14)

figure(14), clf
pcolor(x0(1,2:end-1),y0(2:end-1,1),A), hold on
pcolor(x0(1,2:end-1),y0(2:end-1,1),-R), hold off
colormap(BwdAndFwd), shading flat
axis([-1.5 1.5 -1 1])
xlabel('$x$','Interpreter','latex','FontSize',14)
ylabel('$y$','Interpreter','latex','FontSize',14)
title('rLCSs (bleu) \& aLCSs (red)','Interpreter','latex','FontSize',14)