echo on
%%%%% ch04_exa2.m %%%%%
% Run script and study the effects of each line.

echo off
format compact;
format short;
disp(sprintf('**** Additional Examples for Chapter 04 ****'));


echo off
disp(sprintf('**** Example 4.x Condition number, norm, eigenvalues, eigenvectors*'));
echo on
%%%% Draw a unit circle and A*vector from the unit circle.
tol=100		% sensitivity of the algorithm
A=[1 2; 1 3]
%A=[1 0;0 3]
condA = cond(A)		% condition number
norm(A)*norm(inv(A))	%by definition
[V,D]=eig(A)	% eigenvectors and eigenvalues
t=[-1:.01:1];

%Calculate vectors from the unit circle
x=cos(t*pi);
y=sin(t*pi);

%Multiply these vectors by A, and calculate norm of the resulting vector
echo off
x1=zeros(1,size(x,2));
y1=zeros(1,size(x,2));
n1=zeros(1,size(x,2));
for i=1:size(x,2)
   x1(i)=A(1,:)*[x(i),y(i)]';
   y1(i)=A(2,:)*[x(i),y(i)]';
   n1(i)=norm([x1(i),y1(i)]);
end
echo on
max(n1)	% maximal norm

% The same operation with the inverse of A
B=inv(A)
echo off
x2=zeros(1,size(x,2));
y2=zeros(1,size(x,2));
n2=zeros(1,size(x,2));
for i=1:size(x,2)
   x2(i)=B(1,:)*[x(i),y(i)]';
   y2(i)=B(2,:)*[x(i),y(i)]';
   n2(i)=norm([x2(i),y2(i)]);
end
echo on
max(n2)	% maximal norm
max(n2)*max(n1)	% condition number

% Draw vectors from unit circle x (blue), A*x (green) and inv(A)*x (black)
% Calculate and denote the direction of the eigenvectors (red O and +) and 
% values for the corresponding eigenvalues.
echo off
hold on	
for i=1:size(x,2)
	origVec = [x(i) y(i)]';		% check if vectors have same direction
	multOrthVec = [-y1(i) x1(i)]';
	if (abs(origVec' * multOrthVec) < condA/tol)	% check if this is an eigenvector
		plot(x(i),y(i),'ro');
      plot(x1(i),y1(i),'ro');
      lambda1=x1(i)/x(i)
      eig1=[x(i),y(i)]
   else
		plot(x(i),y(i),'b.');
      plot(x1(i),y1(i),'g.');
   end
   origVec = [x(i) y(i)]';		% check if vectors have same direction
   multOrthVec = [-y2(i) x2(i)]';
   if (abs(origVec' * multOrthVec) < condA/tol)	% check if this is an eigenvector
		plot(x(i),y(i),'r+');
      plot(x2(i),y2(i),'r+');
      lambda2=x(i)/x2(i)
      eig2=[x(i),y(i)]
   else
      plot(x2(i),y2(i),'black.');
   end
end
hold off