Assignment 2 Scripts

P1

% Script P1
% Tests SchurCompSolve

clc
randn('seed',0)

% Generate a strictly diagonally dominant A'

n = 500; 
A = 2*eye(n,n);
for k=1:n-1
   A(k,k+1) = -1;
   A(k+1,k) = -1;
end

% Ignore the fact that this matrix is tridiagonal

range = 50:50:n;
error = [];
for m = range
   c = randn(m,1);
   % Horrible Explicit Method...
   i = n-m+1:n; j = 1:n-m;
   z = (A(i,i) - A(i,j)*inv(A(j,j))*A(j,i))\c;
   % Compare with SchurCompSolver...
   z0 = SchurCompSolver(A,c);
   error = [error norm(z-z0)/norm(z)];
end
% Plot the relative discrepancy between the two methods
semilogy(range,error+eps)
title('SchurCompSolve Behavior (n = 500)')
xlabel('m')
ylabel('Approx Relative Error')

P3

% Script P3
% Tests Collection Solve

A = [   5         6         7         7;...
       11         6        42        23;...
      593       709       639       764;...
      5186     6855      3129      4136;...
     19074    10639     15119      8433;...
   804700    962113    867123   1036747;...
   5293759   6997439   5704410   7540249];

B = [A(:,1)+A(:,2) A(:,3)+A(:,4)];
   
clc
X = CollectionSolve(A,B)
X = CollectionSolve(A/10,B/10)

P4

% Script P4 
% Tests DerDet

clc
randn('seed',0);

%Generate a small example

n = 10;
A = randn(n,n); 
for k=1:n
   A(k,k) = 1 + sum(abs(A(1:k-1,k))) + sum(abs(A(k+1:n,k)));
end
A = A/norm(A);

% Compare DerDet to a crude method that works with
% divided differences: beta = approximately (f(delta)-f(0))/delta

delta = .00000001;
d = det(A);

disp('     k        DerDet(A,k)          (f(delta)-f(0))/delta')
disp('----------------------------------------------------------')
for k=1:n
   beta = DerDet(A,k);
   Atilde = A; Atilde(k,k) = Atilde(k,k) + delta;
   betaTilde = (det(Atilde)-d)/delta;
   disp(sprintf('    %2d   %20.12e      %20.12e  ',k,beta,betaTilde))
end