# Study of Dynamical Behavior and Stability of Iterative Methods for Nonlinear Equation with Applications in Engineering

Authors: Naila Rafiq; Saima Akram; Nazir Ahmad Mir; Mudassir Shams
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this article, we first construct a family of optimal 2-step iterative methods for finding a single root of the nonlinear equation using the procedure of weight function. We then extend these methods for determining all roots simultaneously. Convergence analysis is presented for both cases to show that the order of convergence is 4 in case of the single-root finding method and is 6 for simultaneous determination of all distinct as well as multiple roots of a nonlinear equation. The dynamical behavior is presented to analyze the stability of fixed and critical points of the rational operator of one-point iterative methods. The computational cost, basins of attraction, efficiency, log of the residual, and numerical test examples show that the newly constructed methods are more efficient as compared with the existing methods in the literature.


---

## Body

## 1. Introduction

To solve the nonlinear equation(1)fs=0,is the oldest problem of engineering in general and in mathematics in particular. These nonlinear equations have diverse applications in many areas of science and engineering. To find the roots of (1), we look towards iterative schemes, which can be classified as to approximate single root and all roots of (1). In this article, we are going to work on both types of iterative schemes. A lot of iterative methods of different convergence orders already exist in the literature (see [1–11]) to approximate the roots of (1). Ostrowski [7] defined efficiency index I of these iterative methods in terms of their convergence order k and the number of function evaluations per iteration, say u, i.e.,(2)I=k1/u.An iterative method is said to be optimal according to Kung–Traub conjecture [1] if(3)k=2u−1,holds. The aforementioned methods are used to approximate one root at a time. However, mathematicians are also interested in finding of all roots of (1) simultaneously. This is due to the fact that simultaneous iterative methods are very popular due to their wider region of convergence, are more stable as compared to single-root finding methods, and implemented for parallel computing as well. More detail on single as well as simultaneous determination of all roots can be found in [1, 12–24] and references cited therein.The most famous of the single-root finding method is the classical Newton–Raphson method:(4)si+1=si−fsif′si,i=1,2,….Method (4) requires one evaluation of the function and one of its first derivative to achieve optimal order 2 having efficiency 1.41 using the Traub conjecture. Using Weierstrass’ correction [15],(5)fsif′si=ws=fsi∏j≠ij=1ns−sj,in (4), we get the classical Weierstrass—Dochive method to approximate all roots of nonlinear equation (1) given as(6)si+1=si−fsi∏j≠ij=1nsi−sj.Method (6) has convergence order 2. Later, Aberth-Ehrlich [14] presented the 3rd-order simultaneous method given as(7)si+1=si−11/Nsi−∑j≠ij=1n1/si−sj,where Nsi=f′si/fsi.The main aim of this paper is to construct the family of optimal fourth-order single-root finding methods using the procedure of weight function and then convert them into simultaneous iterative methods for finding all distinct as well as multiple roots of nonlinear equation (1). Using the complex dynamical system, we will be able to choose those values of parameters of iterative methods which give a wider convergence area on initial approximations.


## 2. Construction of the Method and Convergence Analysis

King [25] presented the following optimal fourth-order method (abbreviated as MM1):(8)yi=si−fsif′si,zi=yi−fsi+βfyifsi+β−2fyifyif′si.Chun [26] gave the fourth-order optimal method as (abbreviated as MM2)(9)yi=si−fsif′si,zi=si−fsif′si1+u+2u2,whereu=fyifsi.Cordero et al. [3] proposed the fourth-order optimal method as (abbreviated as MM3)(10)yi=si−fsif′si,vi=si−fsi+fyif′si,zi=vi−fyi22fsi+fyifsi2f′si.Chun [27] introduced the fourth-order optimal method as (abbreviated as MM4)(11)yi=si−fsif′si,zi=si−fsif′si−fsif′sifsifsi−fyi.Here, we propose the following family of iterative methods:(12)yi=si−fsif′si,zi=yi−fsif′siΓu,whereu=fyifsi.For iterative scheme (12), we have the following convergence theorem:Theorem 1.
Letζ∈I be a simple root of a sufficiently differential functionf:I⊆ℜ⟶ℜ in an open interval I. If s0 is sufficiently close to ζ and Γ be a real-valued function satisfying Γ0=0,Γ′0=1,Γ″0=4,andΓ‴0<∞, then the convergence order of the family of iterative method (12) is 4 and satisfies the following error equation:(13)ei+1=5c23−c2c3−16Γ‴0c23ei4+Oei5,where cm=fmζ/m!f′ζ;m≥2.Proof.
Letζ be a simple root of f, and si=ζ+ei. By Taylor’s series expansion of fsi about s=ζ taking fζ=0, we get(14)fsi=f′ζei+c2ei2+c3ei3+c4ei4+Oei5,(15)f′si=f′ζ1+2c2ei+3c3ei2+Oei3.
Dividing (14) by (10), we have(16)fsif′si=ei+c2ei2+2c22−4c3ei3+Oei4.
This gives(17)yi=ζ+c2ei2+−2c22+2c3ei3+⋯u=fyifsi=c2ei+−3c22+2c3ei2+8c23−10c2c3+3c4ei3+⋯.
Thus, using Taylor series, we have(18)Γu=Γ0+uΓ′0+u22!Γ″0+⋯Γu=Γ0+c2Γ′0ei+12c22Γ″0−3Γ′0+2c3Γ′0ei2+⋯.
Now, takingΓq=0,Γ0=1,andΓ″0=4 in equation (12) and on simplification gives(19)ei+1=5c23−c2c3−16Γ‴0c23ei4+Oei5.
Hence, it proves fourth-order convergence.


### 2.1. The Concrete Fourth-Order Family of Methods

We now construct some concrete forms of the family of fourth-order methods from the family of methods (12) by choosing weight function Γu containing an arbitrary real number β as provided in Table 1, satisfying the condition Γ0=0,Γ′0=1,Γ″0=4andΓ‴0<∞ of Theorem 1 with β as a real number. Therefore, we get the following new five families of iterative methods:Method-1 (abbreviated as BB1):(20)yi=si−fsif′si,zi=yi−fsif′siu+2u2+βu4.Method-2 (abbreviated as BB2):(21)yi=si−fsif′si,zi=yi−fsif′siu1−2u+βu2.Method-3 (abbreviated as BB3):(22)yi=si−fsif′si,zi=yi−fsif′siu++2u2+βu3.Method-4 (abbreviated as BB4):(23)yi=si−fsif′si,zi=yi−fsif′siu1−βu2+2u2.Method-5 (abbreviated as BB5):(24)yi=si−fsif′si,zi=yi−fsif′siu+2u21−βu3,whereu=fyifsi.Table 1
Concrete forms of methods.
S. noβΓu,u=fyi/fsiB1ℜΓu=u+2u2+βu4B2ℜΓu=u/1−2u+βu2B3ℜΓu=u+2u2+βu3B4ℜΓu=u/1−βu2+2u2B5ℜΓu=u+2u2/1−βu3


## 2.1. The Concrete Fourth-Order Family of Methods

We now construct some concrete forms of the family of fourth-order methods from the family of methods (12) by choosing weight function Γu containing an arbitrary real number β as provided in Table 1, satisfying the condition Γ0=0,Γ′0=1,Γ″0=4andΓ‴0<∞ of Theorem 1 with β as a real number. Therefore, we get the following new five families of iterative methods:Method-1 (abbreviated as BB1):(20)yi=si−fsif′si,zi=yi−fsif′siu+2u2+βu4.Method-2 (abbreviated as BB2):(21)yi=si−fsif′si,zi=yi−fsif′siu1−2u+βu2.Method-3 (abbreviated as BB3):(22)yi=si−fsif′si,zi=yi−fsif′siu++2u2+βu3.Method-4 (abbreviated as BB4):(23)yi=si−fsif′si,zi=yi−fsif′siu1−βu2+2u2.Method-5 (abbreviated as BB5):(24)yi=si−fsif′si,zi=yi−fsif′siu+2u21−βu3,whereu=fyifsi.Table 1
Concrete forms of methods.
S. noβΓu,u=fyi/fsiB1ℜΓu=u+2u2+βu4B2ℜΓu=u/1−2u+βu2B3ℜΓu=u+2u2+βu3B4ℜΓu=u/1−βu2+2u2B5ℜΓu=u+2u2/1−βu3


## 3. Construction of Simultaneous Methods

Suppose nonlinear equation (1) has n roots. Then, fs and f′s can be approximated as(25)fs=∏j=1ns−sjandf′s=∑k=1n∏j≠kj=1ns−sj.This implies that(26)f′sfs=∑j=1n1s−sj=11/s−si−∑j≠ij=1n1/s−sj.Now, an approximation offsi/f′si is formed by replacing sj with ztj as follows:(27)fsif′si=11/Nsi−∑j≠ij=1n1/si−ztj.Using (27) in (4), we have(28)si+1=si−11/Nsi−∑j≠ij=1n1/si−ztj,fort=1,…,5.Using correctionsztj from BB1 to BB5, we get the following five simultaneous iterative methods for extracting all distinct as well as multiple roots of nonlinear equation (1):(29)si+1=si−σiσi/Nsi−∑j≠ij=1nσj/si−ztj,t=1,…,5,wherez1j=yj−fsjf′sju+2u2+βu4,z2j=yj−fsjf′sju1−2u+βu2,z3j=yj−fsjf′sju+2u2+βu3,z4j=yj−fsjf′sju1−βu2+2u2,z5j=yj−fsjf′sju+2u21−βu3,whereyj=sj−fsjf′sj,u1=fyjfsj.Thus, we constructed new five simultaneous iterative methods (29) abbreviated as M1–M5.


### 3.1. Convergence Analysis

In this section, the convergence analysis of a family of simultaneous methods (M1–M5) is given in the form of the following theorem. Obviously, convergence for the methods (28) will follow from the convergence of the methods (29) from theorem (2) when the multiplicities of the roots are one.Theorem 2.
Letζ1,ζ2,…,ζn be the simple roots of nonlinear equation (1). If s10,s20,s30,…,sn0. If be the initial approximations of the roots, respectively, and sufficiently close to actual roots, the order of convergence of method (58–62) is six.Proof.
Letεi=si−ζiand εi′=si+1−ζi be the errors in approximations si and si+1, respectively. Then, for distinct roots, we have(30)1Nsi=f′sifsi=∑j=1n1si−ζj=1si−ζi+∑j≠ij=1n1si−ζj.
Thus, for multiple roots, we have from (29)(31)si+1=si−σiσi/si−ζi+∑j≠ij=1nσj/si−ζj−∑j≠ij=1nσj/si−ztj,si+1−ζi=si−ζi−σiσi/si−ζi+∑j≠ij=1nσjsi−ztj−si+ζj/si−ζjsi−ztj,this implies,εi′=εi−σiεiσi+εi∑j≠ij=1n−σjztj−ζj/si−ζjsi−ztj,=εi−σiεiσi+εi∑j≠ij=1nEiεj4,whereztj−ζj=εj4Ei=−σjsi−ζjsi−ztj.Thus,εi′=εi2∑j≠ij=1nEiεj4σi+εi∑j≠ij=1nEiεj4.
If it is assumed that absolute values of all errorsεjj=1,2,3,… are of the same order as, say, εj=Oε, then from (30), we have(32)εi′=Oε6.
Thus, (32) shows the convergence order of methods M1–M5 which is six. Hence, the theorem is proved.


## 3.1. Convergence Analysis

In this section, the convergence analysis of a family of simultaneous methods (M1–M5) is given in the form of the following theorem. Obviously, convergence for the methods (28) will follow from the convergence of the methods (29) from theorem (2) when the multiplicities of the roots are one.Theorem 2.
Letζ1,ζ2,…,ζn be the simple roots of nonlinear equation (1). If s10,s20,s30,…,sn0. If be the initial approximations of the roots, respectively, and sufficiently close to actual roots, the order of convergence of method (58–62) is six.Proof.
Letεi=si−ζiand εi′=si+1−ζi be the errors in approximations si and si+1, respectively. Then, for distinct roots, we have(30)1Nsi=f′sifsi=∑j=1n1si−ζj=1si−ζi+∑j≠ij=1n1si−ζj.
Thus, for multiple roots, we have from (29)(31)si+1=si−σiσi/si−ζi+∑j≠ij=1nσj/si−ζj−∑j≠ij=1nσj/si−ztj,si+1−ζi=si−ζi−σiσi/si−ζi+∑j≠ij=1nσjsi−ztj−si+ζj/si−ζjsi−ztj,this implies,εi′=εi−σiεiσi+εi∑j≠ij=1n−σjztj−ζj/si−ζjsi−ztj,=εi−σiεiσi+εi∑j≠ij=1nEiεj4,whereztj−ζj=εj4Ei=−σjsi−ζjsi−ztj.Thus,εi′=εi2∑j≠ij=1nEiεj4σi+εi∑j≠ij=1nEiεj4.
If it is assumed that absolute values of all errorsεjj=1,2,3,… are of the same order as, say, εj=Oε, then from (30), we have(32)εi′=Oε6.
Thus, (32) shows the convergence order of methods M1–M5 which is six. Hence, the theorem is proved.


## 4. Complex Dynamical Study of Families of Iterative Methods

Here, we discuss stability of the family of iterative method (BB1) only in the background contexture of complex dynamics. Rational map arising due to iterative method (BB1) is(33)ℜf=y−fsf′su+2u2+βu4,where y=s−fs/f′s and u=fy/fs.Recalling some basic concepts of this theory, detailed information can be found in [2, 4, 6, 8]. Taking a rational function ℜf:C⟶C, where C denotes the Riemann sphere, the orbit s0∈C defines a set such as orbs=s0,ℜfs0,ℜf2s0,…,ℜfms0,…. A point s0∈C is called a fixed point if ℜfs0=s0. In particular, a fixed point s0 is called the strange fixed point if fsr=0 when s0≠sr. A T-periodic point is defined as a point sT∈C satisfying ℜTsT=sT but ℜ′sT≠sT for t<T. If s0 is the fixed point of ℜ, then it is(i)
Superattracting ifℜf′s0=0(ii)
Attracting ifℜf′s0<1(iii)
Repulsive ifℜf′s0>1(iv)
Neutral ifℜf′s0=1(v)
A strange fixed point if it not associated to any root of nonlinear equation (1)An attracting points∗∈C defines the basin of attraction, ℜs∗, as the set of starting points whose orbit tends to s∗.Furthermore, the implementation of the dynamical plane of the rational operator corresponding to iterative methods divides the complex plane into a mesh of values of real part along thex-axis and imaginary on the y-axis. The initial estimates are depicted in color depending on where their orbit converges, and thus, basins of attraction of the corresponding iterative methods are obtained. The scaling theorem allows the suitable change of the coordinate to reduce dynamics of iteration of general maps to study the specific family of iteration of similar maps.Theorem 3.
Let us take an analytic functionfs and Ts=αs+β be an affine map, with α≠0. Take gs=goTs; then, ToℜgoT−1=ℜfs, i.e., ℜg=ℜf affine conjugate by T (Scaling theorem).As iterative method (33), holds scaling theorem and thus allows the dynamic studies of iterative function (BB1) for the polynomialfs=s−as−b, where a,b∈ℜ. One-point iterative method (BB1) has a universal Julia set if a rational map exists which conjugates by Mobius transformation.Theorem 4.
For a rational mapΩgs arising from (33) applied to fs=s−as−b, where a,b∈ℜ,Ωgs is conjugate via Möbius transformation by Mv=v−a/v−b to(34)Ωgs,β=s4−s6−8s5−27s4+A1−5s6+βs5−24s5−47s4+A2,where A1=−48s3+βs−47s2−24s−5 and A2=−48s3−27s2−8s−1.Proof.
Letfs=s−as−b, where a,b∈ℜ. Möbius transformation is given by(35)Mv=v−av−b,with inverseMv−1=vb−av−1,which is considered as the map from C∪∞. Then, we have(36)Ωgs,β=s4−s6−8s5−27s4+A1−5s6+βs5−24s5−47s4+A2,where A1=−48s3+βs−47s2−24s−5 and A2=−48s3−27s2−8s−1.
Similarly, we can get the following conclusions.Theorem 5.
For a rational mapΩgs arising from (BB2) to (BB5) applied tofs=s−as−b, where a,b∈ℜ,Ωgs is conjugate via Möbius transformation by Mv=v−a/v−b to the following:(37)Ωgs,β=s4s2+β+2s+1β+1s2+2s+1,(38)Ωgs,β=s4−s4−6s3−14s2+β−14s−5β−5s4−14s3−14s2−6s−1,(39)Ωgs,β=s4−s6−8s5+β−27s4+B1β−5s6+2β−24s5+B2,where B1=−48s3+6β−47s2+2β−24s−5 and B2=6β−47s4+4β−48s3+β−27s2−8s−1.(40)Ωgs,β=s4−s4−6s3+βs−14s2−14s−5−5s4+β−14s3−14s2−6s−1.Theorem 6.
For a rational mapΩgsarising from (MM1) to (MM4) applied tofs=s−as−b, where a,b∈ℜ,Ωgs is conjugate via Möbius transformation by Mv=v−a/v−b to the following:(41)Ωgs,β=s4s2+β+2s+2β+12β+1s2+β+2s+1,(42)Ωgs,β=s4s2+4s+55s2+4s+1.(43)Ωgs,β=s4s4+6s3+14s2+14s+44s4+14s3+14s2+6s+1,(44)Ωgs,β=s4s2+2s+22s2+2s+2.The unions of the respective stability functions of all the strange fixed points of rational maps (36)–(41) are graphed in Figures 1(a)–1(f).The fixed points for rational function (36) are s=1,s=0,ands=∞. For the stability of fixed points of iterative method (36), we calculate Ωg′s,β, i.e.,(45)Ωg′s,β=s3s+18A3A42,where A3=−20s4+5β−80s3+−8α2−120s2+5β+80s−20andA4=−5s6+β−24s5−47s4−48s3−27s2−8s−1.It is evident from (36) that s=0 and ∞ are always superattractive fixed points, but stability of other fixed points depends on the value of parameter β which is present here. The operator Ωg′s,β for s=1 gives(46)Ωg′1,β=512β−160.Analyzing (46), as β⟶±∞, we obtained horizontal asymptotes in Ωg′1,β=0.In the following result, we present the stability of the strange fixed point fors=1.Figure 1
(a–f) Unions of the respective stability function of all strange fixed points of the rational map (BB1–BB5 and MM1), respectively.
(a)(b)(c)(d)(e)(f)Theorem 7.
The character of the strange points=1 is as follows:(i)
If512/β−160<1, then s=1 is an attractor, and it can be superattractor for β=1(ii)
When512/β−160=1,s=1 is a parabolic point(iii)
Ifβ≠0 and 512/β−160>1, then s=1 is repulsiveProof.
From (45), we have(47)Ωg′1,β=512β−160.
So,(48)512β−160≤1is equivalent to512≤β−160.
Let us considerβ=a+ib as an arbitrary complex number. Then, a−1602+b2≥5122. Therefore, Ωg′1,β≤1iffβ−160≥512. Finally, if β varies, then Ωg′1,β>1 and s=1 is repulsive, except if β=1andβ=160 for which s=1 is not a fixed point.
The functions where we examine stability of iterative method (36) are given as(49)Ωg′1,β=minΩg′1,β,1.


### 4.1. Analysis of Critical Points

The critical points of (36) satisfy Ωg′s,β=0, i.e., s=0,−1 and s=∞ for β≠0, and(50)cr1=116−16+β−β5β−288165−β1810,cr2=116−16+β−β5β−288165+β1810,cr3=β−1616−β5β−288165−β22,cr4=β−1616−β5β−288165+β22,where β1=−224β+5β2−46085β/5β−288+2685β2/3/5β−288−55β5/2/5β−288, and β2=−2+16−β2/32−30+2β/5+β•/5β−288, where β•=25β−32−16−β2/64+216−β15+β/5. We observe cr1=1/cr2 and cr3=1/cr4. Figure 2 presents the zones of stability of strange fixed points. Fixed points 1 and 0 are represented by black-dotted lines, while critical points −1,cr1,cr2,cr3,cr4 are represented by black, red, blue, green and orange-dotted lines, respectively (see Figure 3).Figure 2
The zones of stability of strange fixed points ofΩ′1,β.Figure 3
Dynamical behavior of strange fixed points and critical points for0≤β≤20 (a) and −20≤β≤0 (b).
(a)(b)Theorem 8.
The only member of the family of iterative methods whose operator is always conjugated to the rational maps2 is the element corresponding to β=0.Proof.
From (51), we represent(51)nns=s4−s6−8s5−27s4−48s3+βs−47s2−24s−5,dds=−5s6+βs5−24s5−47s4−48s3−27s2−8s−1.
The unique value ofβ for which nns=dds is 0.


### 4.2. Parametric Planes

Parametric planes are obtained by takingβ over a mesh of 500×500 values in the complex plane in Reβ×Imβ∈−4,4×−4,4. A critical point is taken as the initial approximation. The method is then iterated until it reaches to the maximum iterations or converges to any fixed point. Taking 10−3 as a tolerance used for stopping criteria fsi. The complex value of the parameter in the complex plane is paint in red if the method converges to any of the roots 0 and∞, black in the other case.In Figures4(a)–4(d) and 5(a)–5(h), red color shows the convergence region, and all those values of parameter β show stable behavior, while those values of the parameter which are taken from the black region (divergence region) show the unstable behavior of the iterative map (36). Stable and unstable behavior is shown in Figures 6 and 7, respectively.Figure 4
(a–d) Parametric planes associated withcr1,cr2,cr3, and cr4, respectively.
(a)(b)(c)(d)Figure 5
(a–h) Details of the parametric planes for cr1, cr2, cr3, and cr4.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
(a–e): Dynamical planes for stable behavior of iterative method (18) for different values of β, i.e., β=0,1,−1,160,  and−20, respectively.
(a)(b)(c)(d)(e)Figure 7
(a–e) Dynamical planes for unstable behavior of the iterative map (18) for different values of β, i.e., β=200,−140,140,500,  and−500, respectively.
(a)(b)(c)(d)(e)


### 4.3. Dynamical Planes

The generations of the dynamical planes are similar to the parametric plane. To execute them, the real and imaginary parts of the starting approximation are represented as two axes over a mesh of250×250 in the complex plane. The stopping criteria are the same as in the parametric plane but assign different colors to indicate to which root the method converges and black in the other case.Let us note that the iterative methods BB2–BB5, MM1–MM4, and BB1 satisfy Cayley’s test [28–31] for all parameter values of β. It can be observed from Figure 8 that iterative methods BB1–BB5 and MM1–MM4 verifying Cayley’s test [28] have the same dynamical properties as Newton’s method. Figures 8(a)–8(e) clearly show quite stable behavior of iterative methods BB1–BB5 over MM1–MM4.Figure 8
(a–i) Dynamical planes of iterative methods (BB1–BB5, MM2, MM4, MM3, and MM1) forβ = 1. (j) Dynamical plane of method MM1 for β = 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)For unstable behavior, the value of parameterβ is chosen from the black region. To generate basins of attractions, we take a square box of −1,1×−1,1∈ℂ. To each root of (1), we assign color to which the corresponding orbit of iterative methods (BB1–BB5 and MM1–MM4) starts and convergences to a fixed point. Take color map as HSV. We use fsi<10−3 as the stopping criteria, and the maximum iteration is taken as 30. We make dark red points if the orbit of the iterative methods does not converge to a root after 30 iterations. Different colors are used for different roots. Iterative methods have different basins of attractions distinguished by their color. In basins, brightness in color represents less number of iterations to achieve the roots of (1). Figures 9 and 10 show the basins of attractions of iterative methods (BB1–BB5 and MM1–MM4) for nonlinear function f1s=s4+4s3−24s2+16s+16 and f2s=sins−1/2sins−2/2sins−2.5/2, respectively. From Figures 9 and 10, divergent regions and brightness in color present that BB1–BB5 are better than MM1–MM4, respectively.Figure 9
(a–i) Basins of attraction of methods BB1–BB5 and MM1–MM4 for the nonlinear equationf1s, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
(a–i) Basins of attraction of methods BB1–BB5 and MM1–MM4 for the nonlinear equationf2s, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 4.1. Analysis of Critical Points

The critical points of (36) satisfy Ωg′s,β=0, i.e., s=0,−1 and s=∞ for β≠0, and(50)cr1=116−16+β−β5β−288165−β1810,cr2=116−16+β−β5β−288165+β1810,cr3=β−1616−β5β−288165−β22,cr4=β−1616−β5β−288165+β22,where β1=−224β+5β2−46085β/5β−288+2685β2/3/5β−288−55β5/2/5β−288, and β2=−2+16−β2/32−30+2β/5+β•/5β−288, where β•=25β−32−16−β2/64+216−β15+β/5. We observe cr1=1/cr2 and cr3=1/cr4. Figure 2 presents the zones of stability of strange fixed points. Fixed points 1 and 0 are represented by black-dotted lines, while critical points −1,cr1,cr2,cr3,cr4 are represented by black, red, blue, green and orange-dotted lines, respectively (see Figure 3).Figure 2
The zones of stability of strange fixed points ofΩ′1,β.Figure 3
Dynamical behavior of strange fixed points and critical points for0≤β≤20 (a) and −20≤β≤0 (b).
(a)(b)Theorem 8.
The only member of the family of iterative methods whose operator is always conjugated to the rational maps2 is the element corresponding to β=0.Proof.
From (51), we represent(51)nns=s4−s6−8s5−27s4−48s3+βs−47s2−24s−5,dds=−5s6+βs5−24s5−47s4−48s3−27s2−8s−1.
The unique value ofβ for which nns=dds is 0.


## 4.2. Parametric Planes

Parametric planes are obtained by takingβ over a mesh of 500×500 values in the complex plane in Reβ×Imβ∈−4,4×−4,4. A critical point is taken as the initial approximation. The method is then iterated until it reaches to the maximum iterations or converges to any fixed point. Taking 10−3 as a tolerance used for stopping criteria fsi. The complex value of the parameter in the complex plane is paint in red if the method converges to any of the roots 0 and∞, black in the other case.In Figures4(a)–4(d) and 5(a)–5(h), red color shows the convergence region, and all those values of parameter β show stable behavior, while those values of the parameter which are taken from the black region (divergence region) show the unstable behavior of the iterative map (36). Stable and unstable behavior is shown in Figures 6 and 7, respectively.Figure 4
(a–d) Parametric planes associated withcr1,cr2,cr3, and cr4, respectively.
(a)(b)(c)(d)Figure 5
(a–h) Details of the parametric planes for cr1, cr2, cr3, and cr4.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 6
(a–e): Dynamical planes for stable behavior of iterative method (18) for different values of β, i.e., β=0,1,−1,160,  and−20, respectively.
(a)(b)(c)(d)(e)Figure 7
(a–e) Dynamical planes for unstable behavior of the iterative map (18) for different values of β, i.e., β=200,−140,140,500,  and−500, respectively.
(a)(b)(c)(d)(e)


## 4.3. Dynamical Planes

The generations of the dynamical planes are similar to the parametric plane. To execute them, the real and imaginary parts of the starting approximation are represented as two axes over a mesh of250×250 in the complex plane. The stopping criteria are the same as in the parametric plane but assign different colors to indicate to which root the method converges and black in the other case.Let us note that the iterative methods BB2–BB5, MM1–MM4, and BB1 satisfy Cayley’s test [28–31] for all parameter values of β. It can be observed from Figure 8 that iterative methods BB1–BB5 and MM1–MM4 verifying Cayley’s test [28] have the same dynamical properties as Newton’s method. Figures 8(a)–8(e) clearly show quite stable behavior of iterative methods BB1–BB5 over MM1–MM4.Figure 8
(a–i) Dynamical planes of iterative methods (BB1–BB5, MM2, MM4, MM3, and MM1) forβ = 1. (j) Dynamical plane of method MM1 for β = 3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)For unstable behavior, the value of parameterβ is chosen from the black region. To generate basins of attractions, we take a square box of −1,1×−1,1∈ℂ. To each root of (1), we assign color to which the corresponding orbit of iterative methods (BB1–BB5 and MM1–MM4) starts and convergences to a fixed point. Take color map as HSV. We use fsi<10−3 as the stopping criteria, and the maximum iteration is taken as 30. We make dark red points if the orbit of the iterative methods does not converge to a root after 30 iterations. Different colors are used for different roots. Iterative methods have different basins of attractions distinguished by their color. In basins, brightness in color represents less number of iterations to achieve the roots of (1). Figures 9 and 10 show the basins of attractions of iterative methods (BB1–BB5 and MM1–MM4) for nonlinear function f1s=s4+4s3−24s2+16s+16 and f2s=sins−1/2sins−2/2sins−2.5/2, respectively. From Figures 9 and 10, divergent regions and brightness in color present that BB1–BB5 are better than MM1–MM4, respectively.Figure 9
(a–i) Basins of attraction of methods BB1–BB5 and MM1–MM4 for the nonlinear equationf1s, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 10
(a–i) Basins of attraction of methods BB1–BB5 and MM1–MM4 for the nonlinear equationf2s, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5. Computational Aspect

Here, we compare the computational efficiency of M. S. Petkovi´c method [32] and the new methods (M1–M5) given by (28). As presented in [32], the efficiency of an iterative method can be estimated using the efficiency index given by(52)ELm=loguD,where D is the computational cost and u is the order of convergence of the iterative method. Using arithmetic operation per iteration with certain weight depending on the execution time of operation, the computational cost D is evaluated. The weights used for division, multiplication, and addition plus subtraction are wd,wm,andwas, respectively. For a given polynomial of degree m and n roots, the number of division, multiplication, addition and subtraction per iteration for all roots is denoted by Dm, Mm, and ASm. The cost of computation can be calculated as(53)D=Dm=wasASm+wmMm+wdDm.Thus, (52) becomes(54)ELm=logrwasASm+wmMm+wdDm.Reducing the number of operations of a complex polynomial of degreem with real and complex roots to operations of real arithmetic, as given in Table 2. Applying (38) and data given in Table 2, we calculate the percentage ratio ρM1−M5,PJ6M [3] which is given by(55)ρPJ6M,M1−M5=ELPJ6MELM1−M5−1×100,ρM1−M5,PJ6M=ELM1−M5ELPJ6M−1×100,where PJ6M is the Petkovic method [32] of order six. Figures 11(a)–11(d) graphically illustrate these percentage ratios. It is evident from Figures 11(a)–11(d) that the newly constructed simultaneous methods (M1–M5) are more efficient as compared with the Petkovic method [32].Table 2
The number of basic operations.
MethodsCOASmMmDmM167m2+Om5m2+Om2m2+OmM267m2+Om2m2+Om2m2+OmM367m2+Om4m2+Om2m2+OmM467m2+Om4m2+Om2m2+OmM567m2+Om4m2+Om2m2+OmPJ6M68m2+Om6m2+Om2m2+OmFigure 11
(a–c) Computational efficiency of methods M1, M3–M5, and M2 w.r.t method PJ6M. (d) Computational efficiency of method PJ6M w.r.t M1–M5.
(a)(b)(c)(d)We also calculate the CPU execution time as all the calculations are done using Maple 18 on Processor Intel(R) Core(TM) i3-3110m CPU@2.4 GHz with 64 bit Operating System. We observe that CPU time of the methods M1–M5 is less than M. S. Petkovi´c methods [32], showing the dominance efficiency of our methods (M1–M5) as compared to them.


## 6. Numerical Results

Here, some numerical examples are considered in order to demonstrate the performance of our family of two-step fourth-order single-root finding methods (BB1–BB5) and six-order simultaneous methods (M1–M5), respectively. We compare our family of optimal fourth-order single-root finding methods (BB1–BB5) with MM1–MM4 methods. Family of simultaneous methods (M1–M5) of convergence order six is compared with M. S. Petkovi´c method [32] of the same order (abbreviated as the PJ6M method). All the computations are performed using CAS Maple 18 with 2500 (64-digit floating-point arithmetic in case of simultaneous methods) significant digits. For single-root finding methods, the stopping criteria are as follows:(i)
fsi<∈,(ii)
sik−α<∈,whereas ei=fsi2 for simultaneous methods.We take∈=10−600 for the single-root finding method and ∈=10−30 for simultaneous determination of all roots of nonlinear equation (1).Numerical test examples from [22, 33–35] are provided in Tables 3–12. In Tables 3, 4, 6, 7, 9, and 11, we present the numerical results for simultaneous determination of all roots while Tables 5, 8, 10, and 12 for single-root finding methods. In all tables, CO represents the convergence order, n, the number of iterations, ρ, the computational order of convergence, and CPU, the computational time in seconds. The value of arbitrary parameter β used in iterative methods BB1–BB5 and M1–M5 is 1.5 for test Examples 1–4. We observe that numerical results of single-root finding methods (BB1–BB5) as well as all-root finding methods (M1–M5) are better than MM1–MM4 and PJ6M, respectively, on the same number of iterations. Figures 12(a)–12(d) represent the residual fall for the iterative methods (M1–M5 and PJ6M), and Figures 12(e)– 12(h) represent the residual fall of iterative methods (BB1–BB5 and MM1–MM4) for given nonlinear functions. Tables 3–12 show that the family of methods BB1–BB5 and M1–M5 are more efficient as compared with MM1–MM4 and PJ6M, respectively.Table 3
Simultaneous finding of all district roots.
f1−1s=s4+4s3−24s2+16s+16ζ1,2=2,ζ3=−4−23,ζ3=2,ζ4=−4+23.s1,20=1.17,s30=−7.4641,s04=−0.5359MethodCOCPUne1e2e3e4PJ6M60.03143.5e − 23.5e − 26.0e − 110.0M1–M560.01545.0e − 135.0e − 130.00.0Table 4
Simultaneous findings of all multiple roots.
f1s=s4+4s3−24s2+16s+16=s−22s2+8s+4ζ1,2=2,ζ3=−4−23,ζ3=2,ζ4=−4+23.s1,20=1.17,s30=−7.4641,s04=−0.5359MethodsCOCPUne1e2e3e4PJ6M60.03243.1e − 103.1e − 103.2e − 60.0M1–M560.01640.00.00.00.0Table 5
Comparison of optimal fourth-order methods.
Methodsi6−αfi6siCPUpf1−1s=s4+4s3−24s2+16s+16,s0=−0.5BB11.4e − 13993.5e − 55940.0324.00BB21.9e − 18251.7e − 72990.0634.00BB31.8e − 18291.5e − 73150.0474.01BB46.4e − 14209.9e − 56760.0314.00BB58.9e − 14203.7e − 56750.0474.01MM14.4e − 12955.3e − 51760.0164.00MM24.9e − 15691.3e − 62720.0314.00MM31.7e − 14584.7e − 58300.0634.01MM41.0e − 14584.0e − 58300.0324.00Table 6
Simultaneous finding of all distinct roots.
f2−1s=sins−1/2sins−2/2sins−2.5/2ζ1=1,ζ2=2,ζ3=2.5,s10=−0.2,s20=1.7,s03=3MethodCOCPUne1e2e3PJ6M60.04740.0194.6e − 39.9e − 10M1–M560.04246.9e − 133.0e − 50.0Table 7
Simultaneous finding of all multiple roots.
f2s=sin3s−1/2sin3s−2/2sin3s−2.5/2ζ1=1,ζ2=2,ζ3=2.5, each with multiplicity 3s10=−0.2,s20=1.7,s30=3MethodCOCPUne1e2e3PJ6M60.09445.0e − 51.7e − 23.3e − 5M1–M560.06240.00.00.0Table 8
Comparison of optimal fourth-order methods.
Methodssi6−αfi6siCPUρf2−1s=sin3s−1/2sin3s−2/2sin3s−2.5/2,s00.5BB15.0e−5734.8e−20031.5164.00BB23.0e−5018.4e−22911.1254.00BB31.5e−5726.8e−22891.1254.00BB48.8e−5062.9e−20211.5314.00BB51.2e−5051.3e−20201.4094.00MM14.9e−4625.1e−18461.4064.00MM21.0e−5632.7e−21451.6574.00MM39.7e−5163.5e−20611.4214.00MM49.7e−5163.5e−0601.5464.00Table 9
Simultaneous finding of all distinct roots.
f3s=2.5/120∗50000∗30000∗600−5s4+2160000s2−6004ζ1=−600,ζ2=−268.328,ζ3=268.328,ζ4=600.s10=−590,s20=−250,s30=255s40=590MethodCOCPUne1e2e3e4PJ6M60.03142.8e − 32.3e − 61.2e − 67.2e − 8M1–M560.01540.01.1e − 301.1e − 300.0Table 10
Comparison of optimal fourth-order methods.
Methodsi6−αfi6siCPUρf3s=2.5/120∗50000∗30000∗600−5s4+2160000s2−6004,s=−590BB15.7e−10503.1e−42080.0314.00BB23.1e−13321.4e−53380.0404.00BB33.6e−13422.5e−53780.0624.00BB41.5e−10601.3e−42500.0454.00BB51.0e−10592.6e−42470.0314.00MM15.4e−9244.7e−37040.0474.00MM21.8e−12296.1e−49261.7344.00MM39.3e−11071.4e−44350.0464.00MM49.3e−11071.4e−44350.0324.00Table 11
Simultaneous finding of all distinct roots.
f4s=s4+11.50s3+47.49s2+83.06325s+51.23266875ζ1=−1.45,ζ2=−2.85,ζ3=−2.85,ζ4=−4.45s10=−1.0,s20=−2.1,s30=−1.8,s40=−3.9MethodCOCPUne 1e 2e 3e 4PJ6M60.04739.0e − 32.0e − 30.42.0e − 3M1–M560.03230.03.6e − 72.4e − 70.0Table 12
Comparison of optimal fourth-order methods.
Methodsi6−αfi6siCPUρf3s=s4+11.50s3+47.49s2+83.06325s+51.23266875,s0=−4.33BB165e − 12393.2e − 49510.0314.00BB21.5e − 15515.1e − 62030.0624.13BB31.3e − 15563.4e − 62230.0474.49BB42.6e − 12561.4e − 50210.0454.49BB55.6e−12571.4e − 50190.0314.01MM17.2e − 11259.1e − 56790.0474.00MM28.9e − 12993.3e − 51920.0354.01MM35.2e − 12998.3e − 51920.0464.11MM45.7e − 12993.3e − 51920.0474.00Figure 12
(a–d) Residual fall of simultaneous iterative methods (M1–M5 and PJ6M). (e–h) Residual fall of iterative methods (BB1–BB5 and MM1–MM4) for nonlinear functionf1s, f2s, f3s, and f4s, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)We also calculate the CPU execution time as all the calculations are done using Maple 18 on Processor Intel (R) Core (TM) i3-3110m CPU@2.4GHz with 64 bit Operating System. We observe from the tables that CPU time of the methods M1–M5 is better than the PJ6M method, showing the efficiency of our methods (M1–M5) as compared to them.


### 6.1. Applications in Engineering

In this section, we discuss the applications in engineering.Example 1.
(for the beam designing model, see [33]).
Consider a problem of beam positioning, resulting in a nonlinear function which is given as(56)f1s=s4+4s3−24s2+16s+16=s−22s2+8s+4.
The exact roots off1s are ζ1,2=2,ζ3=−4−23,ζ3=2,ζ4=−4+23 as shown in Figure 13.
The initial estimates forf3s are taken as(57)s1,20=1.17,s30=−7.4641,s04=−0.5359.
For distinct roots, we used method (28), and for multiple roots, we used method (29).(58)f1−1s=s4+4s3−24s2+16s+16.Figure 13
Exact roots of nonlinear equationf1s.Example 2.
(for determination of all distinct and multiple roots, see [22]). Here, we consider another standard test function for the demonstration of convergence behavior of newly constructed methods.
Consider(59)f2s=sin3s−12sin3s−22sin3s−2.52,with exact roots ζ1=1,ζ2=2,ζ3=2.5 as shown in Figure 14.
The initial guessed values have been taken ass10=−0.2,s20=1.7,s03=3.
For distinct roots, we used method (28), and for multiple roots, we used method (29).(60)f2−1s=sins−12sins−22sins−2.52.Figure 14
Exact roots of nonlinear equationf2s.Example 3.
(for uniform beam design, see [34]).
Figure15 shows a uniform beam subject to a linearly increasing distributed load. The nonlinear equation for the resulting elastic curve is(61)fs=w0120EIL−s5+2L2s3−L4s.
We have to determine the point of maximum deflection which is the value ofs where f′s=0, i.e.,(62)w0120EIL−5s4+6L2s3−L4=0orw0120EIL−5s4+6L2s3−L4where L=600cm,E=50,000KN/cm3,I=30,000cm4,andw0=2.5KN/cm.
The exact roots of (62) are ζ1=−600,ζ2=−268.328,ζ3=268.328,andζ4=600 as shown in Figure 16. The initial estimates for f3s have been taken as(63)s1−0=−590,s20=−250,s30=255,s40=590.Figure 15
(a) Uniform beam design subject to a linearly increasing distributed load. (b) Resulting elastic curve.
(a)(b)Figure 16
Exact roots of nonlinear equationfs.Example 4.
(for continuous stirred tank reactor (CSTR), see [35]). We consider the isothermal stirred tank reactor (CSTR). Items A and R are fed to the reactor at rates of Q and q-Q, respectively. The following complex reaction system develops in the reactor (see [35]):(64)A+R⟶B,B+R⟶C,C+R⟶D,C+R⟶E.
This problem was tested by Douglas (see [36]) in order to construct the simple feedback control systems. In his searching of this system, he designed the following equation for the transfer function of the reactor with a proportional control system:(65)Hc2.98s+2.25s+1.45s+2.852s+4.35=−1,with Hc being the gain of the proportional controller. This control system is stable for values of Hc that yield roots of the transfer function having negative real parts. If we take Hc=0, we get the nonlinear equation:(66)f4s=s4+11.50s3+47.49s2+83.06325s+51.23266875=0.
The transfer function has four negative real roots as shown in Figure17 that are(67)s1=−1.45,s2=−2.85,s3=−2.85,s4=−4.45.
These roots are known as poles of the open-loop transfer function.
The initial guessed values have been taken as(68)s10=−1.0,s20=−2.1,s30=−1.8,s30=−3.9.Figure 17
Exact roots of nonlinear equationf4s.


## 6.1. Applications in Engineering

In this section, we discuss the applications in engineering.Example 1.
(for the beam designing model, see [33]).
Consider a problem of beam positioning, resulting in a nonlinear function which is given as(56)f1s=s4+4s3−24s2+16s+16=s−22s2+8s+4.
The exact roots off1s are ζ1,2=2,ζ3=−4−23,ζ3=2,ζ4=−4+23 as shown in Figure 13.
The initial estimates forf3s are taken as(57)s1,20=1.17,s30=−7.4641,s04=−0.5359.
For distinct roots, we used method (28), and for multiple roots, we used method (29).(58)f1−1s=s4+4s3−24s2+16s+16.Figure 13
Exact roots of nonlinear equationf1s.Example 2.
(for determination of all distinct and multiple roots, see [22]). Here, we consider another standard test function for the demonstration of convergence behavior of newly constructed methods.
Consider(59)f2s=sin3s−12sin3s−22sin3s−2.52,with exact roots ζ1=1,ζ2=2,ζ3=2.5 as shown in Figure 14.
The initial guessed values have been taken ass10=−0.2,s20=1.7,s03=3.
For distinct roots, we used method (28), and for multiple roots, we used method (29).(60)f2−1s=sins−12sins−22sins−2.52.Figure 14
Exact roots of nonlinear equationf2s.Example 3.
(for uniform beam design, see [34]).
Figure15 shows a uniform beam subject to a linearly increasing distributed load. The nonlinear equation for the resulting elastic curve is(61)fs=w0120EIL−s5+2L2s3−L4s.
We have to determine the point of maximum deflection which is the value ofs where f′s=0, i.e.,(62)w0120EIL−5s4+6L2s3−L4=0orw0120EIL−5s4+6L2s3−L4where L=600cm,E=50,000KN/cm3,I=30,000cm4,andw0=2.5KN/cm.
The exact roots of (62) are ζ1=−600,ζ2=−268.328,ζ3=268.328,andζ4=600 as shown in Figure 16. The initial estimates for f3s have been taken as(63)s1−0=−590,s20=−250,s30=255,s40=590.Figure 15
(a) Uniform beam design subject to a linearly increasing distributed load. (b) Resulting elastic curve.
(a)(b)Figure 16
Exact roots of nonlinear equationfs.Example 4.
(for continuous stirred tank reactor (CSTR), see [35]). We consider the isothermal stirred tank reactor (CSTR). Items A and R are fed to the reactor at rates of Q and q-Q, respectively. The following complex reaction system develops in the reactor (see [35]):(64)A+R⟶B,B+R⟶C,C+R⟶D,C+R⟶E.
This problem was tested by Douglas (see [36]) in order to construct the simple feedback control systems. In his searching of this system, he designed the following equation for the transfer function of the reactor with a proportional control system:(65)Hc2.98s+2.25s+1.45s+2.852s+4.35=−1,with Hc being the gain of the proportional controller. This control system is stable for values of Hc that yield roots of the transfer function having negative real parts. If we take Hc=0, we get the nonlinear equation:(66)f4s=s4+11.50s3+47.49s2+83.06325s+51.23266875=0.
The transfer function has four negative real roots as shown in Figure17 that are(67)s1=−1.45,s2=−2.85,s3=−2.85,s4=−4.45.
These roots are known as poles of the open-loop transfer function.
The initial guessed values have been taken as(68)s10=−1.0,s20=−2.1,s30=−1.8,s30=−3.9.Figure 17
Exact roots of nonlinear equationf4s.


## 7. Conclusion

Here, we have developed five families of two-step single-root finding methods of optimal convergence order four and five simultaneous iterative methods of order six, respectively. From Tables3–10 and Figures 9, 10, and 12, we observe that our methods (BB1–BB5 and M1–M5) are superior in terms of efficiency, stability, CPU time, and residual error as compared with the methods MM1–MM4 and PJ6M, respectively [37].



---

*Source: 3524324-2020-07-27.xml*