One more look on visualization of operation of a root-finding algorithm

Many algorithms that iteratively find solution of an equation require tuning. Due to the complex dependence of many algorithm’s elements, it is difficult to know their impact on the work of the algorithm. The article presents a simple root-finding algorithm with self-adaptation that requires tuning, similarly to evolutionary algorithms. Moreover, the use of various iteration processes instead of the standard Picard iteration is presented. In the algorithm’s analysis, visualizations of the dynamics were used. The conducted experiments and the discussion regarding their results allow to understand the influence of tuning on the proposed algorithm. The understanding of the tuning mechanisms can be helpful in using other evolutionary algorithms. Moreover, the presented visualizations show intriguing patterns of potential artistic applications.


Introduction
Methods for determining function's local minimum are based on the value of the function or its gradient (mostly numerically determined). The gradient method determines the direction of particle movement based on the knowledge of the gradient of the objective function (at the point reached in the previous step), and on this basis, it calculates the next position of the particle (Polak 1997)-a typical example of the gradient method is described by the expression (the gradient descent): where −∇ f (x i ) is the negative gradient of f , γ -step size, x i -the current position of the ith particle in a Ddimensional environment and x i -the previous position of the ith particle. Many modifications of this method are described in the literature, and their operation mechanisms  (Klein et al. 2009; Konečný and Richtárik 2017;Senov and Granichin 2017). A group of algorithms solving such problems includes evolutionary algorithms. The analysis of evolutionary algorithm is a very complex task (Gosciniak 2008(Gosciniak , 2017. Some optimization algorithms can have similar behavioural characteristics as evolutionary algorithms (Weise 2009). A particular attention should be paid to the particle swarm optimization (PSO) algorithms (Zhang et al. 2014). The complex nature of particles movement in these algorithms does not allow a precise definition of their effect on the algorithm.
Particle's movement in the PSO algorithm is described by the following equation (Yang and Li 2010): where x i is the current position of the ith particle in a Ddimensional environment, x i -the previous position of the ith particle and v i -the current velocity of the ith particle in a D-dimensional environment that is given by the following formula: where v i is the previous velocity of the ith particle, ω-the inertia weight (ω ∈ [0, 1]), η 1 , η 2 -acceleration constants (η 1 , η 2 ∈ (0, 1]), r 1 , r 2 -random numbers generated uni-formly in the [0, 1] interval, x pb i -the best position of the ith particle and x gb -the global best position of the particles. The inertia weight (ω) and acceleration constants (η 1 , η 2 ) play a very important role in the algorithm-they are responsible for the particle dynamics. The acceleration constants direct the particle towards its own best position or the global best position (Sengupta and Mishra 2014). The particle can be trapped in a local optima for too high values of the acceleration constants or cannot reach the solution for too low values of these constants. Balance of the exploration and exploitation is also controlled by the inertia weight (Bansal et al. 2011). The control rules of inertia weight can be classified as follows: constant (Shi and Eberhart 1998); random ; time varying (linear decreasing Shi and Eberhart 1999, sigmoid increasing/decreasing Malik et al. 2007, simulated annealing Al-Hassan et al. 2006, Sugeno function Lei et al. 2006, exponential decreasing law Chen et al. 2006;Li and Gao 2009, logarithmic decreasing law Gao et al. 2008); adaptive control using feedbacks of the optimization process (best fitness Saber et al. 2006;Shi and Eberhart 2001, fitness of the current and previous iterations Yang et al. 2007, global best and average local best fitness Arumugam and Rao 2008, particle rank Panigrahi et al. 2008, distance to particle and global best positions Qin et al. 2006 and distance to global best position Suresh et al. 2008). The value of inertia weight mostly presented in the literature is within the range of [0.4, 0.9] (see Jordehi and Jasni 2013). Both the inertia weight and the acceleration constants allow to control the particle behaviour within a wider range. The test function (problem being solved) and the selected method of iteration also influence particle's behaviour. The particle dynamics is determined by the adaptation mechanics resulting mostly from particles' cooperation. The method of algorithm parameters selecting is not deterministic-these values are selected by a tuning process depending on a solved problem and the iteration method that is used (similarly to evolutionary algorithms Bansal et al. 2011;Sengupta and Mishra 2014).
Due to the particle motion mechanism and particle's behaviour control parameters, we can use the wording-a particle dynamics. The tuning-selection of the proper values of the parameters-is a very important problem for the algorithm. In most cases, the algorithm tuning process is intuitive and it is based on user experience.
Finding of roots of a given function f , i.e. solving the equation f (x) = 0, is a very important problem in practical applications, e.g. physics (Franklin 2013), electronics (Chun and Kwasinski 2011) or computer graphics (Chen and Ma 2015). In many of these applications, the function is a polynomial one. So, the following question arises: can we solve the polynomial equation by its radicals (i.e. giving the formula for the roots in terms of the polynomial's coefficients, the four algebraic operations: addition, subtraction, multiplication, division and the extraction of roots, i.e. square roots, cube roots, etc.)? The answer to this question is that for polynomials of degree greater than four we cannot find such formulas (Grant and Kleiner 2015). The formulas for quadratic polynomial are known since the Babylonians. For the cubic polynomial the formulas were found by G. Cardano, whereas for the quartics by L. Ferrari both in the XVI century. In 1779 P. Ruffini and in 1826 N.-H. Abel proved the unsolvability by radicals of the polynomial equation of degree greater than four-the Abel-Ruffini theorem (Grant and Kleiner 2015). Thus, to solve a general polynomial equation we need to use some numerical methods, which are iterative methods and define discrete dynamical systems. In the literature, we can find many such methods, e.g. Newton's method (Kalantari 2009) and Halley's method (Kalantari 2009), or even the gradient descent method can be used for root finding (Kotarski and Lisowska 2018).
The analysis of dynamics of dynamical systems is a very important problem (Sayama 2015), because it gives us an understanding of the changes that occur in the system. To analyse the dynamics, we can use many different methods (Broer and Takens 2011). One of such methods is a graphical presentation of the dynamics. The graphical representation provides a lot of intuitive, geometrical insight into the system's dynamics, which would be hard to infer if we had just looked at algebraic equations (Sayama 2015). In the analysis of the complex polynomial root-finding process, the graphical presentation of dynamics is a very popular method and it even got its own name, namely polynomiography (Kalantari 2009). In polynomiography we visualize the root-finding process by using the number of iterations needed to obtain the root of a given polynomial and some colour map. The obtained images show the dynamics of the process. For the visualization and the analysis of algorithm's behaviour, we can use methods similar to these used in polynomiography (Kalantari 2004). The polynomiography visualizes the complex nature of the relationships between parameters.
The images generated using polynomiography very often reveal aesthetic patterns with a very complex structure (Gdawiec 2017). Because these complex patterns are generated with a simple method polynomiography, besides its analytic purposes, it also found application in computeraided design and in the arts (Kalantari 2004). Polynomiography assists the designer in the creation of aesthetic patterns, because the designer must select only some parameters used in the method and the computer generates the pattern, whereas without polynomiography to create the same pattern the designer must spend a lot of time and put a lot of work in the creation process.
The PSO-based gradient-like method is presented in the paper. In this simple root-finding algorithm, tuning parameters such as inertia weight and acceleration constant, and an adaptive mechanism depending on the location of the particle are proposed. This algorithm is used to show the influence of the tuning parameters on particle's behaviour. Moreover, we propose the use of different iteration processes instead of the standard Picard iteration that is used in the root-finding methods and optimization algorithms.
The aim of the article is to show the dynamics of particle motion and the influence of the parameters on particle's behaviour. For this purpose, we use a visualization method and analysis similar to the polynomiography. The used visualization method-presented in the article-can also have an artistic meaning. According to the authors, the understanding of the dynamics of particle's motion is a very important for its use in the group of evolutionary algorithms. The proper selection of particle dynamics may not cause the influence of the environment on the algorithm-as it is realized in the algorithm presented in Gosciniak (2017).
Let's try to summarize. What cannot be found in this article? The article is not a recipe how to effectively adjust the parameters of the algorithm. (Tuning an algorithm is a complex problem to be solved by other optimization algorithms-for instance by a genetic algorithm.) So what is this article about? Its main task is to visualize and analyse the complex nature of particle dynamics resulting from the selection of algorithm's coefficients and its interaction with the environment. However, if we analyse the behaviour of particles in evolutionary algorithms working in multidimensional spaces, a fractal analysis method can be proposed (Gosciniak 2017). What are the benefits of the article? Besides the mentioned visualization and the analysis of the complex nature of particle dynamics, the results of the paper can find application in the generative art. (Another example of this kind of algorithm can be the Newton's method, which is at the forefront in this field.) The rest of the paper is organized as follows. Section 2 introduces a root-finding algorithm that is based on the gradient method. This method will be used to illustrate the influence of the parameters on its behaviour. Next, Sect. 3 introduces the iteration processes known in fixed-point theory for finding the fixed points of different types of functions. Section 4 presents a method of visualization of algorithm's operation. And next, in Sect. 5 based on the obtained polynomiographs we discuss the research results. Finally, in Sect. 6 we give some concluding remarks.

The root-finding algorithm
Let f : R D → R be some function. We want to find zeroes of f , i.e. to solve the following equation To solve (4), we use the following algorithm: where x 0 ∈ R D is a starting position, v 0 = [0, 0, . . . , 0] is a starting velocity, v n+1 is the current velocity of the parti- . . . , v D n+1 ]) and x n is the previous position of the particle (x n = [x 1 n , x 2 n , . . . , x D n ]). The algorithm uses a similar methodology as the PSO algorithm, i.e. it sums the position of the particle x n with its current velocity v n+1 .
The current velocity of the particle is determined by the component of velocity of the previous iteration (inertia) and the component resulting from its current position (acceleration): where k ∈ {1, 2, . . . , D}, v k n+1 -the current velocity of the particle in the direction k, v k n -the previous velocity of the particle in the direction k, ω ∈ [0, 1)-inertia weight, η ∈ (0, 1]-acceleration constant, ε k -the step in the direc- Determining the particle motion towards the root is similar to the gradient method-the marked part of Eq. (6). The acceleration is influenced by the adaptation mechanism resulting from the value of the function in the previous position of the particlef (x n ). (It works effectively at determining a root of the function.) For this reason (when ω = 0) for x n fulfilling the condition f (x n ) = 0, the next particle position is x n+1 = x n . The algorithm's tuning involves the selecting inertia weight (ω) and acceleration constant (η)similarly as in the PSO algorithm. The effect of changes in these parameters, i.e. particle dynamics, will be visualized using the method described in Sect. 4.

Iteration processes
Let T : R D → R D be given by the following formula: Thus, algorithm (5) can be written in the following form In the literature, this type of iteration is called the Picard iteration and it is widely used in many optimization algorithms and in finding fixed points of a given function.
Let us notice that the Mann iteration for α n = 1 reduces to the Picard iteration, the Ishikawa iteration reduces to the Mann iteration when β n = 0 and to the Picard iteration when α n = 1, β n = 0, and the S-iteration reduces to the Picard iteration when α n = 0, or α n = 1 and β n = 0. A review of various iteration processes and their dependencies can be found in Gdawiec and Kotarski (2017).
All the presented iterations used only one mapping, but in the fixed-point theory exist iterations that use several mappings and are used to find common fixed points of the mappings. Let us recall some of them.
Let us notice that the Das-Debata iteration for T 1 = T 2 reduces to the Ishikawa iteration, the Khan-Cho-Abbas iteration reduces to the Agarwal iteration when T 1 = T 2 , and the generalized Agarwal iteration reduces to the Khan-Cho-Abbas iteration when T 1 = T 3 and to the Agarwal iteration when T 1 = T 2 = T 3 .
Having such a variety of different iteration processes, we can use them in our algorithm. In the iterations with a single mapping we use (7) as the mapping, and in the case of iteration with several mappings we use also (7), but with different values of ω and η parameters for each of the mappings.

Visualization of the dynamics
To visualize the dynamics of the proposed algorithm, we will use method that is very similar to the method used in polynomiography (Gdawiec et al. 2015). Because of this similarity, we will call the obtained images polynomiographs, although our algorithm finds roots of arbitrary functions and not only of polynomials as in polynomiography.
In the algorithm, we choose iteration method, e.g. one of the methods presented in Sect. 3, parameters ω, η for a single mapping T or ω 1 , ω 2 , ω 3 and η 1 , η 2 , η 3 for T 1 , T 2 , T 3 (depending on the chosen iteration). Moreover, we choose the maximum number of iterations m which the algorithm should realize, the accuracy of the computations r and a colouring function C : N → {0, 1, . . . , 255} 3 . Then, for each x 0 in the solution space A we use our algorithm. The iterations of the algorithm proceed till the convergence criterion: is satisfied or the maximum number of iterations is reached. When the algorithm ends, we assign a colour to x 0 using colouring function C and the number of performed iterations. The pseudocode of the visualization algorithm is presented in Algorithm 1. In the algorithm, the selected iteration method is denoted as I q , where q is a vector of parameters of the iteration.
The solution space A is contained in a D-dimensional space, so Algorithm 1 will return visualization in this space. If D = 2, then we obtain a single image that can be easily presented on a screen. When D > 2, then it is hard to visualize the result on a 2D surface of the screen. In this case, we make cross section of A with a two-dimensional plane and visualize this cross section.

Discussion on the research results
In this section, we present and discuss the results of visualizing the dynamics of the method introduced in Sect. 2 together with the various iteration methods presented in Sect. 3. In our study, we used functions f i : R 2 → R for i ∈ {1, 2, 3, 4} given by the following formulas: Algorithm 1: Visualization of the dynamics Input: f -function, A ⊂ R D -solution space, m -the maximum number of iterations, I q -iteration method, q ∈ [0, 1] N -parameters of the iteration I q , ω, ω 1 , ω 2 , ω 3 , η, η 1 , η 1 , η 2 , η 3 -parameters defining functions T , T 1 , T 2 , T 3 , C -colouring function, r -accuracy Output: visualization of the dynamics the function f 1 has three roots, namely the function f 2 has four roots, namely the function f 3 has five roots, namely Because R 2 is isomorphic with C and for each [x, y] ∈ R 2 we have x + iy = z ∈ C, then the functions f i for i ∈ {1, 2, 3, 4} can be written in the following form:

Fig. 1 Colour map used in the experiments
The functions (20)-(23) have the same roots as the following functions: Similar approach, to the one presented in this paper, for function z 3 −1 for which visualization of the dynamics using Mann and Ishikawa iteration in the Newton's method can be found in Kotarski et al. (2012).
In all experiments, the same colour map was used to colour the obtained images. The colour map is presented in Fig. 1. The colour map has 256 colours (ordered from left to right). The colour number represents the number of iterations needed to reach the solution. Moreover, in every experiment the following parameters were used: τ = 1.0e−3, m = 256, r = 1.0e−2, image resolution 800 × 800 pixels. The solution spaces depend on the function and are the following: The software used in the research was implemented in the C++ programming language. The experiments were realized on a computer with the Intel Core i5-2520M processor, 4 GB RAM, and Linux 3.16.7-42-desktop openSUSE 13.2 (Harlequin 64-bits, KDE Platform version 4.14.9).
In all experiments, the polynomiographs are generated using Algorithm 1 in 2D space. Thus, looking at the generated images we can read two important characteristics. The first one is the speed of convergence of the algorithm, i.e. the colour of each point gives us information on how many iterations were performed by the algorithm to reach the root. The second characteristic is the dynamics of the algorithm. Low dynamics is in areas where the variation of colours is small, whereas in areas with a large variation of colours the dynamics is high.

The Picard iteration
The behaviour of particles is affected by the speed and inertia. They are controlled by two parameters: the acceleration constant (η) and the inertia weight (ω) for the described above benchmark functions. The visualization of method's dynamics allows to analyse their impact on the algorithm's work.
We start the experiments with the f 1 function given by (16). This function has three roots that, thanks to the symme- Fig. 2 Polynomiographs of f 1 and the Picard iteration for ω = 0 and varying η try similar areas, create the corresponding dynamics of the particle. As was already mentioned, the parameters of the examples have been chosen so that in addition to the visualization of particle's dynamics, the polynomiograph also represents (according to the authors) an aesthetic image.
In Fig. 2 visualizations of the dynamics for the f 1 and the Picard iteration using ω = 0 and varying η are presented. The acceleration constant is the only factor influencing a particle motion-it affects its speed. (The larger the value of η, the greater the speed v n+1 .) Moreover, the generation time of each image is shown in Fig. 2. Looking at the results, we see that the increase in particle's velocity can shorten the creation time of the image. The black colour shows places which have not reached a solution in the given number of 256 iterations- Fig. 2a. The areas of the same colour indicate the same number of iterations required to determine the f (x n ) = 0-they look similar to the contour lines on the map (Fig. 2f). The magnifications of the marked areas in Fig. 2f are presented in Fig. 3. Inertia helps the particle to escape from a not promising area of solution space. Figure 4 presents examples of visualizations obtained for f 1 and the Picard iteration using rather high inertia, i.e. ω = 0.7 and varying η, and also their generation times. The shortest time of the creation of polynomiographs in Fig. 4 is obtained by the polynomiograph from Fig. 4c. Both too low and too high velocities may cause the increase in the time of polynomiograph creation.   Fig. 5c. In this case, both too low and too high inertia may cause the extension of the time. Figures 6 and 7 show examples of visualizations obtained for η = 0.1 and η = 0.3 (respectively) and various values of ω for function f 1 . Looking at the figures, we see that for the higher velocity of particle (acceleration constant increases in the following figures) the high value of inertia weight causes the increase in the creation time of the image.
From all the images obtained for f 1 and the Picard iteration, we see that the dynamics of the method changes in a significant way when we change two parameters of the method (ω, η). A smooth image represents low dynamics of the particles. Moreover, we see that using the proposed method we obtain interesting patterns and more diverse patterns are observed for particles with greater dynamics.
Polynomiographs of the Picard iteration for the f 2 function are shown in Fig. 8. Two groups of areas forming a different dynamics are clearly visible in Fig. 8a. The proper value of ω is responsible for creating a particle dynamics balance-it is evident in Fig. 8c. In this case, we observe a The next figure- Fig. 9-shows polynomiographs of the Picard iteration for the f 3 function. The function f 3 , thanks to the symmetry of roots placement, has four areas that create similar particle dynamics. The increase in the value of the acceleration constant (the particle moves faster-it increases in its dynamics) can shorten the time of creation of the polynomiograph. The dynamics of the particle is limited in the whole area of the polynomiograph-it is shown in Fig. 8b-d. However, areas from which the particle does not reach the solution (black colour) arise.
Two groups of areas with different dynamics are visible in Fig. 10. The images in this figure were created using the f 4 function and the Picard iteration. Increasing the value of the inertia weight increases the dynamics of the particle-it shortens the time of creation of the polynomiograph (Fig. 10a, b). The excessive value of inertia weight can be compensated by a decrease in the value of acceleration constant-consequently, similar particle dynamics is obtained in the whole area (Fig. 10c, d).

The Mann iteration
The Mann iteration is enhanced by the α parameter in relation to the Picard iteration. It gives the possibility to adjust the dynamic of the particle behaviour independently of the ω and η. Figure 11 presents the dynamics of the method for ω = 0, η = 0.025 using the Mann iteration and f 1 . For α = 1, we obtain the Picard iteration (Fig. 11a). From the images, we see that the decrease in the value of the parameter α, which is responsible for a linear combination in the Mann itera- tion, reduces the dynamics of particles. We can also observe that reducing the dynamics of particles can cause the lack of finding a solution (black colour), even for large areas of the solution space. Moreover, we see that the lower the value of α, the greater the generation time.
The acceleration constant (η) affects particle dynamics. Using a wide range of changes in α = 0.2, 0.5, 0.7, 0.9 (which limits the particle dynamics) in Fig. 12, there are no areas where the algorithm does not reach the solution. It is a result of increased particle dynamics by increasing the acceleration constant (η = 0.35).
Relatively large dynamics of particles for f 1 is presented in Figs. 13 and 14 obtained for different values of parameters ω and η defining T and the α parameter in the Mann iteration. From the images, we can observe that the decrease in the value of the parameter α, which reduces the proportion of the part specified by T in the Mann iteration, reduces the dynamic of the particle. The high particle dynamics allows to obtain interesting patterns.
Significantly lower particles dynamics is presented in Fig. 15. The polynomiographs in this figure were obtained for f 1 , ω = 0.3, η = 0.3 and varying α in the Mann iteration. These images are similar to the images obtained for the Picard iteration. The decrease in particle's dynamics smooths out the images-it is also visible in the remaining images. Figure 15a, b and Fig. 15c, d-although different in particle's dynamics-are similar in the colouring.
Tuning the algorithm is a multidimensional problem. It means that an increase in the number of parameters complicates this task. The polynomiographs from Fig. 16 (generated for the f 2 function) show the particle dynamics constituting the intersection of the algorithm tuning space for the α. The increase in the α parameter value increases the dynamics of the particles. Properly selected value of this parameter As in the previous example, in Fig. 17 polynomiographs generated using a varying α parameter are presented, but this time the f 3 function was used. A limitation of the dynamic range observed in Fig. 17c occurs for the shortest time of a polynomiograph creation. Moreover, we can observe that the change in the α parameter by 0.1 affects the dynamics in different way. In Fig. 17a, b the dynamics is almost the   eter is small, which significantly limits the dynamics of the particle. The increase in the value of the acceleration constant slightly reduces the time of the polynomiograph creation- Fig. 18b. A significant reduction in the inertia weight causes a reduction in particle dynamics and reduces the time needed to the polynomiograph creation-it is visible in Fig. 18d.

The Ishikawa and the Das-Debata iterations
In the Ishikawa iteration, an additional parameter β is introduced in relation to the Mann iteration. This parameter, like the parameter α, has an impact on the dynamics of particles. It is also involved in the creation of an additional reference point (the sample in the solution space). In this way, particle motion is a linear combination of its position and a movement of reference point. This gives the possibility to define two different operators (T 1 and T 2 ) for transforming the current position of the particle and the position of the reference point-the Das-Debata iteration. For the differentiation of these operators, the parameters ω 1 , ω 2 and η 1 and η 2 are responsible.
Two types of iterations-Ishikawa and Das-Debata-give a wide range of possibilities for the regulation of the algorithm's work. The use of these iterations causes twisting contour lines. The α parameter, as in the case of the Mann iteration, limits the dynamics of a particle. A similar function as the parameter α is realized by the parameter β in relation to the reference point.  Figures 19 and 20 show the reduction of the dynamics of the particle for f 1 by the α parameter in the Ishikawa iteration. A small value of β parameter causes that the reference point is created close to the active point. The same value of the α parameter is the cause of creation of similar patterns.
The reduction of the dynamics of the reference point creation for f 1 by the α parameter in the Ishikawa iteration is shown in Figs. 21 and 22. A significant decrease in the dynamics of the reference point creation weakens the effect of twisting lines-it is also visible in the direction to the centre of the pattern. For large values of the β parameter, major changes occur on the periphery of the image- Fig. 22. Polynomiographs presented in Fig. 23 for the f 2 function and the Ishikawa iteration show particle dynamics for changing ω and α parameters. Comparing Fig. 23a, b, an increase in dynamics due to an increase in the ω value is observed, which shortens the time of the polynomiograph creation. Then a further increase in dynamics due to an increase in the ω value results in a longer time of a polynomiograph creation- Fig.  23c. Reducing the dynamics of the particle through the appropriate selection of the α parameter reduces the time of the polynomiograph creation. This method of particle dynamics reduction is ineffective and may cause the particle not to reach a solution.  Figure 24 shows the changes in particle dynamics for the function f 3 under the influence of changes in the α and β parameters for the Ishikawa iteration. Comparing the times of polynomiographs creation in Fig. 24b, c with Fig. 24a or d, the effect of change in the α parameter on particle's dynamics is noticeably greater than that of β.
Polynomiographs obtained for the f 4 function using the Ishikawa iteration are presented in Fig. 25. The obtained results confirm the observations from the previous example (Fig. 24), i.e. the change in the α parameter has a much greater effect on the behaviour of the particle than the change in the β parameter. Figures 26,27,28,29,30,31,32,33,34 and 35 present polynomiographs created using the Das-Debata iteration for the f 1 function. In this iteration, the ω 1 inertia weight is responsible for the dynamics of the reference point creation. Small differences of the dynamics of the reference point creation can cause the creation of similar patterns-Figs. 26 and 27. A similar effect is given by the change in the reference point dynamics by determining the η 1 parameter-see Dynamics of particle's movement has a significant impact on the appearance of the polynomiograph. The parameters ω 2 and η 2 are responsible for this dynamics. As it was already mentioned, both the dynamics of the reference point and the dynamics of the particle are responsible for the creation of intriguing images. Figures 30 and 31 show the effect of the change in inertia weight ω 2 on the particle dynamics. High particle dynamics allows to obtain interesting images. Figures 32 and 33 show the effect of the change in the acceleration constant η 2 on the creation of a polynomiograph. In these images, a smaller dynamics is observed.
The images in Figs. 34 and 35 visualize a high dynamics of particle's movement. Due to varying of a large number of parameters, the detailed analysis of their dependence is difficult. The general rule that proper selection of the dynamics of particles reduces the operating time of the algorithm is confirmed. Moreover, from the images we see that the use of Das-Debata iteration allows to create interesting patterns.
Polynomiographs of the f 2 function and the Das-Debata iteration for varying all parameters are presented in Fig. 36. The obtained polynomiographs present a high dynamics of the particle. The shortest time of polynomiograph's creation is obtained for Fig. 36c. In this case, the narrowed range of the particle's dynamic for the entire polynomiograph area is observed. As it was already mentioned, tuning the algorithm cannot be covered by the deterministic rule due to the very complex nature of the relationship between the algorithm's coefficients and the environmental impact.  Figure 37 presents polynomiographs obtained for the f 3 function using the Das-Debata iteration and varying ω 1 parameter. In subsequent images, it is possible to observe the effect of increasing in the value of ω 1 , which is responsible for the dynamics of creating the reference point. Changes in the values of inertia weight of the reference point transformation have a small impact on the dynamics of the particle, as well as changes in the β parameter. Appropriate selection of the ω 1 parameter can shorten the time of the polynomiograph's generation-it is presented in Fig. 37c.
The results of using the Das-Debata iteration for f 4 and a varying η 2 parameter are presented in Fig. 38. Even for low values of the acceleration constant parameter, but with high dynamics of the reference point creation and a high inertia weight of the particle, it is possible to obtain a sufficient particle dynamics for an effective movement in the solution space. The obtained polynomiographs present very similar particle dynamics. The shortest time of polynomiograph's creation is obtained for Fig. 38b.

The Agarwal and the Khan-Cho-Abbas iterations
The generalized Agarwal iteration is more complex than the previously presented iterations. It introduces an additional transformation T 3 (with parameters ω 3 , η 3 ) of particle in relation to the Das-Debata iteration. Figures 39 and 40 present images obtained with the Agarwal iteration (T 1 = T 2 = T 3 ) for f 1 and varying α and β parameters. When we look at these images, we see that they have features of a twisting line, similarly to the Ishikawa iteration. Depending on the dynamics of the particles, it is possible to obtain images similar to the ones obtained with the Ishikawa iteration.
Operator T 3 causes the introduction of an additional dynamics of particles, the contribution of which produces the unique visual effects-the additional twisting lines are visible. The increase in dynamics of particles allows creating images of an artistic meaning. Moreover, we observe that the time of generation of the images (the working time of the algorithm) is longer that in the case of the previous iterations. This is due to the complexity of the iteration and the particle's dynamics.
Polynomiographs of the Agarwal iteration for the f 3 function are shown in Fig. 42. They present the limitation of particle dynamics by reducing the values of the inertia weights ω 1 , ω 2 and ω 3 -limiting the dynamics results in longer polynomiograph creation time.
The last example for the Agarwal iteration- Fig. 43presents images obtained for the f 4 function. The obtained images show the particle dynamics increase by the increasing The generalized Agarwal iteration fulfils the condition T 1 = T 3 . Operators T 2 and T 3 are the same for Figs. 44 and 45, which were obtained for the f 1 function. The dynamics of the particle is influenced by three transformations. Even large change in the transformation T 1 (by using ω 1 or η 1 ), which Polynomiographs of the Khan-Cho-Abbas iteration for fixed α, β and varying ω 1 , η 1 , ω 2 , η 2 , ω 3 , η 3 for the f 2 function are presented in Fig. 48. The ω and η parameters change in the opposite directions, i.e. when ω increases/decreases, then η decreases/increases. It results in a characteristic change in dynamics presented in the polynomiographs by changing the colours. The shortest time of polynomiograph's creation is obtained for Fig. 48b.
A similar change in the parameters is introduced in the examples of using the Khan-Cho-Abbas iteration that are shown in Figs. 49 and 50. The images in Fig. 49 are obtained

Fig. 52
Polynomiographs of f 1 and the generalized Agarwal iteration for α = 0.8, β = 0.3, ω 1 = 0.2, η 1 = 0.5, ω 2 = 0.6, η 2 = 0.7, ω 3 = 0.5 and varying η 3 for the f 3 function, whereas the ones in Fig. 50 for the f 4 function. Shortening the time of the polynomiograph, creation for Figs. 49b and 50b is obtained-these images visualize the reduction of the particle dynamic range. In all cases, the influence of the reference point on the movement of the particle is limited.
The generalized Agarwal iteration introduces T 3 transformation with ω 3 and η 3 parameters. Figures 51 and 52 show the effect of changing in these parameters on the creation of polynomiographs for the f 1 function. Due to the adaptation mechanism on which the particle acceleration depends, the for ω 1 = 0.9 and varying α, β, η 1 , ω 2 , η 2 , ω 3 , η 3 inertia of the particle has a significantly greater effect on the visual effects. This can be also observed in other presented images.
Moreover, the images in Figs. 53 and 54 were obtained for the f 1 function using the generalized Agarwal iteration and varying different parameters. The shortest times are obtained for images in Figs. 53b and 54a. These images illustrate the high particle dynamics, and they look like brush painting images.
Polynomiographs of the generalized Agarwal iteration for the f 2 function and varying all parameters are presented in Fig. 55. Shortening the time of polynomiographs creation occurs in subsequent images. (The shortest time is obtained for Fig. 55d.) The generalized Agarwal iteration gives the greatest possibility of the particle dynamics control.
The last two figures present polynomiographs generated by the generalized Agarwal iteration for the f 3 (Fig. 56) and f 4 (Fig. 57) functions. The generation algorithm is tuned just like in the previous example. We observe a limitation of the Fig. 55 Polynomiographs of the generalized Agarwal iteration for varying α, β, ω 1 , η 1 , ω 2 , η 2 , ω 3 , η 3

Conclusions
The inertia weight and the acceleration constants are presented in many optimization algorithms which require parameters tuning-it is exemplified by the PSO algorithm. The discussion presented in the article allows to show their impact on the work of the proposed root-finding algorithm. Both too large and too small dynamics of the particles will have an adverse effect on the work of the algorithm. Moreover, the presented experiments show that the dependence of particle's dynamics on the parameters used in the root-finding method (inertia weight and acceleration constant) and in the various iteration methods is a non-trivial function.
The visualization of the dynamics is a tool that not only illustrates the behaviour of particles, but also creates attractive patterns of artistic value. Thus, these patterns and the method of their creation can be used in graphics design. Some of the presented images look like brush paintings. Particle dynamics reflects the dynamics of the brush movement. This allows to obtain images that can be determined as art.
The genetic algorithm can be used to tune the algorithm to obtain optimal behaviour due to the large number of parameters. This approach can be realized as a future work.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflicts of interest.
Human participants or animals This article does not contain any studies with human participants or animals performed by the authors.
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://creativecomm ons.org/licenses/by/4.0/.