6
$\begingroup$

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

enter image description here

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.

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – m_goldberg
    Mar 2, 2016 at 15:08
  • $\begingroup$ 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. $\endgroup$
    – Faz
    Mar 2, 2016 at 15:21

2 Answers 2

4
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – Faz
    Mar 2, 2016 at 16:29
  • $\begingroup$ @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 $\endgroup$ Mar 2, 2016 at 16:43
  • $\begingroup$ Note: I'm not using zero anymore as the first point $\endgroup$ Mar 3, 2016 at 2:29
  • $\begingroup$ 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. $\endgroup$
    – Faz
    Mar 3, 2016 at 12:39
2
$\begingroup$

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.
$\endgroup$
4
  • $\begingroup$ 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! $\endgroup$
    – Faz
    Mar 3, 2016 at 12:39
  • $\begingroup$ All the more reason to use a stable evaluation formula, like the one I demonstrated for the divided difference of a polynomial. $\endgroup$ Mar 3, 2016 at 12:43
  • $\begingroup$ 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? $\endgroup$
    – Faz
    Mar 3, 2016 at 14:10
  • $\begingroup$ 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. $\endgroup$ Mar 3, 2016 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.