Assignment 1 Scripts and Functions

x = -20;
sum = 1; k = 1; term = x; next = sum + term;
while next~=sum
   sum = next; k = k+1; term = x*term/k; next = sum + term;
end
format long
sum = next
exact = exp(x)

  function v = Recover(A,v0,k)
% A is an n-by-n nonsingular matrix, v0 is a column n-vector, and
% k is a positive integer.
% Multiplies v0 by A a total of k times and then attempts to recover
% v0 by solving a sequence of k linear systems involving A.

for j=1:k
   v0 = A*v0;
end

for j=1:k
   v0 = A\v0;
end
v = v0;

% P4 The Effect of Condition

n = 100;
% Generate a random Markov matrix...
rand('seed',0)  % ensures that we all solve the same problem
A = rand(n,n);

% Normalize the columns...
for k=1:n 
    A(:,k) = A(:,k)/sum(A(:,k));
end

% Initial vector
v0 = rand(n,1);
v0 = v0/sum(v0);

kappa = cond(A);

clc
disp(sprintf('cond(A) = %10.3e',kappa))
disp('   ')
disp(' k   norm(v0-v1)/norm(v0)')
disp('------------------------------')
for k=1:6
    v1 = Recover(A,v0,k);
    err = norm(v1-v0)/norm(v0);
    disp(sprintf(' %2d     %10.3e',k,err))
end

% P5 Structured Solutions

close all
n = 10;
rand('seed',0)
% Generate a random A...
A = rand(n,n);
for j=1:n
    A(j,j) = 0;
    A(j,j) = -sum(A(:,j));
end

clc
for t=[.5 1 2 4]
    % Plot the error as a function of q...
    F0 = expm(A*t);
    eVec = [];
    for q = 1:16 
        F = Pade(A,t,q);
        AbsErr = norm(F-F0);
        eVec = [eVec AbsErr];
    end
    figure
    semilogy(eVec)
    title(sprintf('Error in Rqq(At) for t = %5.2f',t))
    xlabel('q')
end