
% Iterating a linear 2D map with distinct real eigenvalues
% for contraction/expansion of an e-ball: Ex 2.5, fig 2.12

clear
MTE=0;

map = inline('[A(1,1)*x+A(1,2)*y,A(2,1)*x+A(2,2)*y]','A','x','y');

%%distinct real eig:
a=1.25; b=0.65; A12=0; A21=0; L=0.65;
%%repeated real eig:
%a=0.5; b=0.5; A12=0.5; A21=0.0; L=0.3;   % contracting
%a=1.25; b=1.25; A12=0.5; A21=0.0; L=0.3; % expanding
%%distinct real eig with rotation:
%a=1.25; b=0.65; A12=0.0; A21=0.5; L=0.3;
%%complex eigenvalues 
%a=1.05; b=1.05; A12=-1; A21=1; L=0.35;
%%complex eigenvalues with rotation
%a=0.8; b=0.8; A12=-2; A21=1; L=0.5;

%%matrix times circle equals ellipse:
%a=2; b=1; A12=0.5; A21=1.5; L=0.5;MTE=1;

A = [a A12;A21 b]

[V,D]=eig(A);
eigenvalues=diag(D)
eigenvectors=V

eigenvalues_abs=abs(diag(D))

%fixed point:
xs = 0;
ys = 0;

figure(1)
clf
lw=3;
ms1=20;
ms2=25;

plot(xs,ys,'r+','LineWidth',lw)
axis([-L L -L L]); 
hold on
set(gca,'FontSize',[18]);
xlabel('x_n'); ylabel('y_n');
%title(['Linear map: A=',num2str(A)])
title(['Linear map'])

% build circles around f.p.
Ncircles=1;
e=0.2;
Npts=200;
Npts1=8;
ang=0:2*pi/Npts:2*pi;
ang=ang';
x= xs+e*cos(ang);
y= ys+e*sin(ang);
x1=x(1:Npts/Npts1:Npts);
y1=y(1:Npts/Npts1:Npts);

mycols=jet(Npts1);
 
plot(x,y,'r-','LineWidth',lw)
plot(x1(1),y1(1),'*','MarkerSize',ms1)
for jj=1:Npts1
 plot([x1(jj),x1(jj)]',[y1(jj),y1(jj)]','.','MarkerSize',ms2,'Color',mycols(jj,:))
end
axis equal
axis([-L L -L L]); 
drawnow

for i=1:Ncircles
 fprintf('Press any key to continue...\n');
 pause
 Y=map(A,x,y);
 Y1=map(A,x1,y1);
 x=Y(:,1);
 y=Y(:,2);
 x1=Y1(:,1);
 y1=Y1(:,2);
 if lw>1  
  lw=lw-1;
 end
 plot(x,y,'r-','LineWidth',lw)
 plot(x1(1),y1(1),'*','MarkerSize',ms1)
 for jj=1:Npts1
  plot([x1(jj),x1(jj)]',[y1(jj),y1(jj)]','.','MarkerSize',ms2,'Color',mycols(jj,:))
 end
 axis equal
 axis([-L L -L L]); 
 drawnow
end


if(MTE) %% IF Matrix x Circle = Ellipse:
 [Vt,Dt]=eig(A*transpose(A))
 [V,D]=eig(A)
 plot(e*sqrt(Dt(1,1))*[-1,1]*Vt(1,1)./norm(Vt(:,1)),e*sqrt(Dt(1,1))*[-1,1]*Vt(2,1)./norm(Vt(:,1)),'r-','LineWidth',2)
 plot(e*sqrt(Dt(2,2))*[-1,1]*Vt(1,2)./norm(Vt(:,1)),e*sqrt(Dt(2,2))*[-1,1]*Vt(2,2)./norm(Vt(:,1)),'b-','LineWidth',2)
 e2=e*10;
 plot(e2*sqrt(D(1,1))*[-1,1]*V(1,1)./norm(V(:,1)),e2*sqrt(D(1,1))*[-1,1]*V(2,1)./norm(V(:,1)),'k-','LineWidth',1)
 plot(e2*sqrt(D(2,2))*[-1,1]*V(1,2)./norm(V(:,1)),e2*sqrt(D(2,2))*[-1,1]*V(2,2)./norm(V(:,1)),'k-','LineWidth',1)
 axis([-1 1 -1 1]*0.45)
 drawnow
end
hold off

