format short

clc

% compute and store n iterations of a(n) = a(n-2) - a(n-1)
n = 100;
% the sequence is a matrix of size 1 x (n+1), i.e. a row vector
numericalSeq = zeros(1,n+1); 
% set the first two terms
numericalSeq(1) = 1;
numericalSeq(2) = 0.5*(sqrt(5)-1);

for i = 3:n+1
    numericalSeq(i) = numericalSeq(i-2) - numericalSeq(i-1);
end

for i = 1:n+1
    exactSeq(i) = ((sqrt(5)-1)/2.0)^(i-1);
end

% compute a vector with the (absolute) errors
diffSeq = numericalSeq - exactSeq;

% diff is the error, rel is the relative error
disp('** First 5 terms')
for i = 1:5
    disp(['i=' num2str(i) ...
        ' numeric ', num2str(numericalSeq(i)) ...
        ' exact ' num2str(exactSeq(i)) ...
         ' diff ' num2str(diffSeq(i)) ...
         ' rel ' num2str(diffSeq(i) / exactSeq(i))])

end

% ... un-comment to explore

% disp('** Terms 36...45')
% for i = 36:45
%     disp(['i=' num2str(i) ...
%         ' numeric ', num2str(numericalSeq(i)) ...
%         ' exact ' num2str(exactSeq(i)) ...
%          ' diff ' num2str(diffSeq(i)) ...
%          ' rel ' num2str(diffSeq(i) / exactSeq(i))])
% 
% end
% 
% disp('** Terms 81...85')
% for i = 81:85
%     disp(['i=' num2str(i) ...
%         ' numeric ', num2str(numericalSeq(i)) ...
%         ' exact ' num2str(exactSeq(i)) ...
%          ' diff ' num2str(diffSeq(i)) ...
%          ' rel ' num2str(diffSeq(i) / exactSeq(i))])
% 
% end
% 
% % be aware that log is natural log (ln) here
% % also the we take matrix transposes using ' so that plotData 
% % has the sequences in columns
% plotData = [ log(abs(numericalSeq))' , log(abs(exactSeq))' ];
% % plot the first 80 rows of plotData (for all columns)
% plot( plotData( 1:80 , : ), 'Linewidth', 5 )
% legend( 'numerical', 'exact' )
% xlabel('i')
% ylabel('log( |a_i| )')
% 
% % To understand what is going on, remember that we set 
% % a1 = ( sqrt(5)-1 ) / 2.0
% % but this number is likely not representable
% % 
% % let's compute the exact solution but assuming a small error in the 
% % initial data
% % 
% exactSeqErr = zeros(1,n+1);
% beta = 1e-16;
% alpha = 1-beta;   %% this ensures a_0 = 1 (up to machine epsilon)
% 
% for i = 1:n+1
%     exactSeqErr(i) = ((sqrt(5)-1)/2.0)^(i-1) + ...
%                    beta*((-sqrt(5)-1)/2.0)^(i-1);
% end
% format long
% disp('** these are the first few terms for the exact and comparison solution')
% for i=1:3
%   % note %.16e means 16 decimal places...
%   disp([ num2str(exactSeqErr(i),'%.16f'), ' ', ...
%          num2str(exactSeq(i),'%.16f') ])
% end
% 
% plotData = [ log(abs(numericalSeq))' , log(abs(exactSeq))', ...
%     log(abs(exactSeqErr))' ];
% plot( plotData( 1:80 , : ) , 'Linewidth', 5 )
% legend( 'numerical', 'exact', 'exact with modified initial data' )
% xlabel('i')
% ylabel('log( |a_i| )')


