# Basins of attraction when using Newton's method to solve $x^3-1 = 0$

I created Newton fractal plot for x^3 - 1 and would like to know where did I go wrong. I know I could have just copied the code but that would not have been fun.

x = N[FindDivisions[{-2, 2}, 200], 8];
y = N[FindDivisions[{-2, 2}, 200], 8];
fx = xn^3 - 1;
f1x = 3*xn^2;
r = {};
r2 = {};
r3 = {};

i = 1;
j = 1;
iter = 1;
For[i = 1, i < 200,
xi = x[[i]];
For[j = 1, j < 200,
yi = y[[j]];
xn = xi + I*yi;
For[iter = 1, iter < 7,
xn1 = N[xn - fx/f1x, 8];
xn = xn1;
iter++;
]
Which[Im[xn] < 010, AppendTo[r, {x[[i]], y[[j]]}], Im[xn] > 010,
AppendTo[r2, {x[[i]], y[[j]]}], Im[xn] == 010,
AppendTo[r3, {x[[i]], y[[j]]}]]
j++;
]
i++;]

a = ListPlot[r, PlotStyle -> Red];
b = ListPlot[r2, PlotStyle -> Blue];
c = ListPlot[r3, PlotStyle -> Green];
Show[a, b, c]


As you can see it looks similar but not exactly like the required one. Would like to know what I can do to make it better. Should I increase the number of test cases? Would it help?

• You write: "As you can see it looks similar but not exactly like the required one. ". No, I can't see that as you supply no "required one" to compare it to. – m_goldberg Jun 27 '17 at 1:58
• See: 100053, 15047 – Kuba Jun 27 '17 at 5:40
• Hint. Your use of arbitrary precision arithmetic at low precision is not speeding up your code. Machine arithmetic is always the fastest way to compute. – m_goldberg Jun 27 '17 at 14:28
• More hints: Your For-loop syntax is incorrect -- your iterator increments should be given as the 3rd argument to For. You are a missing the semi-colon needed to separate your Which expression from For-loop that precedes it. – m_goldberg Jun 27 '17 at 15:03
• how to located three roots in the second picture. could you tell me how does it make in matlab? – onk Feb 22 '18 at 9:06

Just to illustrate avoidance of For loop, use of MatrixPlot (requires change of representation):

f[z_] := z - (z^3 - 1)/(3 z^2)
roots = {{1, 0}, {-1/2, Sqrt[3]/2}, {-1/2, -Sqrt[3]/2}};
nt[x_, y_, n_] :=
FixedPoint[f[Complex @@ #] &, {x, y}, n] /. Complex[a_, b_] :> {a, b}
cfn[x_, y_, n_] :=
With[{r = nt[x, y, n]}, Which[Norm[r - {1, 0}] < 10^(-6), 1,
Norm[r - {-1/2., Sqrt[3]/2.}] < 10^(-6), 2,
Norm[r - {-1/2., -Sqrt[3]/2.}] < 10^(-6), 3]]
Legended[MatrixPlot[
Reverse@Transpose@
Table[cfn[i, j, 40], {i, -2, 2, 0.015}, {j, -2, 2, 0.015}],
DataRange -> {{-2, 2}, {-2, 2}},
ColorRules -> {1 -> Red, 2 -> Green, 3 -> Blue},
Epilog -> {Yellow, PointSize[0.02], Point[roots]}],
SwatchLegend[{Red, Green, Blue}, Style[#1 + #2 I, 20] & @@@ roots]]


You need a lot more iterations at a lot more points to get a good image of the basins of attraction. With a faster implementation than yours, I was able to get a reasonably good image when iterating 15 times at 250,000 points (500 x 500 lattice).

I also think your code may have bugs that are distorting the image, so I'm not sure you will get a correct image even if you increase the number of iterations and use more points.

This is the image my code produces.

### Update

I have two sets of code that can generate the above image. The 1st is a repaired and speeded up version of the OP's code. I don't recommend using this code, but it might be instructive to the OP and others to see what is more or less the minimum needed to make the posted code work.

fx[u_] := u^3 - 1;
f1x[u_] := 3 u^2;

n = 500;
y = x = Subdivide[-2., 2., n];
max = 15;
r1 = r2 = r3 = r[];

For[i = 1, i <= n + 1, i++,
xi = x[[i]];
For[j = 1, j <= n + 1, j++,
xn = xi + I*y[[j]];
Quiet @
For[iter = 1, iter < max, iter++,
xn1 = xn - fx[xn]/f1x[xn];
xn = xn1];
imxn = Im[xn];
Which[
imxn == 0., r3 = r[{x[[i]], y[[j]]}, r3],
imxn < 0., r1 = r[{x[[i]], y[[j]]}, r1],
imxn > 0., r2 = r[{x[[i]], y[[j]]}, r2]]]]


Accumulating the deeply nested r expressions is much faster than the AppendTos the OP was using to accumulate basin points. This is a major speed boost.

We must flatten the results afterward, but that too is fast.

pts = List @@@ (Flatten @* r /@ {r1, r2, r3});

a = ListPlot[pts[[1]], PlotStyle -> Blue];
b = ListPlot[pts[[2]], PlotStyle -> Yellow];
c = ListPlot[pts[[3]], PlotStyle -> Red];
Show[a, b, c,
Axes -> None,
AspectRatio -> Automatic,
ImageSize -> Large]


The 2nd code set is the one I recommend. This, in essence, remains a straight-forward translation of OP's code, preserving the the basic procedural approach, but replacing elements of what I consider bad practice with code that is widely accepted as good practice. Note that I do not introduce any global variables. Only the function definitions are made at top-level.

fx[u_] := u^3 - 1;
f1x[u_] := 3 u^2;

With[{n = 500, max = 15},
pts =
Module[{x, b, b1, b2, b3, root, yRoot},
x = Subdivide[-2., 2., n];
b3 = b2 = b1 = b[];
Do[
root = Quiet @ Nest[# - f[#]/df[#] &, x[[i]] + I x[[j]], max];
yRoot = Im[root];
Which[
yRoot == 0., b1 = b[{x[[i]], x[[j]]}, b1],
yRoot < 0., b2 = b[{x[[i]], x[[j]]}, b2],
yRoot > 0., b3 = b[{x[[i]], x[[j]]}, b3]],
{i, n + 1}, {j, n + 1}];
List @@@ (Flatten @* b /@ {b1, b2, b3})]];


Do and Nest are faster than For, and only one Do expression is needed. Further, since the point lattice being constructed is the same in both directions, only one list of divisions is needed.

ListPlot[{pts[[1]], pts[[2]], pts[[3]]},
PlotStyle -> {Red, Blue, Yellow},
Axes -> None,
AspectRatio -> Automatic,
ImageSize -> Large]


Show` is not needed.

• Would you mind sharing it? – user49702 Jun 27 '17 at 9:30
• @user49702. What are you asking me to share? Your "it" is a pronoun with no reference – m_goldberg Jun 27 '17 at 14:31
• Perhaps, the OP is referring to the code you used to produce your plot. – bbgodfrey Jun 27 '17 at 15:13
• @bbgodfrey. Probably. But it is not perfectly clear and I think it should be. I don't want to make an effort only to be told what I posted wasn't what was asked for. (Happened to me more than once already) – m_goldberg Jun 27 '17 at 15:32
• @user49702. I am willing to post my code if you will confirm that is what you are asking me to do. – m_goldberg Jun 27 '17 at 15:40