# On root-ratio multipoint methods for finding multiple zeros of univariate
functions

Miodrag S. Petković 1,, Ljiljana D. Petković 2,111Corresponding author
1Faculty of Electronic Engineering, University of Niš,
A. Medvedeva 14, 18000 Niš, Serbia
2Faculty of Mechanical Engineering, University of Niš,
A. Medvedeva 14, 18000 Niš, Serbia

###### Abstract

Several root-ratio multipoint methods for finding multiple zeros of univariate
functions were recently presented. The characteristic of these methods is that
they deal with $m$-th root of ratio of two functions (hence the name root-
ratio methods), where $m$ is the multiplicity of the sought zero, known in
advance. Some of these methods were presented without any derivation and
motivation, it could be said, out of the blue. In this paper we present an
easy and entirely natural way for constructing root-ratio multipoint iterative
methods starting from multipoint methods for finding simple zeros. In this
way, a vast number of root-ratio multipoint methods for multiple zeros,
existing as well new ones, can be constructed. For demonstration, we derive
four root-ratio methods for multiple zeros. Besides, we study computational
cost of the considered methods and give a comparative analysis that involves
CPU time needed for the extraction of the $m$-th root. This analysis shows
that root-ratio methods are pretty inefficient from the computational point of
view and, thus, not suitable in practice. A general discussion on a practical
need for multipoint methods of very high order is also considered.
AMS Mathematical Subject Classification (2010): 65H05.
Keywords: Solving nonlinear equations; Multipoint methods; Multiple zeros;
Root-ratio approach; Computational efficiency.

††E-mail address<EMAIL_ADDRESS>(L. D. Petković)

##  1 Introduction

In the last ten years a lot of papers were published in the topic of iterative
methods of optimal order for finding multiple zeros of univariate functions.
Among them, several methods were constructed using root-parameter approach
dealing with parameters computed by extracting $m$-th root of some values. In
this paper this class of methods will be called root-ratio methods. We
concentrate on the two issuses: (i) computational cost of these methods and
(ii) demonstration of an easy procedure for natural construction of root-ratio
methods starting from methods for finding simple zeros.

The paper is organized as follows. In Section 2 we present three two-point
root-ratio methods of optimal order four published in [2], [3] and [4]. In
Section 3 we expose three tree-point root-ratio methods of optimal order eight
proposed in [5], [6] and [7]. An easy and fully natural procedure for
constructing root-ratio methods of optimal order, demonstrated by three
examples, is presented in Section 4. The computational efficiency of root-
ratio methods is given in Section 5 applying a comparative analysis that
involves CPU time needed for the extraction of the $m$-th root. This analysis
has shown that root-ratio methods are pretty expensive and not-competitive
with multipoint methods which do not require the extraction of the root even
if the latter ones have lower order of convergence. See Remark 8 at the end of
the paper.

##  2 Two-point fourth order methods

Let $\alpha$ be a zero of the known multiplicity $m\geq 1$ of a differentiable
function $f,$ and $k=0,1,\ldots$ be the iteration index. We will assume that
the chosen initial approximation $x_{0}$ is sufficiently close to the zero
$\alpha$ so that the considered iterative processes are convergent. We start
with two-point root-ratio methods for finding multiple zero of functions of
one variable. Zhou et al. proposed in [2] the following family of iterative
two-step methods for finding multiple zeros:

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=y_{k}-mG(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\quad
u_{k}=\Bigl{(}\dfrac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m}\end{array}\right.\quad(k=0,1,\ldots).$
(1)

Under the conditions $G(0)=0,\ G^{\prime}(0)=1,\ G^{\prime\prime}(0)=4,\
G^{\prime\prime\prime}(0)<+\infty,$ the order of convergence of the methods
(1) is four. The iterative formula was derived correctly but in a complicated
way.

Lee et al. constructed in [3] the following family of fourth order methods:

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\cdot\dfrac{f(x_{k})}{f^{\prime}(x_{k})+\lambda
f(x_{k})},\\\\[13.0pt]
x_{k+1}=y_{k}-mW(u_{k})\cdot\dfrac{f(x_{k})}{f^{\prime}(x_{n})+2\lambda
f(x_{k})},\quad
u_{k}=\Bigl{(}\dfrac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m},\end{array}\right.$ (2)

where

$W(u)=\frac{u(1+(c+2)u+ru^{2})}{1+cu}$

and $\lambda,\ c,\ r$ are arbitrary parameters.

Remark 1. The parameter $\lambda$ in (2) does not improve convergence
characteristics of the methods (2). Indeed, if $m\geq 2$, then both values
$f^{\prime}(x_{k})$ and $f(x_{k})$ are very small in magnitude if the
approximation $x_{k}$ is sufficiently close to the zero $\alpha$ of $f$ and,
consequently, the parameter $\lambda$ does not play any important role so that
one can take $\lambda=0$ in (2) without loss of generality. The derivatives of
$W$ are given by

$W^{\prime}(u)=\frac{1+2(c+2)u+(2c+c^{2}+3r)u^{2}+2cru^{3}}{(1+cu)^{2}},\quad
W^{\prime\prime}(u)=\frac{4+2ru(3+3cu+c^{2}u^{2})}{(1+cu)^{3}}.$

Setting $u=0$ one obtains $W(0)=0,\ \ W^{\prime}(0)=1,\ \
W^{\prime\prime}(0)=4.$ Therefore, the method (2) is essentially a special
case of the more general method (1).

Using root-ratio approach in a slightly different way, Liu et al. derived in
[4] the family of two-point iterative methods

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=y_{k}-mG(t_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\quad
t_{k}=\Bigl{(}\dfrac{f^{\prime}(y_{k})}{f^{\prime}(x_{k})}\Bigr{)}^{\frac{1}{m-1}}\end{array}\right.\quad(k=0,1,\ldots),$
(3)

which converges with the order of convergence four under the conditions

$G(0)=0,\quad G^{\prime}(0)=1,\quad G^{\prime\prime}(0)=\frac{4m}{m-1}.$

The family (3) cannot be applied for finding simple zeros ($m=1$), which is an
obvious disadvantage. This drawback limits its applications in composite
algorithms where a zero-finding method is part of the algorithm in which the
case $m=1$ is possible; actually, this case often appears in practical
problems (for example, solving engineering problems).

All three two-point methods (1), (2) and (3) are optimal in the sense of Kung-
Traub conjecture [8]: The highest order of convergence of an $n$-point method
requiring $n+1$ function evaluations is $2^{n}$.

##  3 Three-point methods of eighth order

In this section we present some recent three-point optimal iterative methods
for finding multiple zeros of the known multiplicity.

Behl et al. proposed in [5] (without derivation and motivation) the family of
three-point methods of optimal order eight:

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-u_{k}S(h_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-u_{k}v_{k}R(h_{k},v_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(4)

where

$u_{k}=\Bigl{(}\frac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m},\quad
h_{k}=\frac{u_{k}}{a_{1}+a_{2}u_{k}}\ (a_{1}\neq 0),\quad
v_{k}=\Bigl{(}\frac{f(z_{k})}{f(y_{k})}\Bigr{)}^{1/m}.$

Under some specific conditions for the function $S$ (at the point 0) and $R$
(at the point (0,0)), the order of convergence of the family of iterative
methods (4) is eight.

Zafar et al. presented in [6] (without derivation and motivation) the
following iterative three-point method

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-mu_{k}H(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-u_{k}v_{k}(B_{1}+B_{2}u_{k})P(v_{k})G(w_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(5)

where $B_{1},B_{2}\in\mbox{\boldmath{$R$}}$ are free parameters, $H,\ P,\ G$
are the weight functions and

$u_{k}=\Bigl{(}\frac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m},\quad
v_{k}=\Bigl{(}\frac{f(z_{k})}{f(y_{k})}\Bigr{)}^{1/m},\quad
w_{k}=\Bigl{(}\frac{f(z_{k})}{f(x_{k})}\Bigr{)}^{1/m}.$ (6)

Zafar et al. proposed in [7] (without derivation and motivation) the family of
three-point methods

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-mu_{k}H(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-u_{k}P(u_{k})G(v_{k})L(w_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(7)

where $u,\ v,\ w$ are defined by (6). The family (7) is obvious generalization
of the family (5). Indeed, comparing (5) and (7), we note that
$P(u)=B_{1}+B_{2}u,\ G(v)=vP(v).$ Clearly, the authors could work immediately
at start with the product of weight functions $P(u)G(v)L(w)$ since the
generalization is quite obvious.

Remark 2. Neither the idea nor the motivation for the construction of the
methods (4), (5) and (7) were presented; the iterative formulas appear out of
the blue, without preliminary explanation/introduction of the basic idea and
derivation procedures. Such non-preamble approach is not in the spirit of the
methodology of scientific work and educational principles; hence it is of a
little significance for readers.

Remark 3. The need for arbitrary parameters in any iterative formula is
discussed in Remark 7. This relates to the parameters $\lambda,c,r$ in (2),
and especially to $a_{1}$ and $a_{2}$ in (4). Since any advantage of using
$a_{1}$ and $a_{2}$ was not proved in [5], it is logical to take $a_{1}=1$ and
$a_{2}=0$ without loss of generality. At first sight, this is a special case
but, in fact, it is a natural choice which simplifies the iterative formula
(4). Consequently, it follows $h_{k}=u_{k}$ and the family (4) reduces to a
simpler iterative method (without $h_{k}$)

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-u_{k}S(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-u_{k}v_{k}R(u_{k},v_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots).$
(8)

##  4 Can we construct them easier?

In this section we show that some of the above-presented methods can be
derived using an easy way, well known in the literature for almost 150 years.
The idea was known to German mathematician E. Schröder, see his paper [9] and
the English translation [10], and goes as follows in the case of Newton’s
method

$x_{k+1}=x_{k}-\frac{f(x_{k})}{f^{\prime}(x_{k})}\quad(k=0,1,\ldots)$

applied to the function $f$ having a zero $\alpha$ of multiplicity $m\geq 1.$

4.1 Schröder’s $f^{1/m}$-approach

Let $F(x)=f(x)^{1/m}.$ Then $\alpha$ is a simple zero of the function $F.$
Applying Newton’s method to $F$, we obtain

$x_{k+1}=x_{k}-\frac{F(x_{k})}{F^{\prime}(x_{k})}=x_{k}-\frac{f(x)^{1/m}}{\frac{1}{m}f^{\prime}(x)f(x)^{1/m-1}}$

and hence

$x_{k+1}=x_{k}-m\frac{f(x_{k})}{f^{\prime}(x_{k})}\quad(k=0,1,\ldots),$

which is well-known Schröder’s method for finding multiple zero of the known
multiplicity $m.$ This useful idea was applied in many papers during the last
70 years.

Schröder’s approach can be applied to multipoint methods for finding a simple
zero in order to construct corresponding multipoint mehods for multiple zeros.
As above, let $f$ be a function having a zero of the known multiplicity $m\geq
1.$ First, we set $F(x)=f(x)^{1/m}$ (Schröder’s $f^{1/m}$-approach) and derive
the following relations:

$\displaystyle\frac{F(x)}{F^{\prime}(x)}$ $\displaystyle=$
$\displaystyle\frac{f(x)^{1/m}}{\frac{1}{m}f^{\prime}(x)f(x)^{1/m-1}}=m\frac{f(x)}{f^{\prime}(x)},$
$\displaystyle\frac{F(y)}{F^{\prime}(x)}$ $\displaystyle=$
$\displaystyle\frac{f(y)^{1/m}}{\frac{1}{m}f^{\prime}(x)f(x)^{1/m-1}}=m\frac{f(x)}{f^{\prime}(x)}\Bigl{(}\frac{f(y)}{f(x)}\Bigr{)}^{1/m}=mu\frac{f(x)}{f^{\prime}(x)},\quad
u=\Bigl{(}\frac{f(y)}{f(x)}\Bigr{)}^{1/m},$
$\displaystyle\frac{F(z)}{F^{\prime}(x)}$ $\displaystyle=$
$\displaystyle\frac{f(z)^{1/m}\cdot\frac{f(y)^{1/m}}{f(y)^{1/m}}}{\frac{1}{m}f^{\prime}(x)f(x)^{1/m-1}}=m\frac{f(x)}{f^{\prime}(x)}\Bigl{(}\frac{f(y)}{f(x)}\Bigr{)}^{1/m}\Bigl{(}\frac{f(z)}{f(y)}\Bigr{)}^{1/m}=muv\frac{f(x)}{f^{\prime}(x)},\quad
v=\Bigl{(}\frac{f(z)}{f(y)}\Bigr{)}^{1/m},$ $\displaystyle\frac{F(y)}{F(x)}$
$\displaystyle=$ $\displaystyle\frac{f(y)^{1/m}}{f(x)^{1/m}}=u,$ (9)
$\displaystyle\frac{F(z)}{F(y)}$ $\displaystyle=$
$\displaystyle\frac{f(z)^{1/m}}{f(y)^{1/m}}=v,$
$\displaystyle\frac{F(z)}{F(x)}$ $\displaystyle=$
$\displaystyle\frac{f(z)^{1/m}}{f(x)^{1/m}}=w=uv.$

4.2 Eighth-order family (4) derived by Schröder’s approach

For demonstration, we will apply the presented Schröder’s approach to the
following three-point family of iterative methods for finding a simple zero of
a function $f,$ proposed in [20]:

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-p(u_{k})\dfrac{f(y_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-q(u_{k},v_{k})\dfrac{f(z_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(10)

where

$u_{k}=\frac{f(y_{k})}{f(x_{k})},\quad v_{k}=\frac{f(z_{k})}{f(y_{k})}$

and $p(u)$ and $q(u,v)$ are functions of one and two variables, respectively.
The following theorem has been proved in [20] (see, also, [12]).

Theorem 1. Let $a,\ b$ and $c$ be arbitrary constants. If $p$ and $q$ are
arbitrary differentiable functions with Taylor’s series of the form

$\displaystyle p(u)$ $\displaystyle=$ $\displaystyle
1+2u+\frac{a}{2}u^{2}+\frac{b}{6}u^{3}+\cdots,$ $\displaystyle q(u,v)$
$\displaystyle=$ $\displaystyle
1+2u+v+\frac{2+a}{2}u^{2}+4uv+\frac{c}{2}v^{2}+\frac{6a+b-24}{6}u^{3}+\cdots,$

then the family of three-point methods (10) is of order eight. It is assumed
that higher-order terms are represented by the dots, and they can take
arbitrary values.

Assume that $f$ has the zero $\alpha$ of multiplicity $m\geq 1,$ known in
advance. We modify (10) by using the above-mentioned Schröder’s strategy with
the function $F(x)=f(x)^{1/m}$ for which $\alpha$ is a simple zero. First,
replacing $f$ with $F$, rewrite (10) in the form

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\dfrac{F(x_{k})}{F^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-P(u_{k})\dfrac{F(y_{k})}{F^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-Q(u_{k},v_{k})\dfrac{F(z_{k})}{F^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(11)

where

$u_{k}=\frac{F(y_{k})}{F(x_{k})},\quad v_{k}=\frac{F(z_{k})}{F(y_{k})}.$

Using the relations (9), the family (11) can be written in the following form:

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-mu_{k}P(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}={\varPhi}(x_{k}):=z_{k}-mu_{k}v_{k}Q(u_{k},v_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(12)

where

$u_{k}=\Bigl{(}\frac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m},\quad
v_{k}=\Bigl{(}\frac{f(z_{k})}{f(y_{k})}\Bigr{)}^{1/m}.\quad$ (13)

The weight functions $P(u)$ and $Q(u,v)$ should be determined in such a way
that the iterative methods defined by (12) have the maximal order of
convergence using only the calculation of $f(x),\;f^{\prime}(x),\;f(y)$ and
$f(z).$ In the ideal case, the order would reach eight, and the family (12)
would be optimal according to Kung-Traub’s hypothesis [8].

We proceed using a standard technique based on Taylor’s series by employing
symbolic computation in computer algebra system Mathematica. Since this
technique was seen many times in existing papers, we present only an outline
of convergence analysis.

We omit the iteration index $k$ and define the errors

$\varepsilon=x-\alpha,\ \ \varepsilon_{y}=y-\alpha,\ \
\varepsilon_{z}=z-\alpha,\ \ \hat{\varepsilon}=\hat{x}-\alpha,$

where $\hat{x}$ is a new approximation $x_{k+1}.$ Introduce

$C_{r}=\frac{m!}{(m+r)!}\frac{f^{(m+r)}(\alpha)}{f^{(m)}(\alpha)}\ \
(r=1,2,\ldots).$

We will use the following development of the function $f$ about the zero
$\alpha$ of multiplicity $m$

$f(x)=\frac{f^{(m)}(\alpha)}{m!}\varepsilon^{m}\Bigl{(}1+C_{1}\varepsilon+C_{2}\varepsilon^{2}+C_{3}\varepsilon^{3}+C_{4}\varepsilon^{4}+C_{5}\varepsilon^{5}+C_{6}\varepsilon^{6}+C_{7}\varepsilon^{7}+C_{8}\varepsilon^{8}+O(\varepsilon^{9})\Bigr{)},$

and a program in Mathematica. As usual, in finding the weight functions $P$
and $Q,$ we represent these functions by their Taylor’s series at the
neighborhood of $u=0$ (for $P$) and $(u,v)=(0,0)$ (for $Q$):

$\displaystyle
P(u)=P(0)+P^{\prime}(0)u+\dfrac{P^{\prime\prime}(0)}{2}u^{2}+\cdots,$
$\displaystyle
Q(u,v)=Q(0,0)+Q_{u0}(0,0)u+Q_{0v}(0,0)v+\frac{1}{2!}\Bigl{(}Q_{uu}(0,0)u^{2}+2Q_{uv}(0,0)uv+Q_{vv}(0,0)v^{2}\Bigr{)}+\cdots$

Subscript indices denote partial derivatives. In the last developments, as
well as in the program, the following notation is used:

$\displaystyle fam=f^{(m)}(\alpha),\ fx=f(x),\ fy=f(y),\ fz=f(z),\
f1x=f^{\prime}(x),$ $\displaystyle e=\varepsilon,\ ey=\varepsilon_{y},\
ez=\varepsilon_{z},\ e1=\hat{\varepsilon},$ $\displaystyle
Pr=\frac{P^{(r)}(u)}{du}\Big{|}_{(u=0)}\ (r=0,1,2),$ $\displaystyle
Q00=Q(0,0),\ Qu0=\frac{\partial\,Q}{\partial\,u}\Bigl{|}_{(u,v)=(0,0)},\
Q0v=\frac{\partial\,Q}{\partial\,v}\Bigl{|}_{(u,v)=(0,0)},$ $\displaystyle
Quu=\frac{\partial^{\,2}Q}{\partial\,u^{2}}\Bigl{|}_{(u,v)=(0,0)},\
Quv=\frac{\partial^{\,2}Q}{\partial\,u\,\partial\,v}\Bigl{|}_{(u,v)=(0,0)},\
Qvv=\frac{\partial^{\,2}Q}{\partial\,v^{2}}\Bigl{|}_{(u,v)=(0,0)}.\ $

The coefficients of Taylor’s developments of the weight functions $P$ and $Q$
are determined using an interactive approach by combining the program realized
in Mathematica (two parts) and the annihilation of coefficients standing at
$\varepsilon$ of lower power.

PART I (Mathematica)

fxx=1+C1*e+C2*e^2+C3*e^3+ C4*e^4+C5*e^5+C6*e^6+ C7*e^7+C8*e^8;
fx=fam/m!*e^m*fxx; fx1=D[fx,e]; newt=Series[fx/fx1,{e,0,8}];
ey=e-m*newt; fyy=1+C1*ey+C2*ey^2+3*ey^3+C4*ey^4;
yx=fyy*Series[1/fxx,{e,0,8}];
yxm=Series[yx^1/m,{e,0,8}];
u=yxm*ey/e;
P=P0+P1*u+P2/2*u^2+P3/6*u^3+P4/24*u^4;
ez=Series[ey-m*u*P*newt//FullSimplify,{e,0,8}]

This programm gives the following OUTCOME:

$\displaystyle\varepsilon_{z}$ $\displaystyle=$
$\displaystyle\sum_{r=2}^{8}S_{r}\varepsilon^{r}$ $\displaystyle=$
$\displaystyle\frac{(C_{1}-C_{1}P_{0})\varepsilon^{2}}{m}$
$\displaystyle+\frac{(-2C_{2}m(-1+P_{0})+C_{1}^{2}(-1+m(-1+P_{0})+3P_{0}-P_{1}))\varepsilon^{3}}{m^{2}}$
$\displaystyle+\frac{1}{2m^{3}}\Bigl{(}-6C_{3}m^{2}(-1+P_{0})+2C_{1}C_{2}m(-4+3m(-1+P_{0})+11P_{0}-4P_{1})$
$\displaystyle+C_{1}^{3}(2-13P_{0}+10P_{1}+m(4-2m(-1+P_{0})-11P_{0}+4P_{1})-P_{2})\Bigr{)}\varepsilon^{4}+\sum_{r=5}^{8}U_{r}+O\bigl{(}e^{9}\bigr{)}.$

To annihilate coefficients by $\varepsilon^{2}$ and $\varepsilon^{3}$, from
the condition $S_{2}=0$ and $S_{3}=0$ we find

$P_{0}=1,\ \ P_{1}=2,\ \ P_{2}\ \mbox{\rm arbitrary},$ (14)

yielding

$\varepsilon_{z}=\frac{-2mC_{1}C_{2}+C_{1}^{3}(9+m-P_{2})}{2m^{3}}\varepsilon^{4}+O\bigl{(}\varepsilon^{5}\bigr{)}.$
(15)

The part II of the program uses previously found entries and serves for
finding additional conditions which provide optimal order eight.

PART II - CONTINUATION (Mathematica)

fzz=1+C1*ez+C2*ez^2; fyy=1+C1*ey+C2*ey^2+ C3*ey^3+C4*ey^4;
zy=fzz*Series[1/fyy{e,0,8}]; zym=Series[zy^(1/m),{e,0,8}];
v=zym*ez/ey;
Q=Q00+Qu0*u+Q0v*v+Quu/2*u^2+Qvv/2*v^2+Quv*u*v;
e1=Series[ez-m*u*v*Q*newt,{e,0,8}]//FullSimplify

The error $\hat{\varepsilon}=\hat{x}-\alpha\ (=e1)$ is given in the form

$\varepsilon_{1}=\sum_{r=4}^{8}T_{r}\varepsilon^{r}+O\bigl{(}\varepsilon^{9}\bigr{)}$

From the conditions $T_{4}=0,\ T_{5}=0,\ T_{6}=0,\ T_{7}=0,$ we find the
coefficients

$Q_{00}=1,\ \,Q_{u0}=2,\ \,Q_{0v}=1,\ \,Q_{uv}=4,\ \,P_{2}\ \mbox{\rm is
arbitrary},\ \,P_{3}=24-6P_{2},\ \,Q_{uu}=P_{2}+2.$ (16)

In addition, we obtain

$\displaystyle T_{8}$ $\displaystyle=$
$\displaystyle-\frac{1}{48m^{7}}\Bigl{\\{}C_{1}\Bigr{[}C_{1}^{2}(m-P_{2}+9)-2mC_{2}\Bigr{]}$
$\displaystyle\times\Bigl{[}C_{1}^{4}\left(-14m^{2}+3Q_{tt}(m-P_{2}+9)^{2}+6mP_{2}-204m+150P_{2}-P_{4}-1054\right)$
$\displaystyle-12mC_{1}^{2}C_{2}\bigl{(}Q_{tt}(m-P_{2}+9)-4m+P_{2}-34\bigr{)}-24m^{2}C_{1}C_{3}+12m^{2}C_{2}^{2}(Q_{tt}-2)\Bigr{]}\Bigr{\\}},$

so that we can write

$\hat{\varepsilon}=T_{8}\varepsilon^{8}+O(\varepsilon^{9}),$

that is,

$x_{k+1}-\alpha={\varPhi}(x_{k})-\alpha=O(\varepsilon_{k}^{8}),\ \
\varepsilon_{k}=x_{k}-\alpha.$

In this way we have proved the following assertion.

Theorem 2. If the initial approximation $x_{0}$ is sufficiently close to the
zero $\alpha$ of $f$ and the conditions (14) and (16) are valid, then the
order of the three-point family (12) is eight.

Note that $T_{8}$ tends to the asymptotic error constant of the family (12)
when $\varepsilon\to 0.$ According to Theorem 2, it follows that the three-
point family (12) is optimal.

The simplest forms of the weight functions $P$ and $Q$ are their truncated
Taylor series

$P(u)=1+2u+\beta\,u^{2}+(4-2\beta)\,u^{3},\quad
Q(u,v)=1+2u+v+4uv+(\beta+1)u^{2},\quad\beta\ \mbox{\rm is arbitrary}.$

Finding a list of particular weight functions $P$ and $Q$ is a routine work
and it is left to the interested reader as exercise.

4.3 Eighth-order family (7) derived by Schröder’s approach

The presented Schröder’s $f^{1/m}$-approach can also be applied for the
derivation of the method (7). Let us start from the following three-point
method for simple zeros

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-u_{k}H(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\\\\[9.0pt]
x_{k+1}=z_{k}-u_{k}P(u_{k})G(v_{k})L(w_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(17)

where

$u_{k}=\frac{f(y_{k})}{f(x_{k})},\quad v_{k}=\frac{f(z_{k})}{f(y_{k})},\quad
w_{k}=\frac{f(z_{k})}{f(x_{k})}\,.$

This method has the order 8 under the conditions

$\displaystyle H(0)=1,\ H^{\prime}(0)=2,\ P^{\prime}(0)=2P(0),\
P^{\prime\prime}(0)=P(0)(2+H^{\prime\prime}(0)),$ $\displaystyle
L^{\prime}(0)=2L(0),\
P^{\prime\prime\prime}(0)=P(0)(H^{\prime\prime\prime}(0)+6H^{\prime\prime}(0)-24),$
(18) $\displaystyle G(0)=0,\ G^{\prime}(0)=\dfrac{1}{L(0)P(0)},\
G^{\prime\prime}(0)=\dfrac{2}{L(0)P(0)}.$

The denotation of weight functions in (17) and (18) is adjusted to the
denotation used in [7].

Applying Schröder’s $f^{1/m}$-approach and (9) we obtain the family of three-
point methods for finding multiple zeros

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-mu_{k}H(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-mu_{k}P(u_{k})G^{*}(v_{k})L(w_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}\end{array}\right.\quad(k=0,1,\ldots),$
(19)

which is equivalent to (7) setting $mG^{*}(v)=G(v)$ (compare (7) and (19)).

Remark 4. Among other methods, the first author of this paper presented the
family of methods (17) in his lecture under the title Multipoint methods for
solving nonlinear equations at International conference Computational Methods
in Applied Mathematics (Berlin, 2012). However, due to the similarity to the
iterative formula (4.91) in the book [12, p 151], the first author did not
publish (17). Other possible sources of (17) are not known to the authors.

4.4 A modification of the family (7)

Since

$w=\Bigl{(}\frac{f(z)}{f(x)}\Bigr{)}^{1/m}=\Bigl{(}\frac{f(y)}{f(x)}\Bigr{)}^{1/m}\Bigl{(}\frac{f(z)}{f(y)}\Bigr{)}^{1/m}=u\cdot
v,$

the weight function $L(w)$ can be omitted in (7). Then we can construct in an
easy way the following family of three-point methods for finding multiple
zeros involving two parametric functions

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
z_{k}=y_{k}-mu_{k}H(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=z_{k}-mu_{k}v_{k}(1+2u_{k}v_{k})P(u_{k})G(v_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})}.\end{array}\right.$
(20)

Making suitable changes in the above program in Mathematica (see § 4.2), we
prove that the method (20) is of order eight under the following conditions:

$\displaystyle H(0)=1,\ \ H^{\prime}(0)=2,\ \ H^{\prime\prime}(0)=m+9,$
$\displaystyle P(0)=1,\ \ P^{\prime}(0)=2,\ \ P^{\prime\prime}(0)=m+11,\ \
P^{\prime\prime\prime}(0)=30+6m+H^{\prime\prime\prime}(0),$ $\displaystyle
G(0)=1,\ \ H^{\prime}(0)=1.$

4.5 Fourth order methods for multiple zeros

Chun constructed in [13] the following two-point family of iterative methods
for finding simple zeros

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=y_{k}-\dfrac{f(y_{k})}{h(u_{k})f^{\prime}(x_{k})},\quad
u_{k}=\dfrac{f(y_{k})}{f(x_{k})},\end{array}\right.$ (21)

where $h(u)$ is the weight function. The method (21) is of fourth order under
the condition $h(0)=1,\ h^{\prime}(0)=-1,\ |h^{\prime\prime}(0)|<+\infty.$ To
generate suitable two-point methods, sometimes it is necessary to develop the
function $h$ into Taylor’s or geometric series.

A slightly more direct approach without altering the weight function, based on
Chun’s idea, was given in [14] in the form

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=y_{k}-p(u_{k})\dfrac{f(y_{k})}{f^{\prime}(x_{k})},\quad
u_{k}=\dfrac{f(y_{k})}{f(x_{k})},\end{array}\right.$ (22)

The family (22) possesses the optimal order four if $p(0)=1,\ p^{\prime}(0)=2$
and $|p^{\prime\prime}(0)|<+\infty.$

Proceeding in the same way as in the case of the family of three-point methods
(10) and using (9), we obtain from (22) the fourth order two-point family of
iterative methods for finding multiple zeros

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-m\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=y_{k}-mu_{k}P(u_{k})\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\quad
u_{k}=\Bigl{(}\dfrac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m},\end{array}\right.$
(23)

where $P(u)$ is the weight function which satisfies

$P(0)=1,\ \ P^{\prime}(0)=2,$ (24)

that is, its truncated Taylor series is

$P(u)=1+2u.$

This result is expected since the iterative formula (23) coincides with the
first two steps of the family (12). Explicit expression of
$\hat{\varepsilon}=\hat{x}-\alpha$ (regarding (23)) is given by (15) (for
$\varepsilon_{z}$).

Several examples of functions that satisfy (24) are listed below.

$\displaystyle P(u)=\frac{1+\beta u}{1+(\beta-2)u}\ \ \mbox{\rm(of King's
type, see \cite[cite]{[\@@bibref{}{king}{}{}]})},\quad
P(u)=\Bigl{(}1+\frac{2u}{r}\Bigr{)}^{r},$ $\displaystyle P(u)=\frac{1+\gamma
u^{2}}{1-2u},\quad P(u)=\frac{1}{1-2u+au^{2}},\quad
P(u)=\frac{u^{2}+(c-2)u-1}{cu-1}.$

where $r\in\mbox{\boldmath{$Q$}}$ and $\gamma,\ a,\ c\in\mbox{\boldmath{$R$}}$
are arbitrary parameters.

4.6 Equivalence of methods

Some words about the equivalence of the methods presented in Section 2 and 3,
and the methods derived by Schröder’s approach in Section 4.

Equivalence 1. In Section 2 it was shown that the method (2) presented in [3]
(assuming slight simplification by taking $\lambda=0$ in (2)) is a special
case of the method (1) proposed by Zhou et al. [2]. Comparing iterative
methods (1) and (23) it is evident that both formulas are equivalent, which is
easily obtained taking $uP(u)\equiv G(u).$ Moreover, the family (23) requires
only two conditions $P(0)=1,\ P^{\prime}(0)=2$ compared to three conditions
$G(0)=0,\ G^{\prime}(0)=1,\ G^{\prime\prime}(0)=4.$ It is important to note
that the derivation of the family (1) and latter convergence analysis are more
complicated than for the family (23). Furthermore, the derivation of (23),
based on Schröder’s approach, is entirely natural and crystal clear.

Equivalence 2. The family (8), obtained by natural choice $a_{1}=1,\;a_{2}=0$
in the family (4) (proposed in [5]), is equivalent to the family (12), which
is evident setting $S(u)\equiv mP(u)$ and $R(u,v)\equiv mQ(u,v).$ The family
(4) (and, consequently, (8)) was presented without derivation and motivation,
while the family (12) was derived on an easy and obvious way using Schröder’s
approach.

Remark 5. Apart from the family (3), other optimal multipoint methods for
finding a simple zero (some of them are presented in the book [12]) can be
transformed by introducing $F(x)=f(x)^{1/m}$ to multipoint methods for
approximating a multiple zero keeping optimal order of convergence. This
subject is left to readers. However, the author of this paper does not expect
new papers in this direction since such methods are pretty expensive. This is
the subject of the next section.

##  5 Root-ratio methods are not competitive

In Section 4 we have demonstration a general scheme for constructing
multipoint methods for multiple zeros using basic iterative formulas for
simple zeros and Schröder’s $f^{1/m}$-idea. We observe that all previously
presented methods deal with real or complex values of the forms

$\Bigl{(}\frac{f(y_{k})}{f(x_{k})}\Bigr{)}^{1/m},\quad\Bigl{(}\frac{f(z_{k})}{f(x_{k})}\Bigr{)}^{1/m},\quad\Bigl{(}\frac{f(z_{k})}{f(y_{k})}\Bigr{)}^{1/m}.$

Computer algebra systems and computation software of digital computers often
meet the problem of finding the $m$-th root for arbitrary $m.$ In the case of
specific values of $m$ they find the principal value of the $m$-root $z^{1/m}$
($m$ is natural number given as numerical entry) among $m$ values of the
sought root in the form

$z^{1/m}=|z|^{1/m}\Bigl{(}\cos\tfrac{\theta}{m}+i\,\sin\tfrac{\theta}{m}\Bigr{)},\
\ \theta={\rm Arg}\,z\in(-\pi,\pi).$ (25)

From (25) it is clear that the computation of the $m$-root consumes a lot of
CPU time. The following test, implemented on PC with Intel-i7 processor and
clock speed 2.8 GHz, has given CPU times (expressed in $\mu sec$) for
different values of $m$ in calculation of
$\bigl{(}(a+ib)/(c+id)\bigr{)}^{1/m}.$ For the authenticity of the test, one
million experiments have been performed taking random numbers for $a,b,c,d$ in
each cycle to eliminate the use of possibly memorized data from previous
cycles. In the real case we set $b=0,\ d=0.$ The average CPU times for one
evaluation are given in Table 1.

$m$ | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
CPU (in $\mu sec$), real case | 5.25 | 22.7 | 32.24 | 31.86 | 32.25 | 32.03 | 33.1
CPU (in $\mu sec$), complex case | 13.77 | 54.82 | 64.42 | 64.07 | 65.83 | 68.05 | 66.3
Table 1: CPU times in the calculation of the $m$-th root of real and complex
numbers

In our experiments we observed that the CPU times for $m\geq 3$ almost do not
change in the real as well as the complex case (see Table 1). In this way we
are able to find reliable ratio of computation times for $m\geq 3$ (multiple
zeros) and $m=1$ (simple zeros):

$\text{\rm Real case:}\quad\frac{{\rm CPU}_{(m\geq 3)}}{{\rm
CPU}_{(m=1)}}\approx 6.2,\qquad\text{\rm Complex case:}\quad\frac{{\rm
CPU}_{(m\geq 3)}}{{\rm CPU}_{(m=1)}}\approx 4.8.$ (26)

Remark 6. For comparison purpose, we performed one million calculations of the
value of a polynomial of degree 20 with real coefficients and complex
argument, both chosen randomly in each cycles (Horner’s scheme was used).
Average CPU times for one evaluation was 116 $\mu sec$, which is only two
times slower than CPU time in calculation of $m$-th root in the case of
complex numbers (see Table 1). Therefore, root-ratio methods are expensive
form a computational point of view.

From (26) and Remark 6 we can draw trustworthy conclusion that root-ratio
multipoint methods, such as (1), (2), (3), (4), (5), (7), (8), (19), (20) and
other non-listed methods (if there exist), are pretty inefficient. As
mentioned in Remark 6, further work on the construction of root-ratio methods
is pointless and does not make any advance in the topic. Combining various
weight functions in order to derive “new methods” is rather a kind of play and
inevitably leads to minor modifications without a proper importance. This fact
was emphasized in [16] but the construction of modest modifications of
original contributions has continued. Different iterative formula does not
mean automatically that an advance in the topic was achieved.

Is there a good alternative method which is efficient and convenient for
applications? The answer is yes, and it is very likely known to many authors
who work in the area of iterative processes. Li et al. proposed in [17] the
following (optimal) two-point method of order four:

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\frac{2m}{m+2}\cdot\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=x_{k}-\dfrac{\frac{1}{2}m(m-2)A_{m}t_{k}-\frac{m^{2}}{2}}{1-R_{m}t_{k}}\cdot\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\quad
t_{k}=\dfrac{f^{\prime}(y_{k})}{f^{\prime}(x_{k})},\
A_{m}=\Bigl{(}\frac{m+2}{m}\Bigr{)}^{m}.\end{array}\right.$ (27)

Zhou et al. [18] later proposed the generalization of the method (27) in the
form

$\left\\{\begin{array}[]{l}y_{k}=x_{k}-\frac{2m}{m+2}\cdot\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\\\\[9.0pt]
x_{k+1}=x_{k}-\phi(t_{k})\cdot\dfrac{f(x_{k})}{f^{\prime}(x_{k})},\quad
t_{k}=\dfrac{f^{\prime}(y_{k})}{f^{\prime}(x_{k})},\end{array}\right.$ (28)

which has the order four under the specific conditions for the weight function
$\phi.$ See, also, [19]. The notion of “convenient” is explained in Remark 8.

We conclude this paper with two remarks of general interest.

Remark 7. The presence of arbitrary parameters in any zero-finding iterative
formula makes sense only if these parameters improve characteristics of
presented methods (such as acceleration of convergence, wider domain of
convergence, more accurate approximations, etc.). Otherwise, from an
algorithmic point of view, free parameters should be chosen so that an
iterative formula is as simple as possible – numerical analysts and
programmers will always choose the simpler formula in such a way that the best
characteristics of the employed methods are maintained. Inserting numerous
useless parameters does not make a method better or more general in the
genuine sense. Unfortunately, many authors construct “novel” iterative
formulas by adding parameters in an artificial way or by varying different
weight functions. In essence, such methods are only modest modifications of
existing methods and offer a little contribution to the topic.

Remark 8. It should be emphasized that very high accuracy of solutions of
nonlinear equations, provided by root-solvers of order eight or more, is not
needed for solving a huge number of practical problems; fourth order methods
(such as (27) and (28)) produce quite satisfactory results in practice. The
question “how many decimals of zero approximations do we really need in
practice?” is equivalent to the question “how many decimals of $\pi$ do we
really need in practice”?” A pretty convincing answer can be found in the
issue of NASA/JPL Edu, March 16, 2016:

$\bullet$ For interplanetary navigation with spacecraft Voyager 1 (launched in
1977, distant from Earth about 22 billion km), Jet Propulsion Labaratory
(California Institute of California, Pasadena, USA) and NASA use very accurate
calculations involving $\pi$ with most 15 decimal digits! Not more! The
distance error is about 5 cm!

$\bullet$ The radius of the visible universe is about 46 billion light years.
To express the circumference of a circle with this radius via the diameter a
hydrogen atom (the simplest atom) we need at most 40 decimal digits of $\pi!$

Acknowledgement. This work was supported by the Serbian Ministry of Education
and Science under Grant 174022.

## References

  * [1]
  * [2] X. Zhou, X. Chen, Y. Song, Families of third and fourth order methods for multiple roots of nonlinear equations, Appl. Math. Comput, 219 (2013), 6030–6038.
  * [3] M. Y. Lee, Y. I. Kim, Á. A. Magrenán, On the dynamics of tri-parametric family of optimal fourth-order multiple-zero finders with a weight function of the principal mth root of a function-to-function ratio, Appl. Math. Comput. 315 (2017), 564–590.
  * [4] B. Liu, X. Zhou, A new family of fourth-order methods for multiple roots of nonlinear equations, Nonlinear Anal. Model. Control, 2013 (18), 143–152.
  * [5] R. Behl, A. Cordero, S. Motsa, J. R. Torregrosa, An eighth-order family of optimal multiple root finders and its dynamics. Numer. Algor. 77 (2018), 1249–1272.
  * [6] F. Zafar, A. Cordero, R. Quratulain, J. R. Torregrosa, Optimal iterative methods for finding multiple roots of nonlinear equations using free parameters. J. Math. Chem. 56 (2018), 1884–1901.
  * [7] F. Zafar, A. Cordero, S. Sultana, J. R. Torregrosa, Optimal iterative methods for finding multiple roots of nonlinear equations using weight functions and dynamics, J. Comput. Appl. Math. 345 (2018), 352–374.
  * [8] H. T. Kung, J. F. Traub, Optimal order of one-point and multipoint iteration, Journal of the ACM 21 (1974), 643–651.
  * [9] E. Schröder, Über unendlich viele Algorithmen zur Auflösung der Gleichungen, Math. Ann. 2 (1870), 317–365.
  * [10] G.W. Stewart, On infinitely many algorithms for solving equations, drum.lib.umd.edu/handle/1903/577 (CS-TR-2990.ps) (University of Maryland Library).
  * [11] J. Džunić, M. S. Petković, L. D. Petković, A family of optimal three-point methods for solving nonlinear equations using two parametric functions, Appl. Math. Comput. 217 (2011), 7612–7619.
  * [12] M. S. Petković, B. Neta, L. D. Petković, J. Džunić, Multipoint Methods for Solving Nonlinear Equations, Elsevier, Amsterdam-Boston-Heidelberg-London-New York, 2013.
  * [13] C. Chun, Some fourth-order iterative methods for solving nonlinear equations, Appl. Math. Comput. 195 (2008), 454–459.
  * [14] M. S. Petković, L. D. Petković, Families of optimal multipoint methods for solving nonlinear equations: a survey, Appl. Anal. Discrete Math. 4 (2010), 1–22.
  * [15] R. F. King, A family of fourth order methods for nonlinear equations, SIAM J. Numer. Anal. 10 (1973), 876–879.
  * [16] M. S. Petković, B. Neta, L. D. Petković, J. Džunić, Multipoint methods for solving nonlinear equations: a survey, Appl. Math. Comput. 226 (2014) 635–660.
  * [17] S. Li, X. Liao, L. Cheng, A new-fourth-order iterative method for finding multiple root of nonlinear equations, Appl. Math. Comput. 215 (2009), 1288–1292.
  * [18] X. Zhou, X. Chen, Y. Song, Construction of higher order methods for multiple roots of nonlinear equations, J. Comput. Appl. Math. 235 (2011), 4199–4206.
  * [19] S. Li, L. Cheng, B. Neta, Some fourth-order nonlinear solvers with closed formulae for multiple roots, Comp. Math. Appl. 59 (2010), 126–135 .
  * [20] J. Džunić, M. S. Petković, L. D. Petković, A family of optimal three-point methods for solving nonlinear equations using two parametric functions, Appl. Math. Comput. 217 (2011) 7612–7619.
  * [21]