How to introduce two successive points inside the FixedPointList for each cycle?

Compiling to get some speedup:

fc = Compile[{{t, _Complex, 1}, {n, _Integer, 0}, {b, _Real, 0}},
  Module[{x, y, resp},
   {x, y} = t;
   resp = {x + ((-1 + x^n) (x + x b - 2 y b))/((y b)^
          n - (x + x b - y b)^n), x}],
  {{x, _Complex, 0}, {y, _Complex, 0}, {resp, _Complex, 0}}, 
  "RuntimeOptions" -> "Speed", CompilationTarget -> "C"]

b = 0.8; n = 3;
Rasterize[DensityPlot[Length[
    FixedPointList[
     fc[#, n, b] &, {# + $MachineEpsilon (1 + I), #} &@(x + I y), 15, 
     SameTest -> (Abs[#1[[1]] - #2[[1]]] <= 10^-2 &)]], {x, -2., 
    2.}, {y, -2., 2.}, Mesh -> False, ColorFunction -> (Hue[#] &), 
   PlotPoints -> 150], ImageResolution -> 72, 
  ImageSize -> 350] // AbsoluteTiming

Mathematica graphics


Code and picture first:

With[{b = 0.8, n = 3},
     f[x_] := x^n - 1; 
     DensityPlot[Length[FixedPointList[{Last[#], 
                        With[{w = #2 + b (#2 - #1)},
                             #2 - f[#2]/Sum[(b #1)^(n - k - 1) w^k, {k, 0, n - 1}]] & @@ #} &,
                        {x + I y + $MachineEpsilon (1 + I), x + I y}, 14, 
                        SameTest -> (Norm[Differences[#2]] <= 1.*^-2 &)]] + 1,
                 {x, -2, 2}, {y, -2, 2}, ColorFunction -> Hue, Mesh -> False,
                 PlotPoints -> 50, PlotRange -> All]]

loopy looking divergence

Notes:

  1. The divided difference of $p(x)=x^n-1$ can be re-expressed in a way that is stable to subtractive cancellation:

$$\frac{p(b)-p(a)}{b-a}=\sum_{j=0}^{n-1}a^{n-j-1}b^j$$

For more general polynomials, Velvel Kahan has derived a numerically stable algorithm based on Horner's method to evaluate the divided difference:

polynomialDividedDifference[poly_, {x_, a_, b_}] /; PolynomialQ[poly, x] :=
          Module[{d = 0, y = 0},
                 Do[y = b y + Coefficient[poly, x, k]; d = a d + y,
                    {k, Exponent[poly, x], 1, -1}];
                 Expand[d]]

This would be useful for the general iterative formula displayed in the OP.

  1. To use FixedPointList[] on two starting values, you can group these two values as a list, and modify the iteration function and SameTest setting to handle the initial list.

This strategy can be illustrated using the simpler NestList[]:

NestList[{Last[#], f @@ #} &, {a, b}, 3]
   {{a, b}, {b, f[a, b]}, {f[a, b], f[b, f[a, b]]},
    {f[b, f[a, b]], f[f[a, b], f[b, f[a, b]]]}}

Of course, this method will undercount the iterates, so you need to add a correction term. (As an analogy, compare the results of Length[Range[n]] and Length[Partition[Range[n], 2, 1]].)

As for the SameTest setting, you can feed only the second argument to the test function to check for convergence. (Recall that each "iterate" is actually the last two iterates of the original iteration function.)

  1. The choice of the other starting point where a tiny perturbation of the starting value was used is of course completely arbitrary.