A Two-Parameter Family of Fourth-Order Iterative Methods with Optimal Convergence for Multiple Zeros

Abstract

We develop a family of fourth-order iterative methods using the weighted harmonic mean of two derivative functions to compute approximate multiple roots of nonlinear equations. They are proved to be optimally convergent in the sense of Kung-Traub’s optimal order. Numerical experiments for various test equations confirm well the validity of convergence and asymptotic error constants for the developed methods.

1. Introduction

A development of new iterative methods locating multiple roots for a given nonlinear equation deserves special attention on both theoretical and numerical interest, although prior knowledge about the multiplicity of the sought zero is required [1]. Traub [2] discussed the theoretical importance of multiple-root finders, although the multiplicity is not known a priori by stating: “since the multiplicity of a zero is often not known a priori, the results are of limited value as far as practical problems are concerned. The study is, however, of considerable theoretical interest and leads to some surprising results.” This motivates our analysis for multiple-root finders to be shown in this paper. In case the multiplicity is not known, interested readers should refer to the methods suggested by Wu and Fu [3] and Yun [4, 5].

Various iterative schemes finding multiple roots of a nonlinear equation with the known multiplicity have been proposed and investigated by many researchers [6–12]. Neta and Johnson [13] presented a fourth-order method extending Jarratt's method. Neta [14] also developed a fourth-order method requiring one-function and three-derivative evaluations per iteration grounded on a Murakami’s method [15]. Shengguo et al. [16] proposed the following fourth-order method which needs evaluations of one function and two derivatives per iteration for 𝑥₀ chosen in a neighborhood of the sought zero 𝛼 of 𝑓(𝑥) with known multiplicity 𝑚≥1 as follows: 𝑥_𝑛+1=𝑥_𝑛−𝛽𝑓^′𝑥_𝑛+𝜙𝑓^′𝑦_𝑛𝑓^′𝑥_𝑛+𝛿𝑓^′𝑦_𝑛𝑓𝑥_𝑛𝑓^′𝑥_𝑛,𝑛=0,1,2,…, (1) where 𝑦_𝑛=𝑥_𝑛−(2𝑚/(𝑚+2))(𝑓(𝑥_𝑛)/𝑓^′(𝑥_𝑛)),  𝛽=−(𝑚²/2), 𝜙=(1/2)(𝑚(𝑚−2)/(𝑚/(𝑚+2))^𝑚), and 𝛿=−(𝑚/(𝑚+2))^−𝑚 with the following error equation: 𝑒_𝑛+1=𝐾₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛,𝑛=0,1,2,…, (2) where 𝐾₄=((−2+2𝑚+2𝑚²+𝑚³)/3𝑚⁴(𝑚+1)³)𝜃³₁ − (1/𝑚(𝑚+1)²(𝑚+2))𝜃₁𝜃₂ + (𝑚/(𝑚+2)³(𝑚+1)(𝑚+3))𝜃₃, 𝑒_𝑛=𝑥_𝑛−𝛼, and 𝜃_𝑗=𝑓^(𝑚+𝑗)(𝛼)/𝑓^(𝑚)(𝛼) for 𝑗=1,2,3.

Based on Jarratt [17] scheme for simple roots, J. R. Sharma and R. Sharma [18] developed the following fourth order of convergent scheme: 𝑥_𝑛+1=𝑥_𝑛𝑓𝑥−𝐴^𝑛𝑓^′𝑥_𝑛𝑓𝑥−𝐵^𝑛𝑓^′𝑦_𝑛−𝐶𝑓(𝑥_𝑛)𝑓^′(𝑦_𝑛)²𝑓(𝑥_𝑛)𝑓^′(𝑥_𝑛)⁻¹,𝑛=0,1,2,…, (3) where 𝐴=(1/8)𝑚(𝑚³−4𝑚+8), 𝐵=−(1/4)𝑚(𝑚−1)(𝑚+2)²(𝑚/(𝑚+2))^𝑚, 𝐶=(1/8)𝑚(𝑚+2)³(𝑚/(𝑚+2))^2𝑚, and 𝑦_𝑛=𝑥_𝑛−(2𝑚/(𝑚+2))(𝑓(𝑥_𝑛)/𝑓^′(𝑥_𝑛)) and derived the error equation below: 𝑒_𝑛+1=−8+12𝑚+14𝑚²+14𝑚³+6𝑚⁴+𝑚⁵3𝑚⁴(𝑚+2)²𝒜³₁−1𝑚𝒜₁𝒜₂+𝑚(𝑚+2)²𝒜₃𝑒⁴_𝑛+𝑂𝑒⁵_𝑛, (4) where 𝒜_𝑗=(𝑚!/(𝑚+𝑗)!)(𝑓^(𝑚+𝑗)(𝛼)/𝑓^(𝑚)(𝛼))=(𝑚!/(𝑚+𝑗)!)𝜃_𝑗 for 𝑗=1,2,3.

The above error equation can be expressed in terms of 𝜃_𝑗(𝑗=1,2,3) as follows: 𝑒_𝑛+1=𝐶₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛, (5) where 𝐶₄=((−8+12𝑚+14𝑚²+14𝑚³+6𝑚⁴+𝑚⁵)/3𝑚⁴(𝑚+1)³(𝑚+2)²)𝜃³₁ − (1/𝑚(𝑚+1)²(𝑚+2))𝜃₁𝜃₂ + (𝑚/(𝑚+2)³(𝑚+1)(𝑚+3))𝜃₃.

We now proceed to develop a new iterative method finding an approximate root 𝛼 of a nonlinear equation 𝑓(𝑥)=0, assuming the multiplicity of 𝛼 is known. To do so, we first suppose that a function 𝑓∶ℂ→ℂ has a multiple root 𝛼 with integer multiplicity 𝑚≥1 and is analytic in a small neighborhood of 𝛼. Then we propose a new iterative method free of second derivatives below with an initial guess 𝑥₀ sufficiently close to 𝛼 as follows: 𝑥_𝑛+1=𝑦_𝑛𝑓𝑥−𝑎^𝑛𝑓^′𝑥_𝑛𝐹𝑦−𝑏^𝑛𝑓^′𝑥_𝑛𝑓𝑥−𝑐^𝑛𝑓^′𝑦_𝑛,𝑛=0,1,2,…, (6) where 𝑦_𝑛=𝑥_𝑛𝑓𝑥−𝛾^𝑛𝑓^′𝑥_𝑛,𝐹𝑦_𝑛=𝑓𝑥_𝑛+𝑦_𝑛−𝑥_𝑛𝜆𝑓^′𝑥_𝑛𝑓^′𝑦_𝑛𝑓^′𝑥_𝑛+𝜌𝑓^′𝑦_𝑛, (7) with 𝑎,𝑏,𝑐,𝛾,𝜆, and 𝜌 as parameters to be chosen for maximal order of convergence [2, 19]. One should note that 𝐹(𝑦_𝑛) is obtained from Taylor expansion of 𝑓(𝑦_𝑛) about 𝑥_𝑛 up to the first-order terms with weighted harmonic mean [20] of 𝑓^′(𝑥_𝑛) and 𝑓^′(𝑦_𝑛).

Theorem 1 shows that proposed method (6) possesses 2 free parameters 𝜆 and 𝜌. A variety of free parameters 𝜆 and 𝜌 give us an advantage that iterative scheme (6) can develop various numerical methods. One can often have a freedom to select best suited parameters 𝜆 and 𝜌 for a sought zero 𝛼. Several interesting choices of 𝜆 and 𝜌 further motivate our current analysis. As seen in Table 1, we consider five kinds of methods Y1, Y2, Y3, Y4, and Y5 list selected parameters (𝜆,𝜌), and the corresponding values (𝑎,𝑏,𝑐), respectively.

Table 1 

Various methods with typical choice of parameters (𝜆,𝜌,𝑎,𝑏,and𝑐).

If 𝜆=−(𝑚((𝑚+2)/𝑚)^𝑚+2/(𝑚²+2𝑚+4)) and 𝜌=−((𝑚+2)/𝑚)^𝑚 are selected, then we obtain 𝑎=0, 𝑏=−(𝑚(𝑚²+2𝑚+4)/2(𝑚+2)), and 𝑐=0, in which case iterative scheme (6) becomes method Y5 mentioned above and blackuces to iterative scheme (1) developed by Shengguo et al. [16].

In this paper, we investigate the optimal convergence of the fourth-order methods for multiple-root finders with known multiplicity in the sense of optimal order claimed by Kung-Traub [21] and derive the error equation. We find that our proposed schemes require one evaluation of the function and two evaluations of first derivative and satisfy the optimal order. In addition, through a variety of numerical experiments we wish to confirm that the proposed methods show well the convergence behavior predicted by the developed theory.

2. Convergence Analysis

In this section, we describe a choice of parameters 𝑎,𝑏, and 𝑐 in terms of 𝜆 and 𝜌 to get fourth-order convergence for our proposed scheme (6).

Theorem 1. Let 𝑓∶ℂ→ℂ have a zero 𝛼 with integer multiplicity 𝑚≥1 and be analytic in a small neighborhood of 𝛼. Let 𝜅=(𝑚/(𝑚+2))^𝑚, 𝜃_𝑗=(𝑓^(𝑚+𝑗)(𝛼)/𝑓^(𝑚)(𝛼))for𝑗∈ℕ. Let 𝑥₀ be an initial guess chosen in a sufficiently small neighborhood of 𝛼. Let 𝑎=(𝑚/(𝑚+2)){𝑚−(1/8)𝑚(𝑚+2)³(1+𝜅𝜌)−(𝑚+2)²(𝑚(−1+2𝜅𝜆)−(𝑚+2)𝜅𝜌)(𝑚+(𝑚+2)𝜅𝜌)²/16𝑚𝜅𝜆}, 𝑏=−((𝑚+2)(𝑚+(𝑚+2)𝜅𝜌)³/16𝜅𝜆), 𝑐=(1/8)𝑚(𝑚+2)³𝜅(1+𝜅𝜌), and 𝛾=2𝑚/(𝑚+2). Let 𝜆,𝜌∈ℝ be two free constant parameters. Then iterative method (6) is of order four and defines a two-parameter family of iterative methods with the following error equation: 𝑒_𝑛+1=𝜓₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛,𝑛=0,1,2,…, (8) where 𝑒_𝑛=𝑥_𝑛−𝛼 and 𝜓₄=8+2𝑚+6𝑚²+4𝑚³+𝑚⁴−24𝜅𝜌/𝑚+(𝑚+2)𝜅𝜌3𝑚⁴(𝑚+1)³()𝜃𝑚+2³¹−1𝑚(𝑚+1)²()𝜃𝑚+2¹𝜃₂+𝑚(𝑚+2)³()𝜃𝑚+1)(𝑚+3³. (9)

Proof. Using Taylor's series expansion about 𝛼, we thave the following relations: 𝑓𝑥_𝑛=𝑓^(𝑚)(𝛼)𝑒𝑚!^𝑛𝑚×1+𝐴₁𝑒_𝑛+𝐴₂𝑒²_𝑛+𝐴₃𝑒³_𝑛+𝐴₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛,𝑓^′𝑥_𝑛=𝑓^(𝑚)(𝛼)()!𝑒𝑚−1^𝑛𝑚−1×1+𝐵₁𝑒_𝑛+𝐵₂𝑒²_𝑛+𝐵₃𝑒³_𝑛+𝐵₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛, (10)(11) where 𝐴_𝑘=(𝑚!/(𝑚+𝑘)!)𝜃_𝑘, 𝐵_𝑘=((𝑚−1)!/(𝑚+𝑘−1)!)𝜃_𝑘, and 𝜃_𝑘=(𝑓^(𝑚+𝑘)(𝛼)/𝑓^(𝑚)(𝛼)) for 𝑘∈ℕ.
Dividing (10) by (11), we obtain 𝑓𝑥_𝑛𝑓^′𝑥_𝑛=1𝑚𝑒_𝑛−𝐾₁𝑒²_𝑛−𝐾₂𝑒³_𝑛+𝐾₃𝑒⁴_𝑛+𝑂𝑒⁵_𝑛, (12) where 𝐾₁=−𝐴₁+𝐵₁, 𝐾₂=−𝐴₂+𝐴₁𝐵₁−𝐵²₁+𝐵₂, and 𝐾₃=−𝐴₃+𝐴₂𝐵₁−𝐴₁𝐵²₁+𝐵³₁+𝐴₁𝐵₂−2𝐵₁𝐵₂+𝐵₃.
Expressing 𝛾=𝑚(1−𝑡) in terms of a new parameter 𝑡 for algebraic simplicity, we get 𝑦_𝑛=𝑥_𝑛𝑓𝑥−𝛾^𝑛𝑓^′𝑥_𝑛=𝛼+𝑡𝑒_𝑛+𝐾₁(1−𝑡)𝑒²_𝑛+𝐾₂(1−𝑡)𝑒³_𝑛+𝐾₃(1−𝑡)𝑒⁴_𝑛+𝑂𝑒⁵_𝑛. (13) Since 𝑓^′(𝑦_𝑛) can be expressed from 𝑓^′(𝑥_𝑛) in (11) with 𝑒_𝑛 substituted by (𝑦_𝑛−𝛼) from (13), we get 𝑓^′𝑦_𝑛=𝑓^(𝑚)(𝛼)()!𝑒𝑚−1^𝑛𝑚−1×𝑡^𝑚−1+𝑡^𝑚−2𝐵₁𝑡²−𝐾₁(𝑚−1)(𝑡−1)𝑒_𝑛+𝑡^𝑚−3×𝐾²₁(𝑚−2)(𝑚−1)(𝑡−1)²−2𝑚𝐵₁𝐾₁(𝑡−1)𝑡²+2𝑡𝐵₂𝑡³−𝐾₂(𝑚−1)(𝑡−1)2×𝑒²_𝑛+𝑡^𝑚−46−𝐾³₁(𝑚−3)(𝑚−2)(𝑚−1)(𝑡−1)³+3𝐵₁𝐾²₁(𝑚−1)𝑚(𝑡−1)²𝑡²−6𝐾₁𝑡()×𝑡−1−𝐾₂(𝑚−2)(𝑚−1)(𝑡−1)+𝐵₂(𝑚+1)𝑡³−6𝑡²𝐵₁𝐾₂𝑚𝑡(𝑡−1)−𝐵₃𝑡⁴+𝐾₃(𝑚−1)(𝑡−1)𝑒³_𝑛+𝑂𝑒⁴_𝑛. (14)
With the aid of symbolic computation of Mathematica [22], we substitute (10)–(14) into proposed method (6) to obtain the error equation as 𝑒_𝑛+1=𝑦_𝑛𝑓𝑥−𝛼−𝑎ⁿ𝑓^′𝑥_𝑛𝐹𝑦−𝑏^𝑛𝑓^′𝑥_𝑛𝑓𝑥−𝑐^𝑛𝑓^′𝑦_𝑛=𝜓₁𝑒_𝑛+𝜓₂𝑒²_𝑛+𝜓₃𝑒³_𝑛+𝜓₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛, (15) where 𝜓₁=𝑡−((𝑎+𝑏+𝑐𝑡^1−𝑚)/𝑚)−(𝑏(𝑡−1)𝑡^𝑚−1𝜆/(1+𝑡^𝑚−1𝜌)), and the coefficient 𝜓_𝑖(𝑖=2,3,4) may depend on parameters 𝑡,𝑎,𝑏,𝑐,𝜆, and 𝜌.
Solving 𝜓₁=0 and 𝜓₂=0 for 𝑎 and 𝑏, respectively, we get after simplifications 𝑎=−𝑏+𝑡𝑚−𝑐𝑡^−𝑚−𝑏𝑚(𝑡−1)𝑡^𝑚−1𝜆1+𝑡^𝑚−1𝜌,𝑡𝑏=^2−2𝑚𝑚𝑡^𝑚+𝑐(𝑡−1)(1+𝑚(𝑡−1)+𝑡)1+𝑡^𝑚−1𝜌²𝑚(𝑡−1)²(())𝜆.1+𝑚𝑡−1+𝑡 (16) Putting 𝜓₃=𝜓₃₁𝜃₁²+𝜓₃₂𝜃₂, we have 𝜓₃₁=𝑡^−1−𝑚−2𝑐(𝑡−1)²𝑡(1+𝑚(𝑡−1)+𝑡)³+𝑚𝑡^𝑚𝑃₁+𝑡^𝑚𝑃₂𝜌2𝑚³(𝑚+1)²(1+𝑚(𝑡−1)+𝑡)𝑡+𝑡^𝑚𝜌, (17)𝜓₃₂=𝑡(𝑚−(𝑚+2)𝑡)𝑚(()),𝑚+1)(𝑚+2)(1+𝑚𝑡−1+𝑡 (18) where 𝑃₁=𝑡[2−3𝑚+𝑚²+(2+𝑚−3𝑚²)𝑡+2(𝑚+1)²𝑡²] and 𝑃₂=2𝑡²(2+𝑡)+𝑚²(𝑡−1)[1+2(𝑡−1)𝑡]+𝑚(1−3𝑡+4𝑡³).
Observe that 𝜓₃₂=0 is satisfied with 𝑡=𝑚/(𝑚+2). Solving 𝜓₃₁=0 for 𝑐, we get 𝑡𝑐=^𝑚−1𝑚𝑃₁+𝑡^𝑚𝑃₂𝜌2(𝑡−1)²(1+𝑚(𝑡−1)+𝑡)³. (19) Substituting 𝑡=𝑚/(𝑚+2) into (16) and (19) with 𝜅=(𝑚/(𝑚+2))^𝑚, we can rearrange these expressions to obtain 𝑚𝑎=𝑚()𝑚+2𝑚−𝑚+2³1+𝜅𝜌8−(𝑚+2)²𝑚(−1+2𝜅𝜆)−(𝑚+2)𝜅𝜌𝑚+(𝑚+2)𝜅𝜌²,16𝑚𝜅𝜆 (20)()()𝑏=−𝑚+2𝑚+𝑚+2𝜅𝜌³,116𝜅𝜆𝑐=8𝑚(𝑚+2)³𝜅1+𝜅𝜌. (21)(22) Calculating by the aid of symbolic computation of Mathematica [22], we arrive at the error equation below: 𝑒_𝑛+1=𝜓₄𝑒⁴_𝑛+𝑂𝑒⁵_𝑛, (23) where 𝜓₄=((8+2𝑚+6𝑚²+4𝑚³+𝑚⁴−(24𝜅𝜌/(𝑚+(𝑚+2)𝜅𝜌))/3𝑚⁴(𝑚+1)³(𝑚+2))𝜃³₁−(1/𝑚(𝑚+1)²(𝑚+2))𝜃₁𝜃₂ + (𝑚/(𝑚+2)³(𝑚+1)(𝑚+3))𝜃₃ with 𝜅=(𝑚/(𝑚+2))^𝑚.

It is interesting to observe that error equation (23) has only one free parameter 𝜌, being independent of 𝜆. Table 1 shows typically chosen parameters 𝜆 and 𝜌 and defines various methods 𝑌_𝑘,(𝑘=1,2,…,5) derived from (6). Method Y5 results in the iterative scheme (1) that Shengguo et al. [16] suggested.

3. Numerical Examples and Conclusion

In this section, we have performed numerical experiments using Mathematica Version 5 program to convince that the optimal order of convergence is four and the computed asymptotic error constant |𝑒_𝑛+1/𝑒⁴_𝑛| agrees well with the theoretical value 𝜂. To achieve the specified sufficient accuracy and to handle small number divisions appearing in asymptotic error constants, we have assigned 300 as the minimum number of digits of precision by the command $𝑀𝑖𝑛𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛=300 and set the error bound 𝜖 to 10⁻²⁵⁰ for |𝑥_𝑛−𝛼|<𝜖. We have chosen the initial values 𝑥₀ close to the sought zero 𝛼 to get fourth-order convergence. Although computed values of 𝑥_𝑛 are truncated to be accurate up to 250 significant digits and the inexact value of 𝛼 is approximated to be accurate enough about up to 400 significant digits (with the command 𝐹𝑖𝑛𝑑𝑅𝑜𝑜𝑡[𝑓[𝑥],{𝑥,𝑥₀}], 𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛𝐺𝑜𝑎𝑙→400, and𝑊𝑜𝑟𝑘𝑖𝑛𝑔𝑃𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛→600), we list them up to 15 significant digits because of the limited space.

As a first example with a double zero 𝛼=√3 and an initial guess 𝑥₀=1.58, we select a test function 𝑓(𝑥)=cos(𝜋𝑥²/6)log(𝑥²−√3𝑥+1). As a second experiment, we take another test function 𝑓(𝑥)=(16+𝑥²)³(log(𝑥²+17))⁴ with a root 𝛼=−4𝑖 of multiplicity 𝑚=7 and with an initial value 𝑥₀=−3.94𝑖.

Taking another test function 𝑓(𝑥)=(1−sin(𝑥²))(log(2𝑥²−𝜋+1))⁴ with a root 𝛼=−𝜋/2 of multiplicity 𝑚=6, we select 𝑥₀=−1.18 as an initial value.

Throughout these examples, we confirm that the order of convergence is four and the computed asymptotic error constant |𝑒_𝑛+1/𝑒⁴_𝑛| approaches well the theoretical value 𝜂. The order of convergence and the asymptotic error constant are clearly shown in Tables 2, 3, and 4 reaching a good agreement with the theory developed in Section 2.

Table 2 

Convergence behavior with 𝑓(𝑥)=cos(𝜋𝑥²/6)log(𝑥²−√3𝑥+1),  (𝑚,𝜆,𝜌)=(2,5,5),  𝛼=√3.

Table 3 

Convergence behavior with 𝑓(𝑥)=(16+𝑥²)³(log(𝑥²+17))⁴,  (𝑚,𝜆,𝜌)=(7,1,−1),  𝛼=−4𝑖.

Table 4 

Convergence behavior with 𝑓(𝑥)=(1−sin(𝑥²))(log(2𝑥²−𝜋+1))⁴,  (𝑚,𝜆,𝜌)=(6,1,−1),  𝛼=−𝜋/2.

The additional test functions 𝑓₁,𝑓₂,…,𝑓₇ listed below further confirm the convergence behavior of our proposed method (6). 𝑓₁(𝑥)=𝑥⁷−𝑥²−7𝑥⁵𝑥+sin𝑥,𝛼=1.3657,𝑚=1,₀𝑓=1.32,₂(𝑥)=𝑒^{𝑥5−𝑥2−7}−17+𝑥²−𝑥⁵,𝛼=−1.16−0.95𝑖,𝑚=2,𝑥₀𝑓=−1.14−0.92𝑖,₃(𝑥)=𝑥⁶−8²log𝑥⁶−7,𝛼=√𝑥2,𝑚=3,₀𝑓=1.39,₄(𝑥)=3−𝑥+𝑥²⁴cot𝑥²√+1,𝛼=1−11𝑖2,𝑚=4,𝑥₀𝑓=0.47−1.79𝑖,₅(𝑥)=𝑥⁴−9𝑥³+4𝑥²−33𝑥−27log(𝑥−8)³,𝛼=9,𝑚=5,𝑥₀𝑓=8.79,₆(𝑥)=log(1−𝜋+𝑥)⁶𝑥,𝛼=𝜋,𝑚=6,₀𝑓=3.09,₇(𝑥)=𝑒^3𝑥+𝑥2−1cos³𝜋𝑥²𝑥18log³−𝑥²+37³,𝛼=−3,𝑚=7,𝑥₀=−2.88. (24)

Table 5 shows the convergence behavior of |𝑥_𝑛−𝛼| among methods S, J, Y1, Y2, Y3, and Y4, where S denotes the method proposed by Shengguo et al. [16], J the method proposed by J. R. Sharma and R. Sharma [18], the methods Y1 to Y4 are described in Table 1. It is certain that proposed method (6) needs one evaluation of the function 𝑓 and two evaluations of the first derivative 𝑓^′ per iteration. Consequently, the corresponding efficiency index [2] is found to be 4^1/3≈1.587, which is optimally consistent with the conjecture of Kung-Traub [21]. For the particularly chosen test functions in these numerical experiments, methods Y1 to Y4 have shown better accuracy than methods S and J.

Table 5 

Comparison of |𝑥_𝑛−𝛼| for high-order iterative methods.

Nevertheless, the favorable performance of proposed scheme (6) is not always expected since no iterative method always shows best accuracy for all the test functions. If we look at the asymptotic error equation 𝜂=𝜂(𝑓,𝛼,𝑝)=lim_𝑛→∞|(𝑥_𝑛+1−𝛼)/(𝑥_𝑛−𝛼)^𝑝| closely, we should note that the computational accuracy is sensitively dependent on the structures of the iterative methods, the sought zeros, convergence orders, the test functions, and good initial values.

It is important to properly choose enough number of precision digits. If 𝑒_𝑘 is small, 𝑒^𝑝_𝑘 gets much smaller, as 𝑘 increases. If the number of precision is small and error bound 𝜖 is not small enough, the term 𝑒^𝑝_𝑘 causes a great loss of significant digits due to magnified round-off errors. This hinders us from verifying 𝑝 and 𝜂 more accurately.

Bold-face numbers in Table 5 refer to the least error until the prescribed error bound is met. This paper has confirmed optimal fourth-order convergence proved the correct error equation for proposed iterative methods (6), using the weighted harmonic mean of two derivatives to find approximate multiple zeros of nonlinear equations. We remark that the error equation of (6) contains only one free parameter 𝜌, being independent of 𝜆.

We still further need to discuss some aspects of root finding for ill-conditioned problems as well as sensitive dependence of zeros on initial values for iterative methods. As is well known, a high-degree (say, degree higher than 20, taking the multiplicity of a zero into account) polynomial is very likely to be ill conditioned. In this case, small changes in the coefficients can greatly alter the zeros. The small changes can occur as a result of rounding process in computing coefficients. Minimal round-off errors may improve the root finding of ill-conditioned problems. Certainly multi-precision arithmetic should be used in conjunction with optimized algorithms reducing round-off errors. High-order methods with the asymptotic error constant of small magnitude are preferred for locating zeros with relatively good accuracy. Locating zeros for ill-conditioned problems is generally believed to be a difficult task.

It is also important to properly choose close initial values near the root for guaranteed convergence of the proposed method. Indeed, initial values are chaotic [23] to the convergence of the root 𝛼. The following statement quoted from [24] is not too much to emphasize the importance of selected initial values: “A point that belongs to the non-convergent region for a particular value of the parameter can be in the convergent region for another parameter value, even though the former might have a higher order of convergence than the second. This then indicates that showing whether a method is better than the other should not be done through solving a function from a randomly chosen initial point and comparing the number of iterations needed to converge to a root.”

Since our current analysis aims on the convergence of the proposed method, initial values [25–27] are selected in a small neighborhood of 𝛼 for guaranteed convergence. Thus the chaotic behavior of 𝑥₀ on the convergence should be separately treated under the different subject in future analysis. On the one hand, future research may be more strengthened with the graphical analysis on the convergence including chaotic fractal basins of attractions. On the other hand, rational approximations [1] provide rich resources of future research on developing new high-order optimal methods for multiple zeros.

Acknowledgments

The first author (Y. I. Kim) was supported by the Research Fund of Dankook University in 2011; the corresponding author (Y. H. Geum) was also supported by the National Research Foundation of Korea funded by the Ministry of Education, Science, and Technology (Project No. 2011-0014638). In addition, the authors would like to give special thanks to anonymous referees for their valuable suggestions and comments on this paper.