function cv
%CV Convergence rates comparison.
%   CV compares convergence rates for several root finding methods.
%
%   See also NEWTON, BISECTION, SECANT, IQI.

%   Francisco J. Beron-Vera
%   $Revision: 1.1 $  $2006/09/21 11:11:32 $


f=@(x) x^3+x-1;
dfdx=@(x) 3*x^2+1;
d2fdx2=@(x) 6*x;

N=54;

xN=newton(f,dfdx,0,N);
xB=bisection(f,0,2,N);
xS=secant(f,0,2,N);
xI=iqi(f,0,1,2,N);

eN=abs(xN-xN(end));
eB=abs(xB-xB(end));
eS=abs(xS-xS(end));
eI=abs(xI-xI(end));
h=loglog(...
   eN(1:end-2),eN(2:end-1),...
   eB(1:end-2),eB(2:end-1),...
   eS(1:end-2),eS(2:end-1),...
   eI(1:end-2),eI(2:end-1),...
   'LineWidth',2); hold on
x=[...
   min([eN(1:end-2)' eB(1:end-2)' eS(1:end-2)' eI(1:end-2)'])...
   max([eN(1:end-2)' eB(1:end-2)' eS(1:end-2)' eI(1:end-2)'])];
S=abs(.5*d2fdx2(xN(end))/dfdx(xN(end)));
h=[h;loglog(x,.5*x,'-.k',x,S*x.^((sqrt(5)+1)/2),'--k',x,S*x.^2,':k')]; hold off
legend(h,...
   {'Newton';...
   'Bisection';...
   'Secant';...
   'IQI';...
   '$e_n = \frac{1}{2}e_{n-1}$';...
   '$e_n = \left|\frac{f''''(r)}{2f''(r)}\right|e_{n-1}^{\frac{\sqrt{5}+1}{2}}$'
   '$e_n = \left|\frac{f''''(r)}{2f''(r)}\right|e_{n-1}^2$'},...
   'FontSize',14,...
   'Interpreter','LaTeX')
xlabel('$e_{n-1}$',...
   'FontSize',14,...
   'Interpreter','LaTeX')
ylabel('$e_n$',...
   'FontSize',14,...
   'Interpreter','LaTeX')

%--------------------------------------------------------------------------
function x=newton(f,dfdx,a,N)
% Newton method.
x=zeros(N,1);
b=a+1;
n=0;
while abs(b-a)>eps*abs(b)
   n=n+1;
   if n>N, break, end
   b=a;
   a=b-f(b)/dfdx(b);
   x(n)=a;
end
x=x(1:n-1);
%--------------------------------------------------------------------------
function x=bisection(f,a,b,N)
% Bisection method.
x=zeros(N,1);
n=0;
while abs(b-a)>2*eps*abs(max(b,1)) || abs(f(b))>realmin
   n=n+1;
   if n>N, break, end
   r=(a+b)/2;
   if f(a)*f(r)<0
      b=r;
   else
      a=r;
   end
   x(n)=r;
end
x=x(1:n-1);
%--------------------------------------------------------------------------
function x=secant(f,a,b,N)
% Secant method.
x=zeros(N,1);
n=0;
while abs(b-a)>2*eps*abs(max(b,1)) || abs(f(b))>realmin
   n=n+1;
   if n>N, break, end
   c=b;
   b=a;
   a=c-f(c)*(c-b)/(f(c)-f(b));
   x(n)=a;
end
x=x(1:n-1);
%--------------------------------------------------------------------------
function x=iqi(f,a,b,c,N)
% Inverse quadratic interpolation method.
x=zeros(N,1);
n=0;
while abs(b-a)>2*eps*abs(max(b,1)) || abs(f(b))>realmin
   n=n+1;
   if n>N, break, end
   d=c;
   c=b;
   b=a;
   A=f(b)/f(c);
   B=f(d)/f(c);
   C=f(d)/f(b);
   a=d-(B*(B-A)*(d-c)+(1-B)*C*(d-b))/((A-1)*(B-1)*(C-1));
   x(n)=a;
end
x=x(1:n-1)