Halley's method

Question

Halley's method uses a quadratic Taylor approximation and results in a fixed point method of order $3$:

$$ x_{n + 1} = x_{n} - {{\rm f}\left(x_{n}\right) \over {\rm f}'\left(x_{n}\right)}\, \left[% 1 - {{\rm f}\left(x_{n}\right){\rm f}''\left(x_{n}\right) \over 2{\rm f}'^{2}\left(x_{n}\right)} \right]^{-1} $$

My original question about finding the cube root of 5 using Halley's method has been solved.

How do I verify numerically that the convergence is cubic?

I know that I have to use the order of convergence formula, but how do I set it up on Maple?

Answer 1

You can estimate the order of the convergence numerically. The procedure estimateorder below is just one of several ways to do it.

First, H is a short and workable (but inefficient) implementation of the method, with no stopping criteria other then a maximal number of iterations. Note that this procedure stores and returns all the iterates, not just the final one. Presumably you already have your own implementation.

restart:

H:=proc(expr,x,x0,maxiter)
  local ddf, df, f, n, X;
  f:=unapply(expr,x);
  df:=D(f);
  ddf:=(D@@2)(f);
  X:=Vector[':-row'](maxiter,[x0]);
  for n from 1 to maxiter-1 do
    X[n+1]:=evalf(X[n]-f(X[n])/df(X[n])/(1-(f(X[n])*ddf(X[n]))/(2*df(X[n])^2)));
  end do;
  X;
end proc:

Let's run that on an example, and make sure to have lots of digits in the results, for comparing.

Digits := 100:
V := H( x^3-5, x, 37.0, 10 ):

Here's the last term in that result,

V[10];

        1.709975946676696989353108872543860109868055110543054924382861707444295920504\

           173216257187010020189002

Now a procedure which estimates the order of convergence.

estimateorder:=proc(S::Vector)
  local d1, d2, dim, i, k, n1, n2, P;
  dim := op(1,S);
  P := Vector[row](dim-3);
  k := dim-3;
  for i from 1 to dim-3 do
     d1, d2 := S[i+2]-S[i+1], S[i+1]-S[i];
     if d1=0.0 or d2=0.0 then
        P[i]:=FAIL; k:=i-1; i:=dim-3; next;
     end if;
     n1, n2 := S[i+3]-S[i+2], S[i+2]-S[i+1];
     P[i] := ln(abs(n1/d1))/ln(abs(n2/d2));
     if not type(P[i],numeric) then
        P[i]:=FAIL; k:=i-1; i:=dim-3; next;
     end if;
     P[i] := evalf[5](P[i]);
  end do;
  P[1..k];
end proc:

estimateorder(V);

          [1.0089, 1.0692, 1.4427, 2.3523, 2.8416, 2.9932, 3.0000]

Notice that the convergence rate may be much worse initially, and that it only tends to cubic.

Another example, for illustration

Digits := 100:
ee:=randpoly(x,degree=7,dense);

             7       6       5       4       3       2     
         -7 x  + 22 x  - 55 x  - 94 x  + 87 x  - 56 x  - 62

V:=H( ee, x, 37.0, 30 ):

V[20]; # the final term

   -1.6254993522230776878001145071381988387744828918946894060940899\

      42338528314967927805393172797627178796

eval(ee,x=%): evalf[5](%); # is it close to being a root?

                                 -71
                        9.7076 10   

convert( estimateorder(V), list);

    [0.99752, 0.99468, 0.98841, 0.97453, 0.94513, 0.89261, 0.85937, 

     1.0931, 1.2070, -2.0048, -2.8923, -1.2156, -0.85385, 0.95828, 

     2.7835, 2.9570, 2.9997, 3.0000]

estimateorder(V)[-1];

                          3.0000