% BEGIN FILE: hlpd034.m
% BEGIN FILE: fxptite2.m
% 17-Sep-2018 
% Want the root (=solution) of  f(x) = 0. 
% Method:  fixed point iteration
% Transform f(x) = 0 in the form x=g(x) and iterate
% x_{n+1}= g(x_n)  which converges if |g'(root)|< 1
% There are many ways to find g(x), but making sure
% |g'(root)|< 1 requires care.
% Trivial choice g(x) = f(x) -x NEED NOT satisfy this.
% Study f(x) = x^2-x-2 with roots 2 and -1 and the choices
% 1. g(x) =x^2-2    2. g(x) = sqrt(x+2)  3. g(x) =1+ 2/x
% 4. g(x) =(x.^2+2)./(2*x -1); 
% f(x) = 0 is then equivalent to g(x) = x
close all; clear;
path(path,'../util')
startup
for ii=1:4
    if ii==1
        myg=@(x)(x.^2-2);
    elseif ii==2
        myg=@(x)(sqrt(x+2));
    elseif        ii==3
        myg=@(x)(1+1./x);
    else
        myg=@(x)(x.^2+2)./(2*x-1);
    end
    figure
    hold on
    exact= 2;
    x=exact-0.5;  % lower bound for plotting range
    ite=[x];
    plot(x:0.02:(x+1), myg(x:0.02:(x+1)), 'k-', ...
        x:0.02:(x+1), (x:0.02:(x+1)), 'b-','LineWidth',1.0  );
    plot([x, x, myg(x)], [x,myg(x), myg(x) ],'r-')
    xlabel('Fixed point iteration')
    grid
    for j=1:10
       x=myg(x);   % fixed point iteration
       plot([x, x, myg(x)], [x,myg(x), myg(x) ],'r-')
       ax=axis;
       text(ax(1)+0.1, 0.95*ax(4)-(j-1)*0.09*(ax(4)-ax(3)),...
           ['x= '   num2str(x,' %4.2f')],'FontSize',14)
       ite=[ite x];
    end
    fprintf('-------\nii= %d,  function g(x) = %s \n',ii, func2str(myg));
    fprintf('Error in odd steps:\n')  
    fprintf('%12.4e',(ite(1:2:end)-exact)');
    fprintf('\n')
end
% END FILE: fxptite2.m
% FILE: hlpd034.m
