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

If we want to draw the attraction basins of an iteration function of the following type $$x_{k+1}=x_k-\frac{f(x_k)}{\frac{f(x_k)-f(w_k)}{x_k-w_k}},$$ where $w_k=x_k+b f(x_k)$, $b\in R-\{0\}$, we can simply apply FixedPointList to find the fixed points and use DensityPlot to draw the final picture as follows:

ClearAll["Global*"]
f[z_] := z^n - 1
Simplify[z - f[z]/((f[z] - f[z + b*f[z]])/(z - (z + b*f[z])))]
b = 0.8; n = 3;
Rasterize[
DensityPlot[
Length[FixedPointList[# + ((-1 + #^n)^2 b)/(#^
n - (# + (-1 + #^n) b)^n) &, x + I y, 15,
SameTest -> (Abs[#1 - #2] <= 10^-2 &)]], {x, -3., 3.}, {y, -3.,
3.},
Mesh -> False, ColorFunction -> (Hue[#] &), PlotPoints -> 250],
ImageResolution -> 72, ImageSize -> 350] // AbsoluteTiming


Here the fixed-point formula could be written as $$\frac{b \left({x_k}^n-1\right)^2}{{x_k}^n-\left(b \left({x_k}^n-1\right)+{x_k}\right)^n}+{x_k}.$$ The above mentioned iterative method needs/receives only one effective point per cycle to gives the next fixed point. But, my inquiry is for the case when two successive points are involved in the iterations. As an illustration for the following Secant-like method

$$x_{k+1}=x_k-\frac{f(x_k)}{\frac{f(w_k)-f(b x_{k-1})}{w_k-b x_{k-1}}},$$

wherein $w_k=x_k+b(x_k-x_{k-1})$. For this case, the fixed-point formula could be written as $$\frac{\left({x_k}^n-1\right) (b {x_k}+ {x_k}- 2 b {x_{k-1}})}{(b {x_{k-1}})^n-(b {x_k}+ {x_k}-b {x_{k-1}})^n}+{x_k}.$$

But the following code does not work:

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


Here, we have to introduce two successive points into the FixedPointList to finally draw the attraction basins, which I failed to do? The interesting point is that #1 and #2 are two automatically saved points in the (Abs[#1 - #2] <= 10^-2 &), but I cannot call them in the fixed-point formula!

And finally in general, how can we deal with iterations which require even more than 2 successive points inside the FixedPointList? I would be very thankful if someone could provide some tips or answers for such cases.

• I suggest you look at NestWhileList, which is more general than FixedPointList. It's 4th argument specifies how many previous values to use in the completion test. Mar 2, 2016 at 15:08
• I tried to use it but I failed again. If you believe that it helps to finally draw the attraction basins with two seeds, please provide it in a complete answer.
– Faz
Mar 2, 2016 at 15:21

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  • Thanks. But I have a query. Why when we choose$b=0.1$, the code seems to break down and does not provide an attraction basin? How did you choose the two initial approximations ($x_{-1},x_0$) for converging to each fixed point? – Faz Mar 2, 2016 at 16:29 • @FazlollahSoleymani Apparently the thing doesn't converge for b < .37 or something near it. I used a fixed$x_{-1} = 0$, you may use whatever you want Mar 2, 2016 at 16:43 • Note: I'm not using zero anymore as the first point Mar 3, 2016 at 2:29 • I got confused! It seems to me that your first response and also the response of JM (below) are better in the sense of final output fractal. There is no need to calculate the second point by another algorithm. Please choose$x_0$as usual and$x_{-1}=x_0+0.001$whenever the iteration is started. After starting, we have two values at each cycle (old one and the new approximation), clearly the new approximation replaces the old one but we keep both. Maybe, in this way, we could obtain a same final fractal for your approaches. – Faz Mar 3, 2016 at 12:39 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]]


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.
• I cast doubt that in theoretic $b$ can be a small value, e.g $b=0.1$, but here we cannot draw its attraction basin! Why? Maybe, due to round-off errors!
– Faz
Mar 3, 2016 at 12:39
• All the more reason to use a stable evaluation formula, like the one I demonstrated for the divided difference of a polynomial. Mar 3, 2016 at 12:43
• What if we use the computations in multiple precision arithmetic, e.g. by applying 32 digits to avoid such cancellation and round off errors so as to draw the attraction basin for small $b$ as well?
– Faz
Mar 3, 2016 at 14:10
• DensityPlot[] supports a WorkingPrecision option, sure. Why not try it out yourself? The only other changes needed in my code would be to use an exact value for $b$ and to replace \$MachineEpsilon` with a high-precision tiny number. Mar 3, 2016 at 14:12