- Email: [email protected]
Newton raphson method, scaling at fractal boundaries and mathematica
Mathl. Comput. Modelling Vol. 21, No. 10, pp. 91-102, 1995
Pergamon
Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/95 $9.50 + 0.00
08957177(95)00073-9
Newton Raphson Method, Scaling at F’ractal Boundaries and MATHEMATICA A. K. BISOI” International Centre for Theoretical Physics 34100 Trieste, Italy (Received
September
1994; accepted
November
1994)
Abstract-The basins of convergence of cubic polynomials having real roots are studied using the Newton Haphson iterative method. The limiting value for the ratio of basin segments for equispaced roots is explained. An algorithm is presented for computing the basin boundaries on the real axis which obviates the necessity of taking recourse to extensive search. MATHEMATICA programs have been developed to help in the above study and also to depict the basins in the complex plane. Keywords-Fmctal basin boundaries, plex dynamics, Julia set.
Newton Taphson,
MATHEMATICA
programming,
Com-
1. INTRODUCTION Fractals are encountered in feedback processes, where iterations are involved, and are characterized by self-similarity. Fractal objects contain their replicas embedded in them, although reduced by a scale factor. If the scale factor is constant over replications, the fractal is strictly self-similar. However, if the scale factor is only asymptotically constant, we say the fractal is quasi-selfsimilar [I]. This scale factor is conceptually related to the fractal dimension. The fractal boundaries between the basins of convergence of roots of polynomials obtained by applying the Newton Raphson (NR) iterative method have been the subject of many studies [2-61. The use of computer graphics for depiction of these boundaries have produced many pictures of awe-inspiring complexity where the phenomenon of self-similarity is visually apparent in qualitative terms. Self-similarity involves repetition of a structure with a change of scale at each repetition. One wonders what is the scale factor involved in the NR method. What is the fractal dimension of the basin boundary? Must one always follow the time-consuming procedure of searching pixel by pixel to locate the basin boundaries? Guided by these questions, we study a family of cubic equations, having real roots and try to explain the scaling involved at the fractal boundaries. Such a study, involving computer experiments, has already been initiated by Becker and DSrfler [2]. We first summarize their result. Becker and Dorfler have studied the cubic polynomial f(X) = (Lr+ 1)X(2 - l), which has zeros at 2 = -1,
0 and 1.
(I)
The plot of the function is given in Figure 1.
The
*Permanent Address: Department of Computer Science and Application, Utkal University, Vani Vihar, Bhubaneswar 751004, India. The author would like to thank the International Centre for Theoretical Physics, Trieste for hospitality and A. Nobile and the members of Scientific Computing Section of ICTP for extending all possible help.
Typeset by 44-W 91 MCM 21/10-G
A. I<. BISOI
92 0.4
-0. I
-0.2 -0.3 -0.4
-0.5
-I
Figure 1. Plot of equation (1). number of alternating basins.
Intervals (-d,
-e)
and (d,e)
derivative f’(x) is zero at x = -d and x = d, where d = l/a. are presented with three broad regions of x; namely (i) x < -d
(ii) - d < x < d
I
0.5
0
and
contain
an infinite
Thus, in the first instance, we
(iii) x > d,
(2)
which are intuitively the most likely regions representing the basins of convergence of roots -1, 0 and 1, respectively. However, actual computation reveals that the basin of convergence of root at x = 0 is not the region from -d to d but a smaller region from -e to e. The point x = e is such that the Newton map
N(x) = returns this point after two iterations. This leads to
$$
Each iteration only changes the sign, i.e., N(e) = -e. e = 0.5%. ‘e
For the cubic under consideration, e = l/d. With 15 digit precision, Becker and Diirfler have identified 12 basin boundaries, gr to gr2, in the x region from -d to -e. These boundaries have been found through extensive search and computation in this region. We shall discuss in Section 2 how extensive search, which is very time-consuming, can be avoided. Between gr and grz, the basins of roots at x = -1 and x = 1 alternate with the width of the basins becoming narrower and narrower. The structure on the positive x axis is symmetric. Becker and Dorfler find that Sn - Sk-1 = 6.0. n+cCJQn+l - gn lim
(5)
They also remark that the above limit is always 6 for a cubic function having equispaced roots on the real axis. An explanation of these facts will be provided in Sections 2 and 3. Generalization to the case of non-equispaced roots is attempted in Section 3. Depiction of basin Boundaries in the complex plane has been dealt with in Section 4. Some concluding remarks are made in Section 5.
93
NewtonRaphsonMethod
0.4
0.2
0
-0.2
-0.4
-0.452
0.4
=x0
0.2
0
-0.2
-0.4
-I
-0.5
0
0.5
I
Figure2(b). The first two tangentsfor g2 < x0 < g3.
2. FINDING THE BASIN BOUNDARIES WITHOUT DOING EXTENSIVE SEARCH Let us consider the cubic equation with roots at -1,0 and 1 again. Figures 2a and 2b sketch the behaviour of the first two tangents obtained through NR iteration starting from x values which lie in the range gi < x0 < g2 and g2 < xc < gs, respectively. These diagrams have been produced through a program in MATHEMATICA [i’] given in Appendix A. Assigning a larger value to the variable iter in the program would result in the depiction of successive tangents converging to the root. This program has been inspired by a program [8,9] to depict the iterations of quadratic map. From -oo to gi = -d, it is the basin of attraction of the root at x = -1. Note that the tangent at f(gi) runs parallel to the x-axis and does not lead to any root. However, it being impossible to represent the irrational number fi in the computer, there is no danger of landing exactly on this starting value. As we move x0 to the right of gi, the point ti, where the tangent intercepts the x axis, moves more and more towards the left. We are in the basin
A. K. BISOI
94
of convergence of the root at x = 1 as long as the point ti lies to the right of -91 = d. As ti moves to the left of x = d, we enter the basin of convergence of the root at x = -1. Therefore, the point gs is such that it is mapped to x = l/t/5 m . a single NR iteration. Hence, the value of g2 is obtained by solving the cubic equation - 39; + 1 = 0.
2&g;
The value for the real root of equation (6) obtained using the MATHEMATICA
function NSolve
agrees with the g2 value of Becker and Dorfler. As we move x0 from g2 to gs, the point ti moves from -gi to -g2 on the positive z axis. This results in the point t2, the point of intersection of the tangent of the second iteration with the x axis, moving from -oo
to gi, (which is the
exclusive basin of the root at x = -1). From ga to g4, it is the basin of the root at x = 1 again. This structure of alternating basins continues an infinite number of times. The value of gs is obtained by solving the cubic equation 29; + 3g2g; - g2 = 0.
(7)
Thus, we have obtained the following recursive algorithm for finding the g’s. Starting with g1=
-5,
(8)
the successive g’s are obtained by finding the real root of the cubic 29; + 3g+1 g” - gi-1 = 0,
(9)
for i = 2,3,. . . , lc. Since the difference between the successive values of the g’s goes on decreasing, at some point i = Ic, it will no longer be possible to distinguish between g&i and gk. This leads to the limiting value
2g; + 39; -
0
gk =
(10)
or gk = 55
= +e.
We are interested in the limiting value of the ratio (gi -gi-i)/(gi+i -gi). of the differences can be replaced by the differential coefficient dgi_r/dgi.
s-1 Calculating dgi_i/dgi tion (II), we obtain
= -3gp
2gs _
1’
(11) For large i, this ratio From equation (9),
(12)
from equation (12) and substituting the limiting value of gi from equaQ-1 .limmdg,=
6
.
(13)
This explains the ratio obtained by Becker and Dijrfler. We can also understand why the basins appear as alternating segments. A MATHEMATICA program implementing the above algorithm for computing the g’s is given in Appendix B together with the output from this program. It has been chosen to calculate 30 number of g’s with a precision of 50 digits. It is to be noted that there exist an infinite number of g’s in the x interval from -d to -e and also in the interval from d to e. The higher the precision we choose, the more is the number of g’s we can identify. It is easy to perform calculations with arbitrary precision using MATHEMATICA. The ratios defined by Becker and DGrfler are also calculated with a precision of 25 digits.
95
Newton Flaphson Method
3. NONEQUISPACED
ROOTS
Consider a cubic equation with distinct real roots at CY,/3 and y such that (I: < p < y. The mapping x’ = (z - /3)/(p - a) relocates the roots at -l,O,l in the transformed variable z’, if the roots Q, p, y are equispaced, i.e., if (y - ,D) = (,0 - CX). Thus, we can understand why the limit given in equation (5) is the same for all cubic functions having equispaced roots on the real axis. The above transformation the mapped variable.
maps the nonequispaced roots to -1,0
and c = (y - @)I(/? - CX)in
Thus, it will suffice to study a cubic of the form f(z)
= (X + 1)x(x - c) to
generalize to the case of nonequispaced roots. Because of the asymmetrical
location of the roots,
the basin boundaries on the left and right hand side will no longer be symmetrically
situated.
Let us denote the basin boundaries on the left and the right hand sides by li and ri with i = 1,2,3,. . . . These points belong to the Julia set, the members of which do not converge to any root. The li and ri can be recursively generated starting from 11 and ~1 given by the solution of f’(z)
= 0; i.e.,
1
1-
=A-B 3
and
= -A+B 3
7.1
where A = c - 1 and B = dm.
’
The real root of the cubic equation N(l)(li) = Ti-1 determines li and similarly N(l)(ri) = li_l determines ri for i = 2,3, . . . , where iV(‘) represents a single iteration of the Newton map N(x) = Calculating
2x3 + (1 - c)s2 3x2 + 2(1-
(14)
c)x -c’
the succeeding l’s and the T’S involves finding the real root of the cubic equation 2x3 - (A + 3w)x2 + 2wAx + WC= 0,
(15)
which, for w = ri_1, returns x = li as the real root and for w = li_1, returns x = ri as the real root. A MATHEMATICA program implementing the above algorithm is given in Appendix C. The program computes the location of basin boundaries for various values of c. The limiting value of li would be such that lk = N(l)(?-k) = Nc2)(lk). This is an equation with a leading power of 6 in the unknown variable, and hence, an analytic solution is not possible in the general case. The ratios defined by (I) _ fj
-
4+1 - 4 &2 _ 1,,1;
p!” 3
=
rj+l rJf2
-
rj rj+l
(16)
are calculated for j = 1 to k - 3 and are given in Appendix C for various values of c. It is noticed that, for the case of nonequispaced roots, p’s with even and odd subscripts tend to separate limiting values. We conjecture that each of the basins of convergence would be scaling down by a constant factor. To verify, we define the ratios as JO _ 3
-1.
lj
.
3+3 - $+2
&+1
-
’
(17)
for j = 1 to k - 3 and observe from Appendix C that the values for sj are suggestive of an asymptotically constant value, s,(c) for each value of c. The variation of so0 with c can be fitted by a quadratic polynomial. It is to be noted that, although the initial values of p and s differ on the left and right hand sides, the asymptotic limits are the same. From the values of the ratio s, it is apparent that the basin segments of roots at -1 and 1 scale down by a factor of the same magnitude.
96
A. K. Brsor
4. DEPICTION
OF BASINS
IN THE COMPLEX
PLANE
To reproduce the depiction of basins in the complex plane using the Pascal programs of Becker and Diirfler, one also has to be aware of machine specific graphic primitives which must be interfaced to these programs. We find that MATHEMATICA
provides a short-cut and a ready
made answer to produce such plots through use of the DensityPlot DensityPlot
function.
A program to
the Mandelbrot set is given in the book by Crandall [lo]. We give in Appendix D
an adaptation of this program to plot the basins pertaining to the NR method in the complex plane.
0.5
0
-0.5
-1 -0.5
0
0.5
1
Figure 3. Basins in complex plane for roots at -1,0
1.5
and 2.
Figure 3 depicts the basins in the complex plane for the cubic with roots at -1, 0 and 2. The white region belongs to the basin of root at 2, the dark region to the root at -1 and the gray region to the root at 0. The picture is symmetric about the x axis. Therefore, in Figures 4 and 5, we plot the diagram for y > 0 only. Figures 4 and 5 depict enlargement of the left and right boundaries of Figure 3. For producing each of these figures, a grid of (200 x 200) points have been plotted. One can improve the resolution of the pictures by choosing a larger number of grid points. Consider the equation N(r)(z) x1 = -0.47664;
= ~1 again. The three roots of this cubic equation are
x2 = 1.39975 + 0.76830 i;
x3 = 1.39975 - 0.76830 i.
We have identified the real root x1 with 12, a point belonging to the Julia set on the negative x-axis. It turns out that 52 and x3 also belong to the Julia set. The location of the point 52 has been marked with a cross (+) in Figure 5 at the germination point of a bud on the right boundary.
Newton Raphson
n
I,,
-0.8
-0.9
Method
”
”
-0.6
-0.7
Figure 4. Upper half of left boundary
0.8
...
8,.
1
97
”
c -0.5
”
”
-0
L
of Figure 3.
.
*
1.2
1.4
Figure 5. Upper half of right boundary
”
of Figure 3.
*I
1.6
A. K. BISOI
98
5. CONCLUDING We have seen that it is not always necessary Analysis
of behaviour
of NR tangents
REMARKS
to make a step by step search for basin boundaries.
leads to the development
of an algorithm
for computation
of the basin boundaries on the real axis. This deterministic algorithm demystifies the occurrence of an infinite number of basin segments on the real axis. The behaviour of the points near the boundary may still appear to be chaotic. However, that is due solely to the finite precision of the computers and not due to any inherent chaos in the dynamical system itself. The execution of the algorithm In computing is felt that determination
requires
much less time compared
these boundaries,
exploring
the iteration
of boundary
points
to the time required
by a search method.
we have made use of only the real root of equation of the complex
roots of this cubic equation
in the complex
plane,
might
as we have demonstrated
(15).
It
lead to the
for one point
on Figure 5. The relation of s, to fractal dimension needs to be explored. It is hoped that this work will stimulate further application of MATHEMATICA to study the dynamics of NR method and search for deterministic algorithms to compute basin boundaries in the complex plane.
REFERENCES 1. S.D. Casey and N.F. Reingold, Self-similar fractal sets: Theory and procedure, IEEE Computer Graphics and Applications 14 (3), 73 (1994). 2. K.H. Becker and M. Dorfler, Dynamical Systems and Fractals, Cambridge University Press, (1989). Computers Math. 3. D.C. Arney and B.T. Robinson, Exhibiting chaos and fractals with a microcomputer, Applic. 19 (3), 1-12 (1990). 4. H.E. Benzinger, S.A. Burns and J.I. Palmore, Chaotic complex dynamics and Newton’s method, Physics Lett. A 119 (9), 441 (1987). 5. H.O. Peitgen, D. Saupe and F.v. Haeseler, Cayley’s problem and Julia sets, The Mathematical Intelligencer 6 (2), 11 (1984). 6. H.O. Peitgen, Editor, Newton’s method and dynamical systems, Acta Applicandae Muthematicae 13 (1988). 7. S. Wolfram, Mathematics, 2 nd edition, Addison-Wesley, (1991). 8. S. Wagon, Mathematics in Action, W.H. Freeman and Company, (1991). 9. A. Shiekh, A Mathematics essay on chaos, Lecture notes, College on Computational Physics, ICTP, Trieste, (1993). 10. R.E. Crandall, Mathematics
for the Sciences, p. 160, Addision-Wesley,
(1991).
Newton
Raphson
APPENDIX (*Program
to produce
Figure
2(a);
Figure
2(b)
A is obtained
tangents[f_,g_,x0_,displaylength_,xrange_:{-l.l,l.l}]:= Block[{z,p,plot, lplot, cplot,start, list, xmin,xmax}, {xmin, xmax} =xrange; gl=-0.57735; g2=-0.46560; eps=.Ol; z=N[xO]; start =z; Do[ { dnl={z, O),pb+ll={z,f[4), {n,l,displaylength,2} list =Table [ p[i], {i,displaylength}
zz[nl=z,z=g[zl1, 1;
lplot =ListPlot[list, PlotStyle -> Thickness[.OOOl], PlotJoined -> True, DisplayFunction -> Identity 1; cplot =Plot[ {y=f[x]}, {x,xmin,xmax}, PlotLabel->FontForm[“Fig. 2(a): The first two tangents for gl Thickness[.OOOl], AxesLabel ->{x,“ “}, DisplayFunction -> Identity 1; th=.02; Show [ {cplot, lplot}, GraphicsI {PointSize[.Ol], Point[{start,O}], Text[“gl”,{gl,-eps},{O,l}], Text[“g2”,{g2,-eps},{O,l}], Text[“-gl”,{-gl,eps},{O,-l}], Text[“-g2”,{-g2,eps},{O,-1}], Text[“tl”,{zz[3],eps},{O,-l}], Text[“t2”,{zz[5],eps},{O,-l}], Text[“=xO”start,{.3,.4},{0,-l}] }], Graphics[{AbsoluteThickness[th],{Dashing[{0.01,0.01}],
Line[{{gl,O},{gl,f[gll}}l,
Line[({g2,O},(g2,f[g21))1, Line[{~-gl,O},~-gl,f[-g~l~~l, Li~[i{-g2.o},ip2,r[821}}1 ) 1 DeiaultFont->{ “Times-Italic”,lO}, DisplayFunction -> $DisplayFunction, PlotRange -> {xrange,{-.5,.5}}, Frame->True 1
1 x0=-0.49 ; iter=5; f[x_]=x ‘3-x (* Define cubic function *) g[x_]=2*xA3/(3*x”2 - 1) (* Define Newton map *) tangents[f,g,xO,iter,{-l.l,l.l}]
99
Method
by changing
x0 value *)
100
A. K. BISOI
APPENDIX
B
PROGRAM (* Calculate the basin boundaries without doing extensive search *) f[x_,w_]:=2*x”3+3*w*xA2-w; (*Define Newton map *) digit=50; k=30; sout=OpenWrite”math/cubic.out”] g[l]=N[-l/Sqrt[3],digit]; d=N[-l/Sqrt[5],digit]; d=“; i=l; space= “ ” ; textd” While[ i
, {j,k-3)
] ;
] ;
OUTPUT d=
-0.44721359549995793928183473374625524708812367192231
gj
pj
-0.57735026918962576450914878050195745564760175127013 7.256874166975170204865840 -0.46560062143367758369557882264206914551351750716552 6.179501149801548318846121 -0.45020147778247548915116189993227676825613611355286 6.029219709826801810191761 6.004851109839777676226105 -0.44770950581291017390380321914867329621880782129811 -0.44729618997943611050770656699237291959172375596111 6.000807997279660568129623 -0.44722735965776555512399767546695435157538351107654 6.000134651751163967293937 -0.44721588948213224738364426683460861182242551262533 6.000022441556852205031045 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
-0.44721397782909460094975366836322578275040527333734 6.000003740248317957732033 -0.44721365922144666848487585782979526457672255774767 6.000000623374409732649248 -0.44721360612020511504562984246847538285813648717433 6.000000103895726346363949 -0.44721359726999910897107968141990472300155464848525 6.000000017315954151919637 -0.44721359579496480016694813774716728620680520843668 6.000000002885992352010467 -0.44721359554912574940908250241415201561224632891307 6.000000000480998725150555 -0.44721359550815257430247960819118973631571735920319 6.000000000080166454186634 -0.44721359550132371178525990484819787143881297117573 6.000000000013361075697630 -0.44721359550018556803240516111586939946974676741518 6.000000000002226845949601 -0.44721359549999587740692979290561567200083082781657 6.000000000000371140991600 -0.44721359549996426230268390993756785678246343518122 6.000000000000061856831933 -0.44721359549995899311864292976882825256150182045676 6.000000000000010309471989 -0.44721359549995811492130276641642540273714007478203 6.000000000000001718245331 -0.44721359549995796855507940585794308751305566624677 6.000000000000000286374222 6.000000000000000047729037 -0.44721359549995794416070884576486968756126315043518 -0.44721359549995794009498041908269098162260018165942 6.000000000000000007954839 -0.44721359549995793941735901463566120268985919937807 6.000000000000000001325807 6.000000000000000000220968 -0.44721359549995793930442211389448957301746818878992 6.000000000000000000036828 -0.44721359549995793928559929710429430140956225609353 -0.44721359549995793928246216097259508947502680163677 6.000000000000000000006138 -0.44721359549995793928193930495064522081927410184747 [ Note:Output has been edited -0.44721359549995793928185216228032024270998207436278 to add symbolsand to format] -0.44721359549995793928183763850193274635843340559163
101
Newton Raphson Method
APPENDIX PROGRAM
(Generalization
C
of program of Appendix
B)
f[x_,w_]:=2xA3-(a+3w)xA2+2w*a*x+w*gama (*Define Newton map*) sout=OpenWrite[ “math/appcl4”]; k=lO; digit=30; space=” “; For[ gama=2,gama<4,gama++, Print[ “gama=” ,gama]; Write[sout,“gama=“,gama]; b=N[Sqrt[l+gama+gama”2],digit]; a=N[gama-l,digit]; r[l]=N[(a+b)/3,digit]; i=l; l[l]=N[(a-b)/S,digit]; While[ i
Write[sout,j,space,l[j],space,r[j]]];
Do[ (* calculate the ratios *) ul=N[l[m+l]-l[m],digit]; rhol[m]=N[ul/(l[m+2]-l[m+l])]; sl[m]=N[ul/(l[m+S]-l[m+2])];
ur=N[r[m+l]-r[m],digit]; rhor[m]=N[ur/(r[m+2]-r[m+l])]; sr[m]=N[ur/(r[m+3]-r[m+2])],
Do[ P rint[ rho l[nI,sPace,rhor[n],space,sl[n],space,sr[n]]; Write[sout,rhol[n],space,rhor[n],space,sl[n],space,sr[n]],
{m,k-3}];
{n,k-3}];
a=.; b=.; 1; (* end For *) Close[sout] OUTPUT gama=2
li
5 6 8 9 10
ri
-0.548583770354863530167205251213 1.215250437021530196833871917880 -0.476640869333724045703691984163 0.92050686281515245538451861863 -0.460538069165090573216401849777 0.90106639772494915960227270705 -0.459209517915383575472653379100 0.89629425599320413258572595171 -0.458877018093209225118648522210.89589288065135084229024118365 -0.458848935287665317313695075310.89579224077708478008583382472 -0.458841891016546856508041845800.89578373734252614997401883867 -0.458841295767567234288724171460.8957816042631496402834120321 -0.458841146448579826619728368140.8957814240139575015944589047 -0.458841133830822855086124200040.8957813787981735177199975455
&L’) ?
p3
s?) ?
SW 3
4.467726126371272 15.16134376614832 54.15139i13369052 61.76337392615517 12.12057131569808 4.073740090509531 48.42950009215275 48.43462730032162 3.995644993188464 11.88947410100077 47.30835199602091 47.41800172691879 11.83997879608157 3.988233737178288 47.20145158850765 47.20155592259313 3.986616226384538 11.8352030079332 47.17825062334561 47.18055754246566 11.83415908235831 3.986459506511234 47.17599041325278 47.17599262296456 3.986425236042335 11.83405790173216 47.17549885967337 47.17554775275868
A. K. BISOI
102
gama=3 i 1 2
li ri -0.535183758487996431039740422490 1.868517091821329764373073755823 -0.482748900024651977097305098177 1.38664403839204525514392598679
3
-0.468787988438663412080404501151.36828202225060524339631698519 -0.468104608587966451955058675321 1.36309882555249380134910541460 -0.467909102935193813729827341301.36284170355894373658578526689 -0.467899374274474753499729975661.3627680850209174923696145195
4 5 6 7 8 9 10
-0.467896588260087871711112638661.3627644209583015340521492272 -0.467896449591436040517087255171.3627633716619036548612952215 -0.4678964098800759272523267310 1.3627633194349988293991589842 -0.4678964079035107365408096364 1.3627633044784694482399673681
pv
(7)
3.755833:31104103 26.24292723181837 PJ 76.7287159694094 92.9683130884962 20.42921161891182 3.542604537491395 71.40924770203088 71.41363478058405 3.495448039508563 20.15851163312566 70.2439801768494 70.40613461060159 20.09584447627069 3.492625640818941 70.17395663611611 70.17401734080785 3.491963560873392 20.0920523862253 70.15760657212744 70.15990731983832 20.09116227850327 3.491923372046278 70.15660956803084 70.15661043294261 3.491913936860466 20.09110823982103 70.15637656822064 70.15640936064984
(Note: Output has been edited to add symbols and to format.)
APPENDIX
D
(* Program to produce Figure 3 *) (* Figures 4 and 5 are obtained by changing the ranges of x and y *) (* The roots are at -1, 0 and gama=2 *) iter[x_,y_,lim_]:=Block[{c,z,ct}, g=l.-gama; gama=2; r[l]=-1.; r[2]=0; r[3]=gama; z=x +I y; Print [z]; ct=o ; it=O; While[(it0.2)&&(Abs[z]>0.2)&&(Abs[z-gama]>0.2), ++it; zz=zA2; z=zz(2z+g)/(3zz+2g*z-gama);
I;
Do[ If[ Abs[z-r[i]]<=.2,
ct=i*lO ,diverge=l],
(4311; Return[ct]
I; DensityPlot[ iter’[x,y,20] , {x,-0.9,1.7} , {y,-l.,l.}, PlotPoints->200, Mesh->False , PlotLabel-> FontForm[ “Fig.S:Basins in complex plane for roots at -1,0 and 2”) {“Helvetica-Bold” ,lO}], AxesLabel-> “x”, “y”
1
Pergamon
Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/95 $9.50 + 0.00
08957177(95)00073-9
Newton Raphson Method, Scaling at F’ractal Boundaries and MATHEMATICA A. K. BISOI” International Centre for Theoretical Physics 34100 Trieste, Italy (Received
September
1994; accepted
November
1994)
Abstract-The basins of convergence of cubic polynomials having real roots are studied using the Newton Haphson iterative method. The limiting value for the ratio of basin segments for equispaced roots is explained. An algorithm is presented for computing the basin boundaries on the real axis which obviates the necessity of taking recourse to extensive search. MATHEMATICA programs have been developed to help in the above study and also to depict the basins in the complex plane. Keywords-Fmctal basin boundaries, plex dynamics, Julia set.
Newton Taphson,
MATHEMATICA
programming,
Com-
1. INTRODUCTION Fractals are encountered in feedback processes, where iterations are involved, and are characterized by self-similarity. Fractal objects contain their replicas embedded in them, although reduced by a scale factor. If the scale factor is constant over replications, the fractal is strictly self-similar. However, if the scale factor is only asymptotically constant, we say the fractal is quasi-selfsimilar [I]. This scale factor is conceptually related to the fractal dimension. The fractal boundaries between the basins of convergence of roots of polynomials obtained by applying the Newton Raphson (NR) iterative method have been the subject of many studies [2-61. The use of computer graphics for depiction of these boundaries have produced many pictures of awe-inspiring complexity where the phenomenon of self-similarity is visually apparent in qualitative terms. Self-similarity involves repetition of a structure with a change of scale at each repetition. One wonders what is the scale factor involved in the NR method. What is the fractal dimension of the basin boundary? Must one always follow the time-consuming procedure of searching pixel by pixel to locate the basin boundaries? Guided by these questions, we study a family of cubic equations, having real roots and try to explain the scaling involved at the fractal boundaries. Such a study, involving computer experiments, has already been initiated by Becker and DSrfler [2]. We first summarize their result. Becker and Dorfler have studied the cubic polynomial f(X) = (Lr+ 1)X(2 - l), which has zeros at 2 = -1,
0 and 1.
(I)
The plot of the function is given in Figure 1.
The
*Permanent Address: Department of Computer Science and Application, Utkal University, Vani Vihar, Bhubaneswar 751004, India. The author would like to thank the International Centre for Theoretical Physics, Trieste for hospitality and A. Nobile and the members of Scientific Computing Section of ICTP for extending all possible help.
Typeset by 44-W 91 MCM 21/10-G
A. I<. BISOI
92 0.4
-0. I
-0.2 -0.3 -0.4
-0.5
-I
Figure 1. Plot of equation (1). number of alternating basins.
Intervals (-d,
-e)
and (d,e)
derivative f’(x) is zero at x = -d and x = d, where d = l/a. are presented with three broad regions of x; namely (i) x < -d
(ii) - d < x < d
I
0.5
0
and
contain
an infinite
Thus, in the first instance, we
(iii) x > d,
(2)
which are intuitively the most likely regions representing the basins of convergence of roots -1, 0 and 1, respectively. However, actual computation reveals that the basin of convergence of root at x = 0 is not the region from -d to d but a smaller region from -e to e. The point x = e is such that the Newton map
N(x) = returns this point after two iterations. This leads to
$$
Each iteration only changes the sign, i.e., N(e) = -e. e = 0.5%. ‘e
For the cubic under consideration, e = l/d. With 15 digit precision, Becker and Diirfler have identified 12 basin boundaries, gr to gr2, in the x region from -d to -e. These boundaries have been found through extensive search and computation in this region. We shall discuss in Section 2 how extensive search, which is very time-consuming, can be avoided. Between gr and grz, the basins of roots at x = -1 and x = 1 alternate with the width of the basins becoming narrower and narrower. The structure on the positive x axis is symmetric. Becker and Dorfler find that Sn - Sk-1 = 6.0. n+cCJQn+l - gn lim
(5)
They also remark that the above limit is always 6 for a cubic function having equispaced roots on the real axis. An explanation of these facts will be provided in Sections 2 and 3. Generalization to the case of non-equispaced roots is attempted in Section 3. Depiction of basin Boundaries in the complex plane has been dealt with in Section 4. Some concluding remarks are made in Section 5.
93
NewtonRaphsonMethod
0.4
0.2
0
-0.2
-0.4
-0.452
0.4
=x0
0.2
0
-0.2
-0.4
-I
-0.5
0
0.5
I
Figure2(b). The first two tangentsfor g2 < x0 < g3.
2. FINDING THE BASIN BOUNDARIES WITHOUT DOING EXTENSIVE SEARCH Let us consider the cubic equation with roots at -1,0 and 1 again. Figures 2a and 2b sketch the behaviour of the first two tangents obtained through NR iteration starting from x values which lie in the range gi < x0 < g2 and g2 < xc < gs, respectively. These diagrams have been produced through a program in MATHEMATICA [i’] given in Appendix A. Assigning a larger value to the variable iter in the program would result in the depiction of successive tangents converging to the root. This program has been inspired by a program [8,9] to depict the iterations of quadratic map. From -oo to gi = -d, it is the basin of attraction of the root at x = -1. Note that the tangent at f(gi) runs parallel to the x-axis and does not lead to any root. However, it being impossible to represent the irrational number fi in the computer, there is no danger of landing exactly on this starting value. As we move x0 to the right of gi, the point ti, where the tangent intercepts the x axis, moves more and more towards the left. We are in the basin
A. K. BISOI
94
of convergence of the root at x = 1 as long as the point ti lies to the right of -91 = d. As ti moves to the left of x = d, we enter the basin of convergence of the root at x = -1. Therefore, the point gs is such that it is mapped to x = l/t/5 m . a single NR iteration. Hence, the value of g2 is obtained by solving the cubic equation - 39; + 1 = 0.
2&g;
The value for the real root of equation (6) obtained using the MATHEMATICA
function NSolve
agrees with the g2 value of Becker and Dorfler. As we move x0 from g2 to gs, the point ti moves from -gi to -g2 on the positive z axis. This results in the point t2, the point of intersection of the tangent of the second iteration with the x axis, moving from -oo
to gi, (which is the
exclusive basin of the root at x = -1). From ga to g4, it is the basin of the root at x = 1 again. This structure of alternating basins continues an infinite number of times. The value of gs is obtained by solving the cubic equation 29; + 3g2g; - g2 = 0.
(7)
Thus, we have obtained the following recursive algorithm for finding the g’s. Starting with g1=
-5,
(8)
the successive g’s are obtained by finding the real root of the cubic 29; + 3g+1 g” - gi-1 = 0,
(9)
for i = 2,3,. . . , lc. Since the difference between the successive values of the g’s goes on decreasing, at some point i = Ic, it will no longer be possible to distinguish between g&i and gk. This leads to the limiting value
2g; + 39; -
0
gk =
(10)
or gk = 55
= +e.
We are interested in the limiting value of the ratio (gi -gi-i)/(gi+i -gi). of the differences can be replaced by the differential coefficient dgi_r/dgi.
s-1 Calculating dgi_i/dgi tion (II), we obtain
= -3gp
2gs _
1’
(11) For large i, this ratio From equation (9),
(12)
from equation (12) and substituting the limiting value of gi from equaQ-1 .limmdg,=
6
.
(13)
This explains the ratio obtained by Becker and Dijrfler. We can also understand why the basins appear as alternating segments. A MATHEMATICA program implementing the above algorithm for computing the g’s is given in Appendix B together with the output from this program. It has been chosen to calculate 30 number of g’s with a precision of 50 digits. It is to be noted that there exist an infinite number of g’s in the x interval from -d to -e and also in the interval from d to e. The higher the precision we choose, the more is the number of g’s we can identify. It is easy to perform calculations with arbitrary precision using MATHEMATICA. The ratios defined by Becker and DGrfler are also calculated with a precision of 25 digits.
95
Newton Flaphson Method
3. NONEQUISPACED
ROOTS
Consider a cubic equation with distinct real roots at CY,/3 and y such that (I: < p < y. The mapping x’ = (z - /3)/(p - a) relocates the roots at -l,O,l in the transformed variable z’, if the roots Q, p, y are equispaced, i.e., if (y - ,D) = (,0 - CX). Thus, we can understand why the limit given in equation (5) is the same for all cubic functions having equispaced roots on the real axis. The above transformation the mapped variable.
maps the nonequispaced roots to -1,0
and c = (y - @)I(/? - CX)in
Thus, it will suffice to study a cubic of the form f(z)
= (X + 1)x(x - c) to
generalize to the case of nonequispaced roots. Because of the asymmetrical
location of the roots,
the basin boundaries on the left and right hand side will no longer be symmetrically
situated.
Let us denote the basin boundaries on the left and the right hand sides by li and ri with i = 1,2,3,. . . . These points belong to the Julia set, the members of which do not converge to any root. The li and ri can be recursively generated starting from 11 and ~1 given by the solution of f’(z)
= 0; i.e.,
1
1-
=A-B 3
and
= -A+B 3
7.1
where A = c - 1 and B = dm.
’
The real root of the cubic equation N(l)(li) = Ti-1 determines li and similarly N(l)(ri) = li_l determines ri for i = 2,3, . . . , where iV(‘) represents a single iteration of the Newton map N(x) = Calculating
2x3 + (1 - c)s2 3x2 + 2(1-
(14)
c)x -c’
the succeeding l’s and the T’S involves finding the real root of the cubic equation 2x3 - (A + 3w)x2 + 2wAx + WC= 0,
(15)
which, for w = ri_1, returns x = li as the real root and for w = li_1, returns x = ri as the real root. A MATHEMATICA program implementing the above algorithm is given in Appendix C. The program computes the location of basin boundaries for various values of c. The limiting value of li would be such that lk = N(l)(?-k) = Nc2)(lk). This is an equation with a leading power of 6 in the unknown variable, and hence, an analytic solution is not possible in the general case. The ratios defined by (I) _ fj
-
4+1 - 4 &2 _ 1,,1;
p!” 3
=
rj+l rJf2
-
rj rj+l
(16)
are calculated for j = 1 to k - 3 and are given in Appendix C for various values of c. It is noticed that, for the case of nonequispaced roots, p’s with even and odd subscripts tend to separate limiting values. We conjecture that each of the basins of convergence would be scaling down by a constant factor. To verify, we define the ratios as JO _ 3
-1.
lj
.
3+3 - $+2
&+1
-
’
(17)
for j = 1 to k - 3 and observe from Appendix C that the values for sj are suggestive of an asymptotically constant value, s,(c) for each value of c. The variation of so0 with c can be fitted by a quadratic polynomial. It is to be noted that, although the initial values of p and s differ on the left and right hand sides, the asymptotic limits are the same. From the values of the ratio s, it is apparent that the basin segments of roots at -1 and 1 scale down by a factor of the same magnitude.
96
A. K. Brsor
4. DEPICTION
OF BASINS
IN THE COMPLEX
PLANE
To reproduce the depiction of basins in the complex plane using the Pascal programs of Becker and Diirfler, one also has to be aware of machine specific graphic primitives which must be interfaced to these programs. We find that MATHEMATICA
provides a short-cut and a ready
made answer to produce such plots through use of the DensityPlot DensityPlot
function.
A program to
the Mandelbrot set is given in the book by Crandall [lo]. We give in Appendix D
an adaptation of this program to plot the basins pertaining to the NR method in the complex plane.
0.5
0
-0.5
-1 -0.5
0
0.5
1
Figure 3. Basins in complex plane for roots at -1,0
1.5
and 2.
Figure 3 depicts the basins in the complex plane for the cubic with roots at -1, 0 and 2. The white region belongs to the basin of root at 2, the dark region to the root at -1 and the gray region to the root at 0. The picture is symmetric about the x axis. Therefore, in Figures 4 and 5, we plot the diagram for y > 0 only. Figures 4 and 5 depict enlargement of the left and right boundaries of Figure 3. For producing each of these figures, a grid of (200 x 200) points have been plotted. One can improve the resolution of the pictures by choosing a larger number of grid points. Consider the equation N(r)(z) x1 = -0.47664;
= ~1 again. The three roots of this cubic equation are
x2 = 1.39975 + 0.76830 i;
x3 = 1.39975 - 0.76830 i.
We have identified the real root x1 with 12, a point belonging to the Julia set on the negative x-axis. It turns out that 52 and x3 also belong to the Julia set. The location of the point 52 has been marked with a cross (+) in Figure 5 at the germination point of a bud on the right boundary.
Newton Raphson
n
I,,
-0.8
-0.9
Method
”
”
-0.6
-0.7
Figure 4. Upper half of left boundary
0.8
...
8,.
1
97
”
c -0.5
”
”
-0
L
of Figure 3.
.
*
1.2
1.4
Figure 5. Upper half of right boundary
”
of Figure 3.
*I
1.6
A. K. BISOI
98
5. CONCLUDING We have seen that it is not always necessary Analysis
of behaviour
of NR tangents
REMARKS
to make a step by step search for basin boundaries.
leads to the development
of an algorithm
for computation
of the basin boundaries on the real axis. This deterministic algorithm demystifies the occurrence of an infinite number of basin segments on the real axis. The behaviour of the points near the boundary may still appear to be chaotic. However, that is due solely to the finite precision of the computers and not due to any inherent chaos in the dynamical system itself. The execution of the algorithm In computing is felt that determination
requires
much less time compared
these boundaries,
exploring
the iteration
of boundary
points
to the time required
by a search method.
we have made use of only the real root of equation of the complex
roots of this cubic equation
in the complex
plane,
might
as we have demonstrated
(15).
It
lead to the
for one point
on Figure 5. The relation of s, to fractal dimension needs to be explored. It is hoped that this work will stimulate further application of MATHEMATICA to study the dynamics of NR method and search for deterministic algorithms to compute basin boundaries in the complex plane.
REFERENCES 1. S.D. Casey and N.F. Reingold, Self-similar fractal sets: Theory and procedure, IEEE Computer Graphics and Applications 14 (3), 73 (1994). 2. K.H. Becker and M. Dorfler, Dynamical Systems and Fractals, Cambridge University Press, (1989). Computers Math. 3. D.C. Arney and B.T. Robinson, Exhibiting chaos and fractals with a microcomputer, Applic. 19 (3), 1-12 (1990). 4. H.E. Benzinger, S.A. Burns and J.I. Palmore, Chaotic complex dynamics and Newton’s method, Physics Lett. A 119 (9), 441 (1987). 5. H.O. Peitgen, D. Saupe and F.v. Haeseler, Cayley’s problem and Julia sets, The Mathematical Intelligencer 6 (2), 11 (1984). 6. H.O. Peitgen, Editor, Newton’s method and dynamical systems, Acta Applicandae Muthematicae 13 (1988). 7. S. Wolfram, Mathematics, 2 nd edition, Addison-Wesley, (1991). 8. S. Wagon, Mathematics in Action, W.H. Freeman and Company, (1991). 9. A. Shiekh, A Mathematics essay on chaos, Lecture notes, College on Computational Physics, ICTP, Trieste, (1993). 10. R.E. Crandall, Mathematics
for the Sciences, p. 160, Addision-Wesley,
(1991).
Newton
Raphson
APPENDIX (*Program
to produce
Figure
2(a);
Figure
2(b)
A is obtained
tangents[f_,g_,x0_,displaylength_,xrange_:{-l.l,l.l}]:= Block[{z,p,plot, lplot, cplot,start, list, xmin,xmax}, {xmin, xmax} =xrange; gl=-0.57735; g2=-0.46560; eps=.Ol; z=N[xO]; start =z; Do[ { dnl={z, O),pb+ll={z,f[4), {n,l,displaylength,2} list =Table [ p[i], {i,displaylength}
zz[nl=z,z=g[zl1, 1;
lplot =ListPlot[list, PlotStyle -> Thickness[.OOOl], PlotJoined -> True, DisplayFunction -> Identity 1; cplot =Plot[ {y=f[x]}, {x,xmin,xmax}, PlotLabel->FontForm[“Fig. 2(a): The first two tangents for gl
Line[{{gl,O},{gl,f[gll}}l,
Line[({g2,O},(g2,f[g21))1, Line[{~-gl,O},~-gl,f[-g~l~~l, Li~[i{-g2.o},ip2,r[821}}1 ) 1 DeiaultFont->{ “Times-Italic”,lO}, DisplayFunction -> $DisplayFunction, PlotRange -> {xrange,{-.5,.5}}, Frame->True 1
1 x0=-0.49 ; iter=5; f[x_]=x ‘3-x (* Define cubic function *) g[x_]=2*xA3/(3*x”2 - 1) (* Define Newton map *) tangents[f,g,xO,iter,{-l.l,l.l}]
99
Method
by changing
x0 value *)
100
A. K. BISOI
APPENDIX
B
PROGRAM (* Calculate the basin boundaries without doing extensive search *) f[x_,w_]:=2*x”3+3*w*xA2-w; (*Define Newton map *) digit=50; k=30; sout=OpenWrite”math/cubic.out”] g[l]=N[-l/Sqrt[3],digit]; d=N[-l/Sqrt[5],digit]; d=“; i=l; space= “ ” ; textd” While[ i
, {j,k-3)
] ;
] ;
OUTPUT d=
-0.44721359549995793928183473374625524708812367192231
gj
pj
-0.57735026918962576450914878050195745564760175127013 7.256874166975170204865840 -0.46560062143367758369557882264206914551351750716552 6.179501149801548318846121 -0.45020147778247548915116189993227676825613611355286 6.029219709826801810191761 6.004851109839777676226105 -0.44770950581291017390380321914867329621880782129811 -0.44729618997943611050770656699237291959172375596111 6.000807997279660568129623 -0.44722735965776555512399767546695435157538351107654 6.000134651751163967293937 -0.44721588948213224738364426683460861182242551262533 6.000022441556852205031045 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
-0.44721397782909460094975366836322578275040527333734 6.000003740248317957732033 -0.44721365922144666848487585782979526457672255774767 6.000000623374409732649248 -0.44721360612020511504562984246847538285813648717433 6.000000103895726346363949 -0.44721359726999910897107968141990472300155464848525 6.000000017315954151919637 -0.44721359579496480016694813774716728620680520843668 6.000000002885992352010467 -0.44721359554912574940908250241415201561224632891307 6.000000000480998725150555 -0.44721359550815257430247960819118973631571735920319 6.000000000080166454186634 -0.44721359550132371178525990484819787143881297117573 6.000000000013361075697630 -0.44721359550018556803240516111586939946974676741518 6.000000000002226845949601 -0.44721359549999587740692979290561567200083082781657 6.000000000000371140991600 -0.44721359549996426230268390993756785678246343518122 6.000000000000061856831933 -0.44721359549995899311864292976882825256150182045676 6.000000000000010309471989 -0.44721359549995811492130276641642540273714007478203 6.000000000000001718245331 -0.44721359549995796855507940585794308751305566624677 6.000000000000000286374222 6.000000000000000047729037 -0.44721359549995794416070884576486968756126315043518 -0.44721359549995794009498041908269098162260018165942 6.000000000000000007954839 -0.44721359549995793941735901463566120268985919937807 6.000000000000000001325807 6.000000000000000000220968 -0.44721359549995793930442211389448957301746818878992 6.000000000000000000036828 -0.44721359549995793928559929710429430140956225609353 -0.44721359549995793928246216097259508947502680163677 6.000000000000000000006138 -0.44721359549995793928193930495064522081927410184747 [ Note:Output has been edited -0.44721359549995793928185216228032024270998207436278 to add symbolsand to format] -0.44721359549995793928183763850193274635843340559163
101
Newton Raphson Method
APPENDIX PROGRAM
(Generalization
C
of program of Appendix
B)
f[x_,w_]:=2xA3-(a+3w)xA2+2w*a*x+w*gama (*Define Newton map*) sout=OpenWrite[ “math/appcl4”]; k=lO; digit=30; space=” “; For[ gama=2,gama<4,gama++, Print[ “gama=” ,gama]; Write[sout,“gama=“,gama]; b=N[Sqrt[l+gama+gama”2],digit]; a=N[gama-l,digit]; r[l]=N[(a+b)/3,digit]; i=l; l[l]=N[(a-b)/S,digit]; While[ i
Write[sout,j,space,l[j],space,r[j]]];
Do[ (* calculate the ratios *) ul=N[l[m+l]-l[m],digit]; rhol[m]=N[ul/(l[m+2]-l[m+l])]; sl[m]=N[ul/(l[m+S]-l[m+2])];
ur=N[r[m+l]-r[m],digit]; rhor[m]=N[ur/(r[m+2]-r[m+l])]; sr[m]=N[ur/(r[m+3]-r[m+2])],
Do[ P rint[ rho l[nI,sPace,rhor[n],space,sl[n],space,sr[n]]; Write[sout,rhol[n],space,rhor[n],space,sl[n],space,sr[n]],
{m,k-3}];
{n,k-3}];
a=.; b=.; 1; (* end For *) Close[sout] OUTPUT gama=2
li
5 6 8 9 10
ri
-0.548583770354863530167205251213 1.215250437021530196833871917880 -0.476640869333724045703691984163 0.92050686281515245538451861863 -0.460538069165090573216401849777 0.90106639772494915960227270705 -0.459209517915383575472653379100 0.89629425599320413258572595171 -0.458877018093209225118648522210.89589288065135084229024118365 -0.458848935287665317313695075310.89579224077708478008583382472 -0.458841891016546856508041845800.89578373734252614997401883867 -0.458841295767567234288724171460.8957816042631496402834120321 -0.458841146448579826619728368140.8957814240139575015944589047 -0.458841133830822855086124200040.8957813787981735177199975455
&L’) ?
p3
s?) ?
SW 3
4.467726126371272 15.16134376614832 54.15139i13369052 61.76337392615517 12.12057131569808 4.073740090509531 48.42950009215275 48.43462730032162 3.995644993188464 11.88947410100077 47.30835199602091 47.41800172691879 11.83997879608157 3.988233737178288 47.20145158850765 47.20155592259313 3.986616226384538 11.8352030079332 47.17825062334561 47.18055754246566 11.83415908235831 3.986459506511234 47.17599041325278 47.17599262296456 3.986425236042335 11.83405790173216 47.17549885967337 47.17554775275868
A. K. BISOI
102
gama=3 i 1 2
li ri -0.535183758487996431039740422490 1.868517091821329764373073755823 -0.482748900024651977097305098177 1.38664403839204525514392598679
3
-0.468787988438663412080404501151.36828202225060524339631698519 -0.468104608587966451955058675321 1.36309882555249380134910541460 -0.467909102935193813729827341301.36284170355894373658578526689 -0.467899374274474753499729975661.3627680850209174923696145195
4 5 6 7 8 9 10
-0.467896588260087871711112638661.3627644209583015340521492272 -0.467896449591436040517087255171.3627633716619036548612952215 -0.4678964098800759272523267310 1.3627633194349988293991589842 -0.4678964079035107365408096364 1.3627633044784694482399673681
pv
(7)
3.755833:31104103 26.24292723181837 PJ 76.7287159694094 92.9683130884962 20.42921161891182 3.542604537491395 71.40924770203088 71.41363478058405 3.495448039508563 20.15851163312566 70.2439801768494 70.40613461060159 20.09584447627069 3.492625640818941 70.17395663611611 70.17401734080785 3.491963560873392 20.0920523862253 70.15760657212744 70.15990731983832 20.09116227850327 3.491923372046278 70.15660956803084 70.15661043294261 3.491913936860466 20.09110823982103 70.15637656822064 70.15640936064984
(Note: Output has been edited to add symbols and to format.)
APPENDIX
D
(* Program to produce Figure 3 *) (* Figures 4 and 5 are obtained by changing the ranges of x and y *) (* The roots are at -1, 0 and gama=2 *) iter[x_,y_,lim_]:=Block[{c,z,ct}, g=l.-gama; gama=2; r[l]=-1.; r[2]=0; r[3]=gama; z=x +I y; Print [z]; ct=o ; it=O; While[(it
I;
Do[ If[ Abs[z-r[i]]<=.2,
ct=i*lO ,diverge=l],
(4311; Return[ct]
I; DensityPlot[ iter’[x,y,20] , {x,-0.9,1.7} , {y,-l.,l.}, PlotPoints->200, Mesh->False , PlotLabel-> FontForm[ “Fig.S:Basins in complex plane for roots at -1,0 and 2”) {“Helvetica-Bold” ,lO}], AxesLabel-> “x”, “y”
1