On the robust Newton’s method with the Mann iteration and the artistic patterns from its dynamics

There are two main aims of this paper. The first one is to show some improvement of the robust Newton’s method (RNM) introduced recently by Kalantari. The RNM is a generalisation of the well-known Newton’s root finding method. Since the base method is undefined at critical points, the RNM allows working also at such points. In this paper, we improve the RNM method by applying the Mann iteration instead of the standard Picard iteration. This leads to an essential decrease in the number of root finding steps without visible destroying the sharp boundaries among the basins of attractions presented in polynomiographs. Furthermore, we investigate visually the dynamics of the RNM with the Mann iteration together with the basins of attraction for varying Mann’s iteration parameter with the help of polynomiographs for several polynomials. The second aim of this paper is to present the intriguing polynomiographs obtained from the dynamics of the RNM with the Mann iteration under various sequences used in this iteration. The obtained polynomiographs differ considerably from the ones obtained with the RNM and are interesting from the artistic perspective. Moreover, they can easily find applications in wallpaper or fabric design.


Introduction
The Newton's method is one of the most commonly known numerical algorithms-that at first was proposed around 1669 by Newton and then modified in 1690 by Raphson-for finding real roots of polynomials. Further, Simpson (1740) extended the method to solve nonlinear equations, and Cayley (1879) used the method for finding complex roots of polynomials. Cayley faced a strange and unpredictable chaotic behaviour-that he could not explain-of Newton's method in the case of a simple equation z 3 − 1 = 0 in the complex plane. The solution to that problem was then found in 1919 by Julia. Julia sets led next Mandelbrot in the 1970s to the discovery of the famous Mandelbrot set and fractals. Newton's method is also applied in optimisation for finding a local extreme of a differentiable function that is a root of its gradient [28]. The Newton's method, depending on a starting point, might be convergent or divergent, has a local quadratic convergence and is undefined for critical points (i.e. points for which the first derivative of a function equals to 0). This method is sensitive to a starting point and finds stable or unstable cyclic values. So then, Newton's method dynamic behaviour is very rich [16,34].
In the literature [15,23,24,30,35] there are known many modifications and improvements of Newton's method. Here, we would like to point out some of the most recent results in this area. In [10] various modifications of Newton's method with arithmetic, harmonic, and contraharmonic means are investigated. In [3,8,14] Newton's method was modified by using fractional derivatives instead of the classical ones. Wang and Tao [33] proposed a new way to construct a selfaccelerating parameter in Newton's method with memory.
The most striking effects that can be observed when finding the roots of the polynomial z 3 − 1 on the complex plane are the fractal boundaries among the basins of attractions caused by chaos and instability. In the case of complex polynomials, the Newton method, producing Newton fractals, is an important model of deterministic chaos [35] similarly as the logistic equation in the real case [22]. Basins of attraction form stunning images that are ones of the most beautiful fractals. Fractal boundaries in the form of braids can be decreased or even, to some extent, eliminated using the damped Newton method [11]. Recently, Kalantari [19] proposed the Robust Newton Method (RNM) that radically smooths the boundaries among basins of attraction that form sharp lines. That result is obtained by precise controlling of the step's length in the Newton method. Unfortunately, the RNM usually needs a very large number of iterations. Kalantari [19] showed that for a given accuracy ε > 0 the number of iterations needed by the RNM to converge to the root or critical point is O(1/ε 2 ). By using this bound, for example, for ε = 0.01, we obtain 1/ε 2 = 10,000. Thus, even for not so accurate computations for some points, we need to perform a large number of iterations, and this number increases with the increase in the accuracy. This can be treated as a drawback of the RNM.
In this paper, we propose the modification of the RNM, which relies on replacing the Picard iteration used in the RNM by the Mann iteration. This modification will decrease the average number of iterations needed to find the solution and, at the same time, the method will preserve its other features. We will analyse numerically different aspects of the proposed modification, e.g. basins of attraction, dynamics, and the average number of iterations. Moreover, by using various sequences of parameters in the Mann iteration, we will generate very complex and intriguing patterns from the modified method dynamics. These patterns might have artistic applications.
The paper is organised as follows. Section 2 presents the description of the RNM. In Sect. 3, modification of the RNM with the Mann iteration is given. Then, in Sect. 4, the plots of average number of iterations (ANI), convergence area index (CAI) and time generation characterising the RNM with Mann iteration are given. Also, polynomiographs presenting basins of attraction and dynamic properties of the considered root finding processes for a given set of polynomials are shown. Section 5 is devoted to artistic aspects of the RNM with Mann iteration related to the generation of very complex and beautiful polynomiographs from dynamics. The last section, Sect. 6, concludes the paper and points out future directions of research.

The robust Newton's method
The robust Newton method (RNM) was proposed to overcome the Newton's method main drawback-the lack of definition at critical points [19]. The RNM guarantees the reduction of the polynomial's modulus in successive iterations and has several differences compared to the classical Newton's method. Firstly, it converges globally, whereas the classical Newton's method converges locally. Secondly, it finds roots and critical points of polynomials, whereas classical Newton's method finds only the roots.
The RNM is based on the following facts: -Instead of finding zeros for polynomial p(z), the equivalent modulus minimisation of F(z) = | p(z)| 2 = p(z) p(z) is considered. Observe that minimisation of | p(z)| and F(z) are also equivalent. -In the classical Newton's method applied to F(z) at every iteration, a descent direction is chosen according to the geometric modulus principle (GMP) [18], and the step size is taken on the base of the modulus reduction theorem (MRT) [19]. The GMP gives a complete characterisation of all descent and ascent directions for | p(z 0 )| at an arbitrary point z 0 and additionally reveals a surprising geometric pattern showing partition of a unit disc, centred at z 0 , into sectors of ascent and descent directions. The GMP for polynomials [18] says that if a non-constant polynomial p(z) at z 0 equals to zero, then every direction at z 0 is an ascent direction for | p(z 0 )|. If p(z 0 ) = 0 then the cones of ascent and descent directions at z 0 divide the unit disc centred at z 0 into alternating ascent and descent sectors of equal angle π/k, where k ≥ 1 is the smallest index with p (k) (z 0 ) = 0. Such typical disc partitionings, except for possible rotations, for k = 0, 1, 2, 3 are presented in Fig. 1. Next, the MRT assures a reduction in the modulus of p(z) when moving from z 0 to z 1 by the amount proportional to | p(z 0 ) p (z 0 )| 2 . Summing up, the GMP enables to choose consequently a descent direction, whereas the MRT guarantees the decrease of F(z) along the chosen descent direction at every iteration step of the RNM. Here, it is worth to add that a random direction at z 0 has a fifty per cent chance to be a descent (or ascent) direction-what stresses the great importance of the GMP.
In the sequel, we recall the formulas defining the RNM, following [19].
Assume that p(z) = 0, z ∈ C and let us define [19]: Moreover, let us define: Additionally, let θ k be the angle given by the following formula: Now, the RNM for the starting point z 0 ∈ C is defined as: where and i denotes the imaginary unit, i.e. i = √ −1, and The term (u k /|u k |)e iθ k is called as the normalised robust Newton direction at z i [19] and C k (z i )/3 is the step size. Because the RNM finds roots and critical points, the stopping criterion has some other form than the one for the classical Newton's method, and is given by the following condition: where ε > 0 is the accuracy.
To summarise the RNM, in Algorithm 1, we present its pseudocode, and in Algorithm 2, the pseudocode for the N p (z) computation.

Algorithm 1: The Robust Newton Method.
Input: p -complex polynomial; z 0 ∈ C -the starting point; M -the maximum number of iterations; ε > 0 -accuracy. Output: A root or critical point of p. If the RNM does not find the solution in M iterations, the algorithm throws an exception.
Throw "No solution found."

The robust Newton's method with the Mann iteration
The problem of finding a fixed point of a given mapping is well-studied in the literature [1,6]. In these studies, one of the directions is using iterative approximation methods for finding the fixed points. The most famous method is the Picard iteration [27]: where T is a given mapping for which we search fixed points and z 0 ∈ C is a starting point. Another known iteration is the Mann iteration (defined in 1953 by William Robert Mann) [21]: (14) where α i ∈ (0, 1]. Depending on the type of operator T (contractive, expanding, etc.) and the space (Banach, hyperbolic, CAT(0), etc.) the sequence α i has to fulfil some additional conditions in order (14) to be convergent to a fixed point.
Let us notice that the Mann iteration for α i = 1 for all i reduces to the Picard iteration. In the literature, one can find many other iterations. A review of 18 different iterations and their dependencies can be found in [13].
At this point, we combine the RNM with the Mann iteration. We do so by using N p as T in (14) and obtaining: From now on, we will call the RNM with the Mann iteration, shortly, as M-RNM.
Now, let us write (15) in the following form: By looking at this formula, we can observe that at the (i +1)-th iteration the M-RNM moves z i into the direction given by vector N p (z i ) − z i scaled by α i . Let us note that N p (z i ) moves z i towards the solution assuring reduction in the modulus of p, so N p (z i ) − z i points in the same direction. Now, everything depends on α i . If α i < 0, then the point moves in the direction opposite to the modulus reduction. If α i > 0, the point moves towards the modulus reduction, but, depending on the value of α i , this point can move faster or slower, i.e. the modulus reduction is larger or smaller or if the value α i is too high one can obtain an increase in the modulus. After including (10) into (16), we obtain the following: From this form of the M-RNM, it is easily seen that we obtained the RNM with an additional damping factor α i . If we join α i with the term representing the step size, i.e. C k (z i )/3, we see that α i changes the step size in the original RNM. To see how α i affects the root finding process of the RNM let us consider the following example. We search for the roots of p(z) = z 3 − 1 and we use α i = const = α. We presented the orbits for the starting point z 0 = 0.9 + i in Fig. 2 for various values of α. We fixed the maximum number of iterations to 500. According to our previous observation, the orbit for α = 1 presents the orbit for the RNM introduced in [19]. When we look at the various orbits, we can see that for α = 1, the orbit forms a smooth curve. If we increase the value of α, we see that the curve becomes a polyline. In the plot presented in Fig. 2a, we can see that for α > 1.0, we obtain a greater modulus reduction and the length of the corresponding orbit is shorter (see Table 1) than in the case of α = 1. However, if the value of α is too high, then the orbit's points jump from the neighbourhood of one root to the neighbourhood of the other root and, as a consequence, this method does not find the solution (see Fig. 2b). Therefore, from this example, we see that by a proper selection of α we can make (15) faster than the original method proposed in [19].
The selection of α is not an easy task, because the best value of α is different for different polynomials and there is no obvious dependency between α, polynomial p and the orbit's length. In each step, we could try to find the best value of α by minimising the modulus, i.e.
This is a similar equation as in the gradient descent method with the exact line search [31]. By using (18), we can always find α that assures modulus reduction. This is since for α = 1, we get the RNM which reduces the modulus [19]. Therefore, if, in the current step, we do not find α = 1 that reduces the modulus better than the RNM, then we use RNM (α = 1). This method may give the best value of α, but it is computationally too expensive in practice. Instead, we can use a method similar to the backtracking line search known in optimisation [31]. The proposed method is presented in Algorithm 3.

Algorithm 3:
Finding of the best value of α in the Mann iteration for (i + 1)-th step.
Input: p -complex polynomial; z i -root approximation from the i-th step; α -the starting value; s ∈ (0, 1) -the scaling factor. Output: The value of α for (i + 1)-th step.
Using the same example as in the case of the constant α, i.e. the contour plot of p(z) = z 3 − 1, let us see how the orbits look like when we use Algorithm 3. In this example, we take 19.65 as the starting value of α, and we generate orbits for various values of s. The obtained orbits are presented in Fig. 3. The starting value of α is the same as for the constant α used in Fig. 2b, where the orbit has jumped between neighbourhoods of the roots. After using Algorithm 3, we see that the orbits do not jump. Instead, they tend towards the root. For various values of s, we can see that the first step is the same, but the orbits differ from the second step. The lengths of the orbits are gathered in Table 2. From the data, we can see that the orbits for s equal to 0.1, 0.3, 0.5 and 0.7 are much shorter than in the case of the RNM (orbit's length equal to 80). Only for s = 0.9 we obtain longer orbit, but the method has converged in contrast to the constant α value equal to 19.65.
By introducing the Mann iteration in the RNM, we can ask whether the new method is convergent or not. The answer is affirmative when we assume that α i is selected so that in each iteration we get modulus reduction. We can make such an assumption, because if in Table 2 Lengths of the orbits presented in Fig. 3 s Orbit's length 0. the current iteration we do not find α = 1 that reduces the modulus, then we can use α = 1, i.e. the standard RNM for which it has been proven that it reduces the modulus [19]. In this situation, by using the M-RNM for any starting point, we obtain a decreasing sequence bounded below, so it is convergent.

Numerical results
In this section, we present the numerical results obtained with the help of the method proposed in Sect. 3 for α i = const = α and for α i computed by using Algorithm 3. The results will be divided into two categories. The first category consists of images showing basins of attraction and dynamics, the so-called polynomiography [17]. In the second category, we will present various numerical measures, such as the average number of iterations [4], convergence area index [5] and polynomiograph's generation time [12]. Nowadays, polynomiography is an essential part of modern analysis of the quality of the root-finding methods [26]. Depending on the selected colouring method, we can visualise various aspects of the considered root finding method with the help of polynomiography. There are two standard methods of colouring polynomiographs.
In the first method, called basins of attraction, each root has its distinct colour, but other than some fixed colour, usually black, that represents the case in which the method has not converged to any root. Then, for the last point z i+1 , we search for the closest root (by using the modulus metric) and colour the starting point with the colour of the root if the number of performed iterations was less than M, and with the fixed colour (black) otherwise. For the RNM and M-RNM, we add one more colour to colour the critical points. The algorithm for generating basins of attraction is presented in Algorithm 4.

Algorithm 4: Polynomiograph generationbasins of attraction.
Input: p -complex polynomial; {r 1 , r 2 , . . . , r n } -the roots of p; A ⊂ C -area; M -the maximum number of iterations; ε > 0 -accuracy; α i ∈ R + -the parameters sequence in the Mann iteration; colours [1..n] -colours for the roots; c cr -the colour for critical points; c nc -the colour for non-convergent points. Output: Polynomiograph for the area A.
is Root = f alse In the second method, the so-called iteration colouring, which depicts the speed of convergence and the dynamics of the considered method, the starting point is coloured by linearly mapping the number of performed iterations on the colour in the given colour map. The algorithm for generating this type of polynomiographs is presented in Algorithm 5.
As we mentioned at the beginning of this section, we used three numerical measures in our comparisons. The first measure is the average number of iterations (ANI) [4]. To compute ANI, we need a polynomiograph generated with the iteration colouring. In this polynomio- colour z 0 with colours [ j] graph, each point corresponds to the number of iterations needed to find a root. As the name suggests, ANI is computed as the average number of iterations in the given polynomiograph. The second measure we use is the convergence area index (CAI) [5]. This measure can be computed from basins of attraction or a polynomiograph generated with the iteration colouring. CAI is the ratio of the number of convergent points N c to the number of all points N in the considered polynomiograph, i.e.
From the definition of CAI we can see that it is between 0 and 1, and it gives us information about the percentage of polynomiograph's area that has converged to roots. The last considered measure is the time of generation of the polynomiograph [12]. It gives us information about the real time of computations.
In the numerical examples presented in this section, we used the following three polynomials, which are the most commonly used in the literature: To generate the polynomiographs, we used the following parameters: A = [−3, 3] 2 , M = 250, ε = 0.001, and image resolution of 800 × 800 pixels. The colour map used for the iteration colouring is presented  Fig. 4. In the basins of attraction, the lack of convergence is depicted by the black colour, whereas the critical points by yellow colour and the roots by any other ones. In the experiments with the constant α parameter used in the Mann iteration, the α was taken from (0, 50] with the step equal to 0.01. In case of the variable α parameter, i.e. computed with Algorithm 3, several values for the starting α were investigated and the s parameter was taken from (0, 1) with the step equal to 0.01. All experiments were performed on the computer with the following specifications: Intel i5-9600K (@ 3.70 GHz), 32 GB DDR4 RAM, NVIDIA GeForce GTX 1660 Ti graphics card with 6 GB GDDR6 SDRAM and Windows 10 (64 bit). The software for polynomiograph's generation was implemented in the C++ programming language with OpenGL (Open Graphics Library) and GLSL (OpenGL Shading Language). The computations of polynomiographs were performed on the graphics card with the use of shaders written in GLSL. Fig. 5, the basins of attraction for p(z) = z 3 − 1 for different values of α are presented. Since this polynomial has three roots and one critical point, we can see four corresponding basins. Let us note that for α = 1.0, we cope with the Picard's iteration. Thus, we obtained the polynomiograph for the RNM. The boundaries of the basins are sharp and not fractal-like in the case of classical Newton's method. Moreover, in the centre of this polynomiograph (point 0 + 0i), there is a yellow dot. This is the critical point of p. Furthermore, we can see that the decrease in α causes that fewer points converge to the roots. Indeed, when α tends to zero, the black colour floods from the proximity of the critical point until the total vanishing of the red, green, and blue areas. The yellow dot in the centre stays forever. On the other hand, the increase in α causes something different and the dynamics of changes is also different. We can observe that starting from about α = 10.4, the basins' regularity is broken at the borders and tends to be larger. From α = 17.78 the number of convergent points dramatically changes and, finally, these points disappear. Instead, some points converge to the critical point (yellow areas). From about α = 34.0 the points that converge to the critical point start to appear randomly.
In Fig. 6, the dynamics for the same polynomial and the same values of α are presented. In such a way, we can observe how fast the proposed method converges for the plotted starting points according to the colour map from Fig. 4. However, because we have chosen the colour map that is not monotonic, and we did it to emphasize the dynamics and clearly see the differences, it is a little bit difficult to see the increase in the speed of convergence for appropriate points without looking at the colour map. However, yet one can observe that, as in the case of basins of attraction, for the same values of α the same points are non-convergent (black areas). Additionally, for α = 1.0 one, can see that the closer the point is to one of the roots, the faster the convergence is (what is in accordance with intuition). When we tend with α to zero, we see that the colours tend to move to the right on the colour map (see Fig. 4). It means that the overall convergence becomes slower. On the other hand, when we increase α we move to the left on the colour map, so the overall convergence becomes faster. However, from about α = 10.4, the convergence becomes slower again. This is the visual confirmation of the facts discussed in Sect. 3 that depending on the value of α the speed of convergence changes, and we see here that this is done globally for all points.
To make the considerations complete, in Fig. 7, the plots of the commonly used measures: ANI, CAI and generation time (in seconds) are presented. From Fig. 7a, one can see that the average number of iterations is the smallest for α = 9.32 and equals to 5.63, which is considerably smaller than in the case of the RNM for which the value of ANI is equal to 84 the polynomiograph much faster (using less number of iterations). As one can expect, all these plots confirm the visual results presented in Figs. 5 and 6.

Variable α sequence
In this subsection, we analyse the M-RMN algorithm for the case of non-constant α sequence. We take the starting values of α as 20 and 50, and a discrete scale values s from the interval (0, 1). As a result of the experiments the images of basins of attractions, presented in Fig. 8, the images of dynamics, shown in Fig. 9, for z 3 − 1 and numerical characteristics of ANI, CAI and time generation, given in Fig. 10, were obtained. From the sets of the obtained images, the most representable ones were chosen. In the images in Fig. 8a, e dominates the colour black, what means non-convergence and denotes very bad behaviour of the M-RNM algorithm. The remaining images, from Fig. 8b-d, show that the M-RMN method is convergent. In these images, mostly three colours are clearly seen: red, green, and blue, all denoting basins of attractions, and small yellow areas related to the critical point.
The obtained basins look different compared to those generated using the RNM or by the classical Newton's method. They have complicated shapes with rotational symmetry, split into several bigger and smaller areas.
Their structure is not fractal with typical braids as in the classical Newton's method and different to the three triangles creating basins as in the case of the RNM algorithm. Thus then, the M-RMN in the considered cases of variable α sequences is less robust in comparison with the RMN or M-RNM with constant α sequences.
For the same choice of α and s parameters as for the basins of attraction, in Fig. 9, dynamics of the M-RNM root finding process with variable α sequence is presented. In Fig. 9a, e dominates the colour black what means the lack of convergence, whereas in Fig. 9b-d bright bluish colour is mostly seen what means, according to the used colour map (Fig. 4), a very low number of iterations, i.e. very high speed of convergence.
In these experiments, we obtained the following best results: ANI = 5.54, CAI = 1.

Constant α sequence
In Fig. 11, the polynomiographs' sequence showing the evolution of the five basins (four for the roots and one for the critical point 0 + 0i) of z 4 − 1 depending on α is given. The basins of the roots are coloured by four colours: red, green, blue, and magenta. Black colour denotes the lack of convergence, whereas yellowconvergence to the critical point. In Fig. 11c- Further information concerning the speed of convergence of the method can be extracted from the poly-nomiographs presented in Fig. 12 that show dynamics of the considered root-finding process for the same sequence of α values as in Fig. 11. Figure 12  on the diagonals. If α is decreasing from 1.0 more and more nonconvergency appears (Fig. 12a, b). When α is increasing from 1.0 up to about 32.0 the colours in Fig. 12e-h are moving to the left part of the colour

Variable α sequence
In this subsection, we show some examples of the M-RMN for variable α sequence. As in Sect. 4.1.2, we fixed the starting values of α as 20 and 50 since they are the most representative. In Fig. 14, the examples of basins of attraction for the considered polynomial are presented for these starting values of α and different values of parameter s. As in the previous polynomial case, the values of s were chosen to give the most diversity of this presentation. It is worth mentioning that in the case of starting α = 20 for all values of s the polynomiograph looks the same, whereas for starting α = 50, it changes from the one with much non-converging points (black areas, Fig. 14b) via the lack of them (Fig. 14c) and back again (Fig. 14d).
Similarly as in the previous subsection, in Fig. 15, the dynamics polynomiographs are presented for different starting values of α and s (the same as in Fig 14). For starting α = 20, the dynamics remains the same for all values of s. On the other hand, for starting α = 50 we can see the diversity of convergence for different values of s.
Finally, in Fig. 16, the plots of ANI, CAI and generation time are presented for the three different starting values of α. As one can see for starting α equal to 20 and 30 the number of iterations is constant and so the convergence area index and time of generation. In the case of starting α = 50, we can observe some peaks in the plots. These peaks represent the polynomiographs with many non-convergent points. The higher the peak, the more the non-convergent points. Moreover, we can observe that the best choice of s is about 0.4 since it assures the shortest generation time and the best convergence for all tested starting values of α. In more detail, for the variable α sequence the shortest generation times were obtained as follows: for α = 20, the time was oscillating around 0.033 s, for α = 30, the time was oscillating around 0.067 s, and for α = 50, the shortest time was 0.040 s.

p(z) = z 5 + z 4.3.1 Constant α sequence
In the last example, we will consider the polynomial p(z) = z 5 + z. It has five distinct roots and four critical points, i.e. −0.4728 + 0.4728i, −0.4728 − 0.4728i, 0.4728 − 0.4728i, 0.4728 + 0.4728i. The basins of attraction obtained for this polynomial are presented in Fig. 17. The roots are coloured by red, green, blue, cyan, and magenta colours. Firstly, let us look at the basins of attraction for the RNM with the Picard iteration (Fig. 17c). We can easily observe that many of the points have not converged to any root or critical point. Most of the converging points are converging to one of the roots, namely to 0 + 0i (cyan colour). For the Picard iteration, we can increase the number of converging points by simply increasing the value of M obtaining convergence of all points in the considered area for M = 1255, which is a high value. If we decrease the value of α in the Mann iteration below 1.0, then we can observe that the M-RNM converges to only one of the roots and that the basins for this root shrunk. Thus, the method behaves worse than the RNM. If we increase the value of α above 1.0, then we see, the opposite behaviour, i.e. the non-convergent points for α = 1.0 start to converge to the four roots other than 0 + 0i. For α = 3.0, we see almost full coverage of the area with the convergent points (see Fig. 17e), but for α = 10.0, we obtain the situation in which all points converge (see Fig. 17f). We obtain this by using M = 250 in the Mann iteration, which is considerably smaller than M = 1255 for the Picard iteration. By increasing α further, we can observe that some points lose convergence to the root represented by the cyan colour. For α = 22.0, we obtain the situation in which none of the points converges to the 0 + 0i root. By increasing α further, we see that some chaotic behaviour occurs. We still see some visible basins of attraction for the roots, but they are much smaller than for the lower values of α. Most of the points converge to a random root.
When we look at the dynamics presented in Fig. 18, we can observe a very interesting behaviour of the M-RNM. Similarly to the basins of attraction, we start with the polynomiograph for the RNM with the Picard iteration (Fig. 18c). From this polynomiograph, we can observe that very fast convergence is obtained near the 0 + 0i root. Moreover, the farther away from the root, the more iterations the method needs to perform to reach this root. When we decrease the value of α in the Mann iteration below 1.0, the speed of convergence is also decreasing, i.e. we need to perform more iterations to find the root. The number of iterationsdepicted by colour-changes radially from the 0 + 0i root. If we take values of α greater than 1.0, then we can observe a considerable reduction in the number of performed iterations. For instance, in the central part, we see a more dark pink colour, which means that to reach the root, we need fewer iterations than in the case of α = 1.0. In the same time, we can observe that the dynamics increases. From Fig. 18f, we see that the M-RNM method needs a very small number of iterations to find the roots (we see many points with dark pink and light violet colours). However, when we increase the value of α further, we see that in the cross-like area in the centre of the polynomiograph, the number of iterations increases and the method loses its convergence (black colour represents the maximal number of iterations), see Fig. 18h. After a certain threshold, the M-RNM starts to behave increasingly chaotic. In Fig. 18i-o, we see this behaviour. Moreover, we can observe that for the points that converge a similar number of iterations is needed to find the roots, i.e. we see light violet and pink colour, which means that we need less than half of the maximal number of iterations (M = 250). The plots of ANI, CAI and polynomiograph generation time (in seconds) for p(z) = z 5 +z are presented in Fig. 19. The values of these measures for the RNM with the Picard iteration (α = 1.0 in the Mann iteration) are equal to 229.77, 0.20 and 1.02 s, respectively. When we look at the ANI plot (Fig. 19a), we see that for values of α less than 1.0 the value of ANI increases its value obtaining the maximal value of 250 for α < 0.07. For values of α greater than 1.0, we see that ANI decreases with the increase in α. The local minimal value of 25.23 is attained at α = 15.59. Then, ANI increases obtaining the local maximal value of 197.11 for α = 39.46. For α > 39.46, we can observe another decrease in ANI. The minimal value of 23.92 (which is the global minimum) is attained at α = 48.71. By comparing the results for α > 1.0 with the result obtained for the RNM (229.77) we see that for all the values of α we get a lower value of ANI. Now, let us look at the CAI plot presented in Fig. 19b. From the plot, we see that for α < 1.0 we obtain results worse than for the RNM.
The lowest value of CAI equal to 0.0 is obtained for α < 0.07. If we change α in the opposite direction, i.e. we increase its value above 1.0, we see that the value increases violently obtaining the maximal value of 1.0 for α = 4.16. The value 1.0 of CAI does not change until α = 17.79, when it starts to decrease. The local minimum of 0.31 is attained at α = 39.46. Then, again the value of CAI starts to increase obtaining the value near 1.0 (≈ 0.99) for α > 40.7. Additionally, when we compare the results for α > 1.0, then we see that the obtained results for the M-RNM are better than in the case of the RNM (0.20). In the last plot, presented in Fig. 19c, we see the polynomiographs' generation times. The overall tendency observed in this plot is very similar to the one obtained for ANI (Fig. 19a). One interesting thing that we can observe by comparing the plots for ANI (Fig. 19a) and times (Fig. 19c) is that the value of α for which we get the lowest value best value of α that assures good convergence and performance, we need to consider all three measures (ANI, CAI, time) simultaneously.

Variable α sequence
In this subsection, we present some examples obtained for the M-RMN with a variable α sequence. The  Fig. 20. For the starting α = 20 (see Fig. 20a, b), the basins of attraction for the majority of the s values look like the basins presented in Fig. 20a. The method converges to all roots, and there are no visible non-converging points. The loss of convergence of the method appears for the s values near 1, where the method loses its convergence in the cross-like region, see Fig. 20b. For the other starting α value (see Fig. 20d, e), the basins of attraction look different, i.e. we do not see the cross-like region. However, also in this case, for the majority of the s values, the basins look like the ones presented in Fig. 20c. For values of s near 1 and for some single values of s in the rest of the (0, 1) interval, e.g. for s = 0.36, we loose convergence of the method, see Fig. 20d, e. The polynomiographs showing the dynamics and the speed of convergence for the same values of the parameters (starting α and s)-like in the case of basins of attraction-are presented in Fig. 21. When we look at starting α = 20, we see that the dynamics and the speed of convergence change only in the cross-like region. In the other regions, the polynomiographs do not change. The speed of convergence, compared to the RNM (Fig. 18c), is faster for the majority of the s values. Only for values of s near 1 the speed in the cross-like region is worse. For starting α = 50, we see that for the polynomiographs with good convergence (Fig. 21c) the speed of convergence is very fast. For the points that converge in the other cases (Fig. 21d,  e), the speed of convergence is faster or comparable to the RNM. Moreover, when we compare the polynomiographs for starting α = 50, and the one for a constant α = 50 (Fig. 18o), we see that in the variable case there is no chaotic behaviour that is visible for the constant case.
The plots of ANI, CAI and the generation time for three starting values of α, namely 20, 30 and 50, are presented in Fig. 22. From the ANI plot (Fig. 22a), we can see that in all three cases the method obtains larger values of ANI for s near 0 and 1. Moreover, for starting α values equal to 30 and 50 we observe some peaks, which correspond to the exemplary basins of attractions presented in Fig. 20d. When we compare the obtained results with the result for the RNM, i.e. 229.77, we see that for the variable α sequence, we obtain much lower values of ANI. Only in the peaks visible for starting α = 50, we can see the values close to 229.77.

(c) (d) (e)
The minimal values of ANI-for the three considered starting α values-equal to 16.298, 9.789 and 3.671 were attained for s equal to 0.45, 0.3 and 0.18, respectively. By looking at the plots of CAI (Fig. 22b), we can observe that for most of the s values the value of CAI is equal to 1.0. The values less than 1.0 are only obtained for s close to 1 and in the visible peaks that appear in the same points as the peaks in the ANI plot. By comparing these results with the value of CAI obtained for RNM, i.e. 0.20, we can see that the convergence of M-RNM with the variable α sequence is much better. From the last plot (Fig. 22c), i.e. the generation time plot, we can observe a very similar tendency as in the case of the ANI plot. The only difference is that for the values of s near 1 we can observe a very large time difference compared to the other peaks. The best times in all three considered cases are much shorter than in the case of the RNM (1.

Comments
We only consider positive values of parameter α in the M-RNM algorithm because for negative values of α convergence cannot be obtained. For α < 1, when α is getting smaller and smaller, the M-RNM behaves increasingly badly-more and more points of nonconvergence appear in polynomiographs. The M-RNM behaves better than the RNM for α > 1, i.e. more points converge, and the number of iterations is smaller. The optimal value of α is such that the ANI measure is as low as possible, and the CAI index simultaneously is as close as possible to 1. Such cases occur for z 3 − 1 and z 4 − 1 when α attain 9.32 and 20.47 for which ANI measure attain global minimum 5.63 and 5.41, respectively. For other polynomials, the optimal value of α can be taken as the one for which ANI measure attains local minimum and CAI index is close to 1. Such a situation occurs for polynomial z 5 +z. That example shows that because of mutual nonlinear dependency between parameters ANI, CAI and time versus α it is not possible to formulate general recommendations on how to choose optimal α to achieve the best improvement of the M-RNM algorithm. However, the best values of α for a specific polynomial could be determined experimentally via analysing of ANI, CAI and time plots.
From the examples for the variable α given by Algorithm 3, we can observe that also, in this case, we cannot give general recommendations on the choice of the starting α and s values. The relationship between these parameters is also nonlinear, but one should avoid the values of s near 0 and 1 because for these values the method obtains worse results than for the other values in (0, 1).
Similarly, like in the case of the RNM, the boundaries between basins of attraction for M-RNM are sharp in a wide range of α values for all considered polynomials, and this method is stable. However, after some threshold value (different for each of the polynomials), the basins of attraction look chaotic and the method losses its stability.

Artistic patterns from dynamics of the M-RNM
Dynamics of discrete dynamical systems, besides its essential meaning in the analysis of these systems, can also have artistic applications. Graphical presentation of the phase portraits of many discrete dynamical systems can reveal artistic patterns of aesthetic value. We can find many examples of such patterns in the literature, e.g. Ouyang et al. [25] presented spiral patterns that are coloured by the dynamics of the dynamical system that is compatible with some symmetry group, Lu et al. [20] and Chung et al. [9] proposed methods for wallpaper patterns creation from dynamics, and patterns from the dynamics of combined root-finding methods were presented by Gdawiec [12].
When we look at the dynamics of polynomiographs generated by the RNM, e.g. Figs. 6c, 12c or 18c, we can see rather boring patterns-from the artistic point of view. However, when we look at the patterns generated by the Mann iteration in conjunction with the RNM method (see, for example, Fig. 12), we can see that very interesting patterns can emerge. The images presented in Sects. 4.1.1, 4.2.1 and 4.3.1 were obtained for constant sequences of α i . In the Mann iteration, we are not limited to only constant sequences, especially if we are interested in artistic applications. Therefore, in the rest of this section, we will introduce some nonconstant sequences that will generate very interesting patterns from the dynamics of the M-RNM.
Before presenting examples of artistic patterns, we need to consider two technical subjects that greatly influence the appeal of the patterns presented in the polynomiographs. Let us look at Fig. 23a which was obtained for p(z) = z 5 + z, α = 17.7, M = 200, A = [−3, 3] 2 and ε = 0.001. When we look closely at the cross pattern, we see many unpleasant artefacts that look like jagged shape, noise or chaos. This unpleasant effect is caused by the fact that the image has a discrete nature. We sample the signal with a high frequency, i.e. we discretize the space that contains many tiny details that cannot be captured by the discretization process. This is a well-known problem in computer graphics, and it is called aliasing effect [2]. To reduce this effect, one can use one of the many methods of anti-aliasing. The simplest anti-aliasing method is supersampling [2]. In this method, we render the image at a higher resolution than the image to be displayed, usually 2, 4 or 8 times higher. And then, we average the colours in non-overlapping blocks. The block size depends on the resolution's magnification factor, i.e. if f is the factor, then the block size is equal to f × f . The result of using anti-aliasing method to the pattern from Fig. 23a is presented in Fig. 23b, c. We see that the appeal of the pattern in the areas with high frequencies is much better. Now, the pattern does not look noisy or chaotic.
In all examples presented in this section, we will use this anti-aliasing method.
The other technical subject that we need to consider when we want to use polynomiographs generated by the proposed M-RNM is choosing a proper colour map. In Sect. 4 in the polynomiographs with iteration colouring, we used only one colour map. However, when we think about artistic applications, we can use various colour maps-even for the same polynomiograph-to obtain very interesting results. In Fig. 24, we see an example of the same polynomiograph ( p(z) = z 3 − 1, α = 17, M = 200, A = [−3, 3] 2 and ε = 0.001), but coloured with various colour maps. From the images, we can see that by a proper choice of the colour map we can, for example, emphasize different parts of the same pattern (e.g. see Fig. 24b, e), show banding effect, or hide it (e.g. see Fig. 24a, f). Moreover, the proper colours can create the depth, movement, mood and harmony of an image [29,32]. Therefore, in the all presented examples in this section, we use various colour maps that greatly increase the polynomiographs' artistic values.
The first sequence that we will consider is very simple, and it is based on the switching of several constant values. Let m ∈ N. Then, we define α i in the following way: where q 0 , q 1 , . . . , q m−1 ∈ R.
As an example of the use of (20) let us take m = 3 and three constant values In Fig. 26, we see the obtained polynomiographs. From these polynomiographs, we can observe that although we used parameters which led to non-attractive patterns, their combinations in (20) allow us to obtain a variety of interesting patterns. Each of the obtained patterns has three visible braids, but they significantly differ in details.
Another way of defining α i for artistic applications is the use of linear interpolation between two given values, i.e.
where q 0 , q 1 ∈ R and M is the maximum number of iterations. Let us consider two values −17.0 and 30.0 and the following parameters for the polynomiographs: p(z) = z 4 + z 2 − 1, A = [−3, 3] 2 , M = 50 and ε = 0.1. Similarly like in the previous example, let us start with the polynomiographs generated for the two considered values, see Fig. 27. These two polynomiographs present two completely different patterns. In Fig. 27a, we see a very simple and small pattern in the centre of the polynomiograph, whereas in Fig. 27b, the pattern occupies the whole area and it is rich in details. When we modify α i by using the linear combination of these two values: (a) q 0 = −17.0, q 1 = 30.0, (b) q 0 = 30.0, q 1 = −17.0, then we obtain the polynomiographs presented in Fig. 28. We see that both polynomiographs differ from the ones presented in Fig. 27. However, we can notice that some characteristic structures were inherited from the original patterns, i.e. in the pattern from Fig. 28a, we see the small pattern from Fig. 27a, and in the pattern from Fig. 28b, we see that the overall shape is similar to the pattern from Fig. 27b.
The linear interpolation applied in (21) can be used in many different ways to define α i . For instance, let us notice that when we generate a polynomiograph the When we look at the templates in Figs. 31e and 12k, we see that they have a fourfold symmetry and four axes of symmetry (two diagonal, one horizontal and one vertical). Thus, using this kind of templates, we can create a wallpaper pattern of p4m symmetry (Fig. 32a, b) by repeating the template on a square lattice. If we place template from Fig. 24c on a square lattice, then we see that it has no translational symmetry nor glide symmetry. However, when between two adjacent cells in the square lattice we add a mirror reflection in the horizontal direction, then we obtain a pmm symmetry pattern (Fig. 32c). The template in Fig. 33a has a twofold symmetry and two axes of symmetry (one horizontal and one vertical). These types of symmetry of the template will not give the p4m symmetry if we repeat it on a square lattice. Instead, we obtain a pmm symmetry pattern (Fig. 32d). A twofold symmetry is visible in the templates from Fig. 33b, c. These templates lack the horizontal and vertical mirror symmetries present in the template from Fig. 33a, so we cannot repeat them on a square lattice to obtain a pmm symmetry. However, we can use the mirror symmetries in horizontal and vertical directions between two adjacent cells in the square lattice to obtain a cmm symmetry patterns (Fig. 32e, f).  Fig. 33d), (c) p(z) = z 4 + z 2 − 1, A = [−3, 3] 2 , M = 50, ε = 0.1 and α i is given by (22) with T = 0.789, (see Fig. 29a).
The template from Fig. 26a cannot be directly used to obtain a frieze symmetry pattern, because it lacks translational or glide symmetry after repeating it in a horizontal direction. Nevertheless, if we rotate it by 90 degrees and repeat this template in a horizontal direction, then we obtain a pm1 symmetry pattern (Fig. 34a). When we look at the templates in Figs. 29a and 33d, we see that they have a twofold symmetry and two axes of symmetry (one horizontal and one vertical). Now, by repeating such templates in a horizontal direction, we obtain patterns of pmm symmetry (Fig. 34b, c).

Conclusions
The RNM is a stable global root-finding process that needs a large number of iterations to achieve a root or a critical point of any polynomial. No chaos and no fractal structures are observable on the boundaries of basins when solving polynomial equations on the com- plex plane, contrary to the classical Newton's method.
In this paper, we modified the RNM by replacing the standard Picard iteration with the Mann iteration. The proposed M-RNM root-finding algorithm significantly improves the RNM algorithm. The M-RNM is much faster than the RNM. Its speed of convergence depends on the parameter α and the considered polynomial. The best value of α, which assures good convergence and performance, can be determined experimentally for a given polynomial. The M-RNM algorithm does not lose good properties of the RNM such as global convergence, stability, and robustness even for relatively large deviations from the best value of α. As the parameter α increases and moves far away from the best value, the patterns become increasingly dynamic and eventually become more and more ran-dom and chaotic. Also, for the proposed variable α sequence, we can observe the good behaviour of the M-RNM, obtaining better results than the RNM. The appearance of the polynomiographs generated with the M-RNM can be easily modified by using various colour palettes. However, the artistic aspects of the images can be additionally achieved with the help of specially selected sequences of α i as in Sect. 5. Furthermore, these images have been successfully used as templates for wallpaper symmetry or frieze symmetry patterns giving nicely looking examples of applied art.
Encouraged by the improvement of the speed of the RNM algorithm by using the Mann iteration in place of Picard's iteration, we are thinking about checking whether the use of other types of iterations known from the literature and collected, e.g. in [13], can lead to an even better improvement of the RNM algorithm. It also seems likely to find other intriguing patterns by using other types of iterations that may be considered aesthetic. Moreover, one can try to embed a proper template created with the M-RNM to the fundamental region of the spiral patterns introduced in [25]. These tasks will be the subject of our further research.
Funding This research received no external funding.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.