% V02 of generalized_Fibonacci, updated 20/09/28 by Wei-Cheng Wang.
% to demonstrate stability and instability of recurrence formula
% pay attention at line 49 for the textbook case.

clear;

format long;

% n = 73 is near overflow, stay below it.
n=50; a=zeros(1,n); ar=zeros(1,n); a_true=zeros(1,n); 

% Need to enforce c2 = 0 only for these cases
% The case in the slides:
% a(1) = 1; a(2) = 1/3 ; b = 13/3; c=-4/3; force_c2_zero = true;
%
% The case in textbook:
%          a(1) = 1/3; a(2) = 1/9; b = 10/3; c=-1; force_c2_zero = true; 
% same as above, in principle (a(2) = 1/9)
% a0 = 1; a(1) = 1/3; b = 10/3; c=-1; a(2) = b*a(1)+c*a0; force_c2_zero = true; 
% more cases with c2=0
% a(1) = 1/5; a(2) = 1/10; b = 9/2; c = -2; force_c2_zero = true;

% Other cases cases, c2 neq 0
% a(1) = 1; a(2) = 1/4; b = 10/3; c=-1; force_c2_zero = false; %only minor difference from above case
% a(1) = 1/5; a(2) = 1/5 ; b = 9/2; c = -2; force_c2_zero = false;
  a(1) = 1; a(2) = 1.1; b = 1; c = 1; force_c2_zero = false;

% Some cases that can computed without error for n = 70 (why?):
% a(1) =   1; a(2) = 1/2 ; b = 9/2; c = -2; force_c2_zero = false;
% a(1) = 1; a(2) = 1/2; b = 5/2; c = -1; force_c2_zero = false;

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

% 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_1 = -2*c/(b + sqrt(b^2+4*c)) %alternative formula for b>0, |c|<<b 
lambda_2 = (b + sqrt(b^2+4*c))/2; 
lambda = [lambda_1, lambda_2]

% Instead of typing detailed formula of lambda_1, lambda_2 to get c1, c2
% use computed lambda_1, lambda_2 in the formula of c1, c2. 
% It is clearer this way and easier to debug.
c1 = (a0*lambda_2 - a(1))/(lambda_2-lambda_1); 
c2 = (a(1) - a0*lambda_1)/(lambda_2-lambda_1);
c1_c2 = [c1 c2]
% c2= 0  %% forcing c = 0 makes a difference for first 2 cases
if ( force_c2_zero == true ) 
 c2 = 0
end

for i = 1:n
  a_true(i) = c1*lambda_1^i + c2*lambda_2^i;
end
an_given_by_solution_formula = a_true(n)

% double check if the above calculation is correct.
check_errors_in_a1_a2 = [ a(1)-a_true(1)  a(2)-a_true(2) ]

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

% history of true error and relative error Vs 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) )