% Ex1stability.m
% Example of stability vs instability in a sequence
% In the following examples, the sequence y(n) = int_0^1 x^n exp(-x) dx
% is considered. For all n, y(n)>0 and y(n) -> 0 as n increases.
% INSTABILITY
% The following sequence is unstable because the error in
% y(n-1) is magnified (by multiplying by n). The first few
% terms are ok, but the error grows very quickly.
y(1) = 1 - exp(-1);
fprintf(' %d %25.15e\n',0,y(1))
for n=2:30
y(n) = -exp(-1) + (n-1)*y(n-1);
fprintf(' %d %25.15e\n',n-1,y(n))
end
% STABILITY
% The following sequence is stable because the error in
% y(n) is reduced (by dividing by n). To compute a certain
% value of y(n), say y(20), start with a much higher n value,
% say, y(100)=0. This value is certainly wrong, so it has a
% big error. But as iterations proceed, the error is reduced
% so that after 80 iterations, the error is essentially zero.
for n = [100 70 50 40 30 25]
disp(' ')
clear y
y = zeros(1,n);
y(n) = 1;
factor = 1;
%fprintf(' %d %25.15e\n',100,y(100))
for k=n:-1:21
y(k-1)=(y(k)+exp(-1))/k;
factor = factor/k;
% fprintf(' %d %25.15e\n',k-1,y(k-1))
end
fprintf('With n=%d, y(20) is %20.15f \n',n,y(20))
fprintf('The diminishment factor is %15.6e\n \n',factor)
end