A2 Scripts

P1

% Script P1: Tests NearestRepeat

randn('seed',0);
n = 8;

clc
A = randn(n,n);
[B_opt,dist] = NearestRepeat(A)

P2

% P2: Tests BreakDownStep  

clc
randn('seed',0)
n = 100;
disp(' Example   alpha_min      BreakDownStep     sigma_min')
disp('-------------------------------------------------------')
for eg = 1:20
    A = randn(n,n); z = sum(abs(A)); A = A + 2*diag(z);
    u = randn(n,1);
    v = randn(n,1);
    [alpha_min,BreakDownStep] = BreakLU(A,u,v);
    B = A + alpha_min*u*v';
    sigma_min = min(svd(B(1:BreakDownStep,1:BreakDownStep)));
    disp(sprintf(' %2d    %20.15f   %3d            %10.3e',eg,alpha_min,BreakDownStep,sigma_min))
end

P3

% Script P3: Tests TriProdShiftSol

n = 8;
randn('seed',0);
S = triu(randn(n,n));
T = triu(randn(n,n));
lambda = 10;
b = randn(n,1);
x = TriProdShiftSol(S,T,lambda,b);
x0 = (S*T - lambda*eye(n,n))\b;
clc
disp('   TriProdShiftSol        (S*T-lambda*I)\b')
disp('--------------------------------------------------')
for k=1:n
    disp(sprintf('%20.15f   %20.15f',x(k),x0(k)))
end

P4

% Script P4; Tests RayThruPolygon

% Polygon vertices...
x = [0.8249   0.7604   1.2581   1.3963   1.9954   2.4194   2.5023   3.2028   3.3963   3.4608 ...
     2.5760   2.4194   2.0968   2.0323   1.4885   1.5714   1.1014   1.1751   0.6037   0.3180];

y = [ -0.3801  -1.1871   -0.8246   -1.3977   -0.6257   -1.6082   -0.7076   -0.0175   -0.9766    0.8246  ...  
       1.5965  -0.2164    1.1871   -0.0526   -0.6257    0.3333    0.7661   -0.1111    0.3567   -0.2865];
A = [x;y];
v = [0;1];
close all
for theta = linspace(-pi/2,pi/3,10);
    figure
    w = [cos(theta);sin(theta)];
    plot([x x(1)],[y y(1)],[v(1) v(1)+5*w(1) ],[v(2) v(2)+5*w(2)],'r',v(1),v(2),'*r')
    axis([-1 4 -2 2])
    Hits = RayThruPolygon(A,v,w);
    [m,k] = size(Hits);
    if k>0
        hold on
        for j=1:k
            plot(Hits(1,j),Hits(2,j),'*g')
        end
        m = n/2;
        d = sum(sqrt((Hits(1,1:2:k)-Hits(1,2:2:k)).^2 + (Hits(2,1:2:k)-Hits(2,2:2:k)).^2 ));
        title(sprintf('Intersection distance = %10.6f',d))
        hold off
    else
        title(sprintf('Intersection distance = %10.6f',0))
    end
    pause
end