% V02 of generalized_Fibonacci, written by Wei-Cheng Wang.
% to demonstrate stability and instability of recurrence formula

format long;

n = 40;
% n = 73;

a = zeros(1,n);
ar = zeros(1,n);
a_true = zeros(1,n);

a(1) = 1/5; a(2) = 1/10; b = 9/2; c = -2;
a(1) = 1/5; a(2) = 1/5 ; b = 9/2; c = -2;
a(1) =   1; a(2) = 1/2 ; b = 9/2; c = -2;
% a(1) = 1; a(2) = 1/2; b = 5/2; c = -1;
% a(1) = 1; a(2) = 1.01; b = 1; c = 1;
% a(1) = 1; a(2) = 1/3 ; b = 13/3; c=-4/3; % This is case A
% a(1) = 1; a(2) = 1/4; b = 13/3; c=-4/3;

for i = 3:n
  a(i) = b*a(i-1) + c*a(i-2);
end

% true solution: a_true(i) = c1*lambda_1^i + c2*lambda_2^i;
% c1, c2 can be obtained from a(1), a(2)
% or from a(0) and a(1).
% Here is the formula using a(0) and a(1).

a0 = ( a(2)-b*a(1) )/c ;  
lambda_1 = (b - sqrt(b^2+4*c))/2; 
lambda_2 = (b + sqrt(b^2+4*c))/2; 
c1 = (a0*lambda_2 - a(1))/(lambda_2-lambda_1); % need not plug lambda_1, 
                                               % lambda_2 into the formula
c2 = (a(1) - a0*lambda_1)/(lambda_2-lambda_1); % for c1, c1. It is cleaner 
                                               % this way and easier to debug.
lambda = [lambda_1, lambda_2]
c1
c2
% c2 = 0  % only for case A, overwrite c2, make sure c2 = 0 to get correct result
for i = 1:n
  a_true(i) = c1*lambda_1^i + c2*lambda_2^i;
end
% double check if the above calculation is correct.
check_a1 = a(1)-a_true(1)
check_a2 = a(2)-a_true(2)

% relative error for a(n)
check_an_rel_err = ( a(n)-a_true(n) )/a_true(n) 

rel_err = abs( (a-a_true)./a_true );
figure(1), plot(1:n,abs(a-a_true),'x'), title('true error')
figure(2), plot(1:n,rel_err,'x'), title('relative error')

% ar(n) = a_true(n);
% ar(n-1) = a_true(n-1);
ar(n) = a(n);
ar(n-1) = a(n-1);
for i = n-2:-1:1
  ar(i) = ( ar(i+2)-b*ar(i+1) )/c ;  
end

% check the results for backward recursion to a(1)
check_a1_backwards = ( ar(1)-a_true(1) )
check_a1_backwards_rel_err = abs( ( ar(1)-a_true(1) )/a_true(1) )