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

format long;

n = 10;
% n = 20;
% n = 40;
% n = 73; % near overflow, stay below it.

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

  a(1) = 1; a(2) = 1/3 ; b = 13/3; c=-4/3; % The case in the slides
%       a(1) = 1/3; a(2) = 1/9; b = 10/3; c=-1; % The case in textbook
% a0 = 1; a(1) = 1/3;             b = 10/3; c=-1; a(2) = b*a(1)+c*a0; %a(2)=1/9
% and more cases
% 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/4; b = 13/3; c=-4/3;

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 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_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) 

% 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) )