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

% Here we test out Newton's method:

% f(x) = x^5 - 1/2
>> type f.m
function value = f(x)
value = x.^5 - 0.5;
end

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

% here is a simple implementation, which actually stops based on
% an error in y, not in x (we'll talk about a way of doing that later)
>> type newt.m
function x = newt(x, err)
while  abs(f(x)) > err
    x = x - f(x)/df(x);
    i = i + 1;
end

% We see f has a root at around 0.9
>> x = 0:0.01:1;
>> plot(x,f(x),x,0)

% we see that it converges to (approximately) the same place
% when we start it in different places

>> newt(0.5, 1e-10)
ans = 0.87055056329613

>> newt(2, 1e-10)
ans = 0.87055056329613

>> newt(20, 1e-10)
ans = 0.87055056330047

% if we use a smaller tolerance they agree

>> newt(20, 1e-15)
ans = 0.87055056329612

>> newt(2, 1e-15)
ans = 0.87055056329612

>> newt(0.01, 1e-15)
ans = 0.87055056329612

>> newt(-1, 1e-15)

ans = 0.87055056329612

% Now we change newt.m so that it outputs a vector of iterates:

>> type newt.m

function s = newt(x, err)
i = 1;
s = [];
s(i) = x;
while  abs(f(x)) > err
    x = x - f(x)/df(x);
    i = i + 1;
    s(i) = x;
end

% And the same for bisect.m:
>> type bisect.m

function s = bisect(a, b, ep)
i = 1;
s = [];
c=(a+b)/2;
s(i) = c;
while (b-c >= ep)
    if sign(f(b))*sign(f(c)) <= 0
        a=c;
    else
        b=c;
    end
    c=(a+b)/2;
    i = i + 1;
    s(i) = c;
end

% We can see how Newton's method converges in this case in a graph
>> xn = newt(0.5, 1e-16);
>> plot(1:11,xn,'x-')

>> xb = bisect(0,1,1e-8)

% We plot both together:
>> plot(1:11,xn,'x-', 1:27,xb,'.-r')

% The value of the root:
>> r = (1/2)^(1/5)
r = 0.87055056329612

% Here we plot the absolute value of the error in both.  It's hard
% to see what's going on where it converges...
>> plot(1:11,abs(xn - r),'x-', 1:27,abs(xb - r),'.-r')

>> ne = abs(xn - r);
>> be = abs(xb - r);

% If we plot a subrange, it doesn't help much...
>> sv = 7;
>> plot(sv:11,ne(sv:11),'x-', sv:27,be(sv:27),'.-r')

% Differences become clear in a semilog plot:
>> plot(1:11,log10(ne),'x-', 1:27,log10(be),'.-r')

% Let's look at other initializations of Netwon:

>> ne = abs(newt(0.6,1e-15) - r);
>> ne2 = abs(newt(-1,1e-15) - r);

% And we can compare these in a semi-log plot:
>> plot(1:8,log10(ne),'x-', 1:27,log10(be),'.-r', 1:34, log10(ne2), '.-g');

% And a semilog-log plot:
>> plot(1:8,log10(abs(log10(ne))),'x-', 1:27,log10(abs(log10(be))),'.-r', 1:34, log10(abs(log10(ne2))), '.-g');