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