% ROUNDOFF_ERROR is an m-file to demonstrate how roundoff errors
% can add up in numerical calculations.
%
% Consider the algebraic equations:
% phi^(n+1) = phi^(n-1)-phi^n
% and note that phi^(n+2) = phi^(n+1) * phi.
% We can calculate the powers of phi in two ways:
% by taking phi to the n-th power, or
% iteratively.
% 
% The roundoff errors grow differently in both methods.
%
% This routine outputs a table with four rows. The first row contains n,
% the second the solution calculated by raising phi to the nth power, 
% the third contains the iterative solution and the forth the relative
% error (difference between both solutions divided by magnitude of the 
% power solution).
%
% The results are placed in the variable R in the workspace.

% Initialize some variables
phi=(sqrt(5)-1)/2;

power=0:50;

% Take the nth power in vectorized form:
P=phi.^power;

% Now use the recursive algorithm in a "for" loop
Pr(1)=1;   % =phi^0
Pr(2)=phi; % =phi^1
for i=3:51
   Pr(i)=Pr(i-2)-Pr(i-1);
end

% Calculate the relative differences between both solutions
Pe=(P-Pr)./P;

% Change the display format to "long"
format long

% Put it together and display on the screen as a table
R=[power' P' Pr' Pe']

format short
