For a complex polynomial $P(z)$, define $T$ by $T(z)=z-\frac{P(z)}{P^{\prime}(z)}$, where $P^{\prime}(z)$ is the first derivative of $P(z)$. We consider Newton method for finding the roots of polynomial $P(z)$, which is given by the recursive formula

z_{n+1}=T(z_{n}), \ n\in\mathbb{N},\ \ \ \ (16)

where $T$ is defined as above. Similarly, we set a new method (and called it S iteration) where $T$ is defined as above, by

z_{n+1}=(1-\alpha_{n})Tz_{n}+\alpha_{n}Ty_{n}, n\in\mathbb{N}, (17)

Now If the sequence $\{z_n\}_n\in\mathbb{N}$ (the orbit of the point $z_1$) converges to a root $z^{\ast}$ of the polynomial $P$, then we say that $z_1$ is attracted by $z^{\ast}$. The attraction basin of the root $z^{\ast}$ of the polynomial $P$ is the set of all starting points $z_1$ which are attracted by $z^{\ast}$. If the test polynomial is

P(z) = z^5 + z^4 + z^3 + z^2 + z,  (18)

For this test polynomials $P(z)$ above, we consider square domain is centered at the origin in the complex plane as follows


If $\alpha_n=\beta_{n}= 0.8$ for all $n \in\mathbb{N}$ in S iteration method (17), To generate the basins of attraction and to study the dynamics of methods (16) and (17), we divide the domains $D$ into 250 × 250 grids. If the sequence $\{z_n\}_n\in\mathbb{N}$ attempts a root of polynomial $P$ with accuracy of $10^−4$ in number of iterations $n \leq 13$, then the converging point $z_0$ is colored in a color assigned to $n$; otherwise, the point is colored in white. Figure 1a, b (see the image in attach or in the link for clear picture https://link.springer.com/article/10.1007/s00500-020-05038-9) represents the obtained basins of attraction by Newton method (16) and S-iteration method (17) for the polynomial $P$ defined in (18).

Now the attached pictures are generated by matlab and I want to generate the same pictures in mathematica.

enter image description here

  • 1
    $\begingroup$ People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find the meta Q&A, How to copy code from Mathematica so it looks good on this site, helpful $\endgroup$
    – Michael E2
    Commented Feb 7, 2022 at 0:12
  • 1
    $\begingroup$ JuliaSetPlot[z - #/D[#, z] &[z^5 + z^4 + z^3 + z^2 + z], z]? $\endgroup$
    – Michael E2
    Commented Feb 7, 2022 at 0:18

1 Answer 1


For Newton's method, this could look as follows:

p = z^5 + z^4 + z^3 + z^2 + z
znew = Together[z - p/D[p, z]]

cf = With[{
    p = p,
    num = N@HornerForm[Numerator[znew]],
    denom = N@HornerForm[Denominator[znew]]
   Compile[{{z0, _Complex}, {maxiter, _Integer}, {\[Epsilon], _Real}},
    Block[{iter, z},
     z = z0;
     iter = 0;
     While[Abs[z] > \[Epsilon] && iter < maxiter,
      z = num/denom
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True

n = 1000;
gridpts = Tuples[Subdivide[-10., 10., n - 1], 2] . {1., 1. I};
result = cf[gridpts, 13, 1. 10^-4];

ArrayPlot[Partition[N[result], n],
 ColorFunction -> "SunsetColors"

enter image description here

  • $\begingroup$ How to write the code for S iteration (17)? $\endgroup$
    – Junaid
    Commented Feb 7, 2022 at 10:23
  • $\begingroup$ Well, I did mean to get you started, not to solve the problem for you... $\endgroup$ Commented Feb 7, 2022 at 10:24

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.