%% Annotated Matlab transcript -- 2/22
%% CS 323 - Doug DeCarlo

%% Looking into Newton's method...

% We're still looking at f(x) = x^5 - 1/2
>> type f.m
function value = f(x)
value = x.^5 - 0.5;
end

% Plot error in iterates of Newton's method

>> xn = newt(100, 1e-16);
>> r = (1/2)^(1/5);

% in normal plot - hard to see what's going on
>> plot(1:length(xn),r-xn, '.')

% in semilog plot - you can see when it starts to converge quadratically
>> plot(1:length(xn),log10(abs(r-xn)), '.')

% in semiloglog plot, the quadratically convergent area is a line
% (as the function is x^(2^i))
>> plot(1:length(xn),log10(abs(log10(abs(r-xn)))), '.')

% here are the 1st and 2nd derivatives of f(x):

>> type df.m
function value = df(x)
value = 5*x.^4;
end

>> type ddf.m
function value = ddf(x)
value = 20*x.^3;
end

% we can look at the value of M (which we evaluate at the root)

>> -ddf(r)/(2*df(r))

ans = -2.29739670999407

>> -1/(ddf(r)/(2*df(r)))

ans =  -0.43527528164806

% and we expect that a guess between the root and here will converge
% (also below the root, but in this example, it approaches from above)

>> r + 1/(ddf(r)/(2*df(r)))

ans = 1.30582584494419

% and if we look at the graph, it seems to be converging well around the
% 20th and 21st iterations, and this is indeed in the right area
 
>> xn(21)

ans = 1.18641606662474

% we can see how much M varies here -- so this gives us an idea of
% how approximate statements that use it really are; it varies by around
% a factor of 2 in this area

>> x = 0.5:0.01:1.5;
>> plot(x,ddf(x)./(2*df(x)))



%% Now we look at the error estimates: x(n+1) - xn

% we use the diff command in Matlab which finds differences between
% adjacent elements in a vector

>> diff([1,5,6,7,15])
ans =  4     1     1     8

% Now we can compare the error r-xn with the error estimate
>> plot(1:length(xn), r-xn, 1:(length(xn)-1), diff(xn))

% and we can compare them better in a semilog plot
% here, we see how they agree very closely in the convergent region

>> plot(1:length(xn), log(abs(r-xn)), '.-', 1:(length(xn)-1), log(abs(diff(xn))), 'x-');