Optimal Fourth , Eighth and Sixteenth Order Methods by Using Divided Difference Techniques and Their Basins of Attraction and Its Application

Abstract: The principal objective of this work is to propose a fourth, eighth and sixteenth order scheme for solving a nonlinear equation. In terms of computational cost, per iteration, the fourth order method uses two evaluations of the function and one evaluation of the first derivative; the eighth order method uses three evaluations of the function and one evaluation of the first derivative; and sixteenth order method uses four evaluations of the function and one evaluation of the first derivative. So these all the methods have satisfied the Kung-Traub optimality conjecture. In addition, the theoretical convergence properties of our schemes are fully explored with the help of the main theorem that demonstrates the convergence order. The performance and effectiveness of our optimal iteration functions are compared with the existing competitors on some standard academic problems. The conjugacy maps of the presented method and other existing eighth order methods are discussed, and their basins of attraction are also given to demonstrate their dynamical behavior in the complex plane. We apply the new scheme to find the optimal launch angle in a projectile motion problem and Planck’s radiation law problem as an application.


Introduction
One of the most frequent problems in engineering, scientific computing and applied mathematics, in general, is the problem of solving a nonlinear equation f (x) = 0.In most of the cases, whenever real problems are faced, such as weather forecasting, accurate positioning of satellite systems in the desired orbit, measurement of earthquake magnitudes and other high-level engineering problems, only approximate solutions may get resolved.However, only in rare cases, it is possible to solve the governing equations exactly.The most familiar method of solving non linear equation is Newton's iteration method.The local order of convergence of Newton's method is two and it is an optimal method with two function evaluations per iterative step.
In the past decade, several higher order iterative methods have been developed and analyzed for solving nonlinear equations that improve classical methods such as Newton's method, Chebyshev method, Halley's iteration method, etc.As the order of convergence increases, so does the number of function evaluations per step.Hence, a new index to determine the efficiency called the efficiency index is introduced in [1] to measure the balance between these quantities.Kung-Traub [2] conjectured that the order of convergence of any multi-point without memory method with d function evaluations cannot exceed the bound 2 d−1 , the optimal order.Thus the optimal order for three evaluations per iteration would be four, four evaluations per iteration would be eight, and so on.Recently, some fourth and eighth order optimal iterative methods have been developed (see [3][4][5][6][7][8][9][10][11][12][13][14] and references therein).A more extensive list of references as well as a survey on the progress made in the class of multi-point methods is found in the recent book by Petkovic et al. [11].
This paper is organized as follows.An optimal fourth, eighth and sixteenth order methods are developed by using divided difference techniques in Section 2. In Section 3, convergence order is analyzed.In Section 4, tested numerical examples to compare the proposed methods with other known optimal methods.The problem of Projectile motion is discussed in Section 5 where the presented methods are applied on this problem with some existing ones.In Section 6, we obtain the conjugacy maps of these methods to make a comparison from dynamical point of view.In Section 7, the proposed methods are studied in the complex plane using basins of attraction.Section 8 gives concluding remarks.

Design of an Optimal Fourth, Eighth and Sixteenth Order Methods
Definition 1 ([15]).If the sequence {x n } tends to a limit x * in such a way that lim n→∞ x n+1 − x * (x n − x * ) p = C for p ≥ 1, then the order of convergence of the sequence is said to be p, and C is known as the asymptotic error constant.If p = 1, p = 2 or p = 3, the convergence is said to be linear, quadratic or cubic, respectively.Let e n = x n − x * , then the relation ( is called the error equation.The value of p is called the order of convergence of the method.

Definition 2 ([1]
).The Efficiency Index is given by where d is the total number of new function evaluations (the values of f and its derivatives) per iteration.
Let x n+1 = ψ(x n ) define an Iterative Function (IF).Let x n+1 be determined by new information at x n , φ 1 (x n ), ..., φ i (x n ), i ≥ 1.No old information is reused.Thus, Then ψ is called a multipoint IF without memory.The Newton (also called Newton-Raphson) IF (2 nd NR) is given by The 2 nd NR IF is one-point IF with two function evaluations and it satisfies the Kung-Traub conjecture with d = 2. Further, EI 2 nd NR = 1.414.

An Optimal Fourth Order Method
We attempt to get a new optimal fourth order IF as follows, let us consider two step Newton's method The above one is having fourth order convergence with four function evaluations.But, this is not an optimal method.To get an optimal, need to reduce a function and preserve the same convergence order, and so we estimate f (ψ 2 nd NR (x)) by the following polynomial On implementing the above conditions on Equation (6), we obtain three unknowns a 0 , a 1 and a 2 .Let us define the divided differences From conditions, we get , respectively, by using divided difference techniques.Now, we have the estimation Finally, we propose a new optimal fourth order method as The efficiency of the method ( 7) is EI 4thYM = 1.587.

An Optimal Eighth Order Method
Next, we attempt to get a new optimal eighth order IF as following way The above one is having eighth order convergence with five function evaluations.But, this is not an optimal method.To get an optimal, need to reduce a function and preserve the same convergence order, and so we estimate f (ψ 4 th YM (x)) by the following polynomial which satisfies On implementing the above conditions on (8), we obtain four linear equations with four unknowns b 0 , b 1 , b 2 and b 3 .From conditions, we get b 0 = f (x) and b 1 = f (x).To find b 2 and b 3 , we solve the following equations: Thus by applying divided differences, the above equations simplifies to Solving Equations ( 9) and ( 14), we have Further, using Equation ( 11), we have the estimation Finally, we propose a new optimal eighth order method as The efficiency of the method ( 12) is EI 8thYM = 1.682.Remark that the method is seems a particular case of the method of Khan et al. [16], they used weight function to develop their methods.Whereas we used finite difference scheme to develop proposed methods.We can say the methods 4 th YM and 8 th YM are reconstructed of Khan et al. [16] methods.

An Optimal Sixteenth Order Method
Next, we attempt to get a new optimal sixteenth order IF as following way The above one is having eighth order convergence with five function evaluations.However, this is not an optimal method.To get an optimal, need to reduce a function and preserve the same convergence order, and so we estimate f (ψ 8 th YM (x)) by the following polynomial which satisfies On implementing the above conditions on (13), we obtain four linear equations with four unknowns c 0 , c 1 , c 2 and c 3 .From conditions, we get c 0 = f (x) and c 1 = f (x).To find c 2 , c 3 and c 4 , we solve the following equations: , Further, using Equation ( 15), we have the estimation Finally, we propose a new optimal sixteenth order method as The efficiency of the method ( 16) is EI 16thYM = 1.741.

Convergence Analysis
In this section, we prove the convergence analysis of proposed IFs with help of Mathematica software.
Theorem 1.Let f : D ⊂ R → R be a sufficiently smooth function having continuous derivatives.If f (x) has a simple root x * in the open interval D and x 0 chosen in sufficiently small neighborhood of x * , then the method 4 th YM IFs (7) is of local fourth order convergence, and the 8 th YM IFs (12) is of local eighth order convergence. and Thus, Expanding f (ψ 2 nd NR (x)) about x * by Taylor's method, we have Using Equations ( 17)-( 20) in divided difference techniques.We have Substituting Equations ( 18)-( 21) into Equation ( 7), we obtain, after simplifications, Expanding f (ψ 4 th YM (x)) about x * by Taylor's method, we have Now, Substituting Equations ( 19)-( 21), ( 23) and ( 24) into Equation ( 12), we obtain, after simplifications, Hence from Equations ( 22) and ( 25), we concluded that the convergence order of the 4 th YM and 8 th YM are four and eight respectively.
The following theorem is given without proof, which can be worked out with the help of Mathematica.Theorem 2. Let f : D ⊂ R → R be a sufficiently smooth function having continuous derivatives.If f (x) has a simple root x * in the open interval D and x 0 chosen in sufficiently small neighborhood of x * , then the method (16) is of local sixteenth order convergence and and it satisfies the error equation 16 + O(e 17 ).

Numerical Examples
In this section, numerical results on some test functions are compared for the new methods 4 th YM, 8 th YM and 16 th YM with some existing eighth order methods and Newton's method.Numerical computations have been carried out in the MATLAB software with 500 significant digits.We have used the stopping criteria for the iterative process satisfying error = |x N − x N−1 | < , where = 10 −50 and N is the number of iterations required for convergence.The computational order of convergence is given by ( [17]) We consider the following iterative methods for solving nonlinear equations for the purpose of comparison: ψ 4 th SB , a method proposed by Sharma et al. [18]: ψ 4 th CLND , a method proposed by Chun et al. [19]: ψ 4 th SJ , a method proposed by Singh et al. [20]: ψ 8 th KT , a method proposed by Kung-Traub [2]: , , ψ 8 th PNPD , a method proposed by Petkovic et al. [11] , a method proposed by Sharma et al. [12] , a method proposed by Sharma et al. [13] ψ 8 th CFGT , a method proposed by Cordero et al. [6] Table 1 shows the efficiency indices of the new methods with some known methods.

Methods p d EI
The following test functions and their simple zeros for our study are given below [10]: x * = 0.4099920179891371316... Table 2, shows that corresponding results for f 1 − f 6 .We observe that proposed method 4 th YM is converge in a lesser or equal number of iterations and with least error when compare to compared methods.Note that 4 th SB and 4 th SJ methods are getting diverge in f 5 function.Hence, the proposed method 4 th YM can be considered competent enough to existing other compared equivalent methods.
Also, from Tables 3-5 are shows the corresponding results for f 1 − f 6 .The computational order of convergence agrees with the theoretical order of convergence in all the functions.Note that 8 th PNPD method is getting diverge in f 5 function and all other compared methods are converges with least error.Also, function f 1 having least error in 8 th CFGT, function f 2 having least error in 8 th CTV, functions f 3 and f 4 having least error in 8 th YM, function f 5 having least error in 8 th SA2, function f 6 having least error in 8 th CFGT.The proposed 16 th YM method converges less number of iteration with least error in all the tested functions.Hence, the 16 th YM can be considered competent enough to existing other compared equivalent methods.

Projectile Motion Problem
We consider the classical projectile problem [21,22] in which a projectile is launched from a tower of height h > 0, with initial speed v and at an angle θ with respect to the horizontal onto a hill, which is defined by the function ω, called the impact function which is dependent on the horizontal distance, x.We wish to find the optimal launch angle θ m which maximizes the horizontal distance.In our calculations, we neglect air resistances.
The path function y = P(x) that describes the motion of the projectile is given by When the projectile hits the hill, there is a value of x for which P(x) = ω(x) for each value of x.We wish to find the value of θ that maximize x.
Differentiating Equation (37) implicitly w.r.t.θ, we have Setting dx dθ = 0 in Equation (38), we have or An enveloping parabola is a path that encloses and intersects all possible paths.Henelsmith [23] derived an enveloping parabola by maximizing the height of the projectile fora given horizontal distance x, which will give the path that encloses all possible paths.Let w = tan θ, then Equation (36) becomes Differentiating Equation (41) w.r.t.w and setting y = 0, Henelsmith obtained so that the enveloping parabola defined by The solution to the projectile problem requires first finding x m which satisfies ρ(x) = ω(x) and solving for θ m in Equation ( 40) because we want to find the point at which the enveloping parabola ρ intersects the impact function ω, and then find θ that corresponds to this point on the enveloping parabola.We choose a linear impact function ω(x) = 0.4x with h = 10 and v = 20.We let g = 9.8.Then we apply our IFs starting from x 0 = 30 to solve the non-linear equation whose root is given by x m = 36.102990117..... and Figure 1 shows the intersection of the path function, the enveloping parabola and the linear impact function for this application.The approximate solutions are calculated correct to 500 significant figures.The stopping criterion |x N − x N−1 | < , where = 10 −50 is used.Table 6 shows that proposed method 16 th YM is converging better than other compared methods.Also, we observe that the computational order of convergence agrees with the theoretical order of convergence.

Planck's Radiation Law Problem
We consider the following Planck's radiation law problem found in [24]: which calculates the energy density within an isothermal blackbody.Here, λ is the wavelength of the radiation, T is the absolute temperature of the blackbody, k is Boltzmann's constant, h is the Planck's constant and c is the speed of light.Suppose, we would like to determine wavelength λ which corresponds to maximum energy density ϕ(λ).From (44), we get It can be checked that a maxima for ϕ occurs when B = 0, that is, when (ch/λkT)e ch/λkT e ch/λkT − 1 = 5.
Here putting x = ch/λkT, the above equation becomes The aim is to find a root of the equation f (x) = 0. Obviously, one of the root x = 0 is not taken for discussion.As argued in [24], the left-hand side of (45) is zero for x = 5 and e −5 ≈ 6.74 × 10 −3 .Hence, it is expected that another root of the equation f (x) = 0 might occur near x = 5.The approximate root of the Equation ( 46) is given by x * ≈ 4.96511423174427630369 with x 0 = 3.Consequently, the wavelength of radiation (λ) corresponding to which the energy density is maximum is approximated as λ ≈ ch (kT)4.96511423174427630369 .
Table 7 shows that proposed method 16 th YM is converging better than other compared methods.Also, we observe that the computational order of convergence agrees with the theoretical order of convergence.Hereafter, we will study the optimal fourth and eighth order methods along with Newton's method.

Corresponding Conjugacy Maps for Quadratic Polynomials
In this section, we discuss the rational map R p arising from 2 nd NR, proposed methods 4 th YM and 8 th YM applied to a generic polynomial with simple roots.Theorem 3. (2 nd NR) [18] For a rational map R p (z) arising from Newton's method (4) Theorem 4. (4 th YM) For a rational map R p (z) arising from Proposed Method (7)  (z−1) , which may be considered as map from C ∪ {∞}.We then have Theorem 5. (8 th YM) For a rational map R p (z) arising from Proposed Method (12) , which may be considered as map from C ∪ {∞}.We then have Remark 1.The methods (29)-( 35) are given without proof, which can be worked out with the help of Mathematica.
Remark 2. All the maps obtained above are of the form S(z) = z p R(z), where R(z) is either unity or a rational function and p is the order of the method.

Basins of Attraction
The study of dynamical behavior of the rational function associated to an iterative method gives important information about convergence and stability of the method.The basic definitions and dynamical concepts of rational function which can found in [4,25].
We take a square R × R = [−2, 2] × [−2, 2] of 256 × 256 points and we apply our iterative methods starting in every z (0) in the square.If the sequence generated by the iterative method attempts a zero z * j of the polynomial with a tolerance | f (z (k) )| < ×10 −4 and a maximum of 100 iterations, we decide that z (0) is in the basin of attraction of this zero.If the iterative method starting in z (0) reaches a zero in N iterations (N ≤ 100), then we mark this point z (0) with colors if |z (N) − z * j | < ×10 −4 .If N > 50, we conclude that the starting point has diverged and we assign a dark blue color.Let N D be a number of diverging points and we count the number of starting points which converge in 1, 2, 3, 4, 5 or above 5 iterations.In the following, we describe the basins of attraction for Newton's method and some higher order Newton type methods for finding complex roots of polynomials p Figures 2 and 3 shows the polynomiographs of the methods for the polynomial p 1 (z).We can see that the methods 2 nd NR, 4 th YM, 8 th SA2 and 8 th YM performed very nicely.The methods 4 th SB, 4 th SJ, 8 th KT, 8 th LW, 8 th PNPD, 8 th SA1, 8 th CFGT and 8 th CTV are shows some chaotic behavior near the boundary points.The method 4 th CLND have sensitive to the choice of initial guess in this case.
Figures 2 and 4 shows the polynomiographs of the methods for the polynomial p 2 (z).We can see that the methods 2 nd NR, 4 th YM, 8 th SA2 and 8 th YM performed very nicely.The methods 4 th SB, 8 th KT, 8 th LW and 8 th CTV are shows some chaotic behavior near the boundary points.The methods 4 th CLND, 4 th SJ, 8 th PNPD, 8 th SA1, and 8 th CFGT have sensitive to the choice of initial guess in this case.
Figures 2 and 5 shows the polynomiographs of the methods for the polynomial p 3 (z).We can see that the methods 4 th YM, 8 th SA2 and 8 th YM are shows some chaotic behavior near the boundary points.The methods 2 nd NR, 4 th SB, 4 th CLND, 4 th SJ, 8 th KT, 8 th LW, 8 th PNPD, 8 th SA1, 8 th CFGT and 8 th CTV have sensitive to the choice of initial guess in this case.In Tables 8-10, we classify the number of converging and diverging grid points for each iterative method.We note that a point z 0 belongs to the Julia set if and only if the dynamics in a neighborhood of z 0 displays sensitive dependence on the initial conditions, so that nearby initial conditions lead to wildly different behavior after a number of iterations.For this reason, some of the methods are getting divergent points.The common boundaries of these basins of attraction constitute the Julia set of the iteration function.It is clear that one has to use quantitative measures to distinguish between the methods, since we have a different conclusion when just viewing the basins of attraction.
In order to summarize the results, we have compared mean number of iteration and total number of functional evaluations (TNFE) for each polynomials and each methods in Table 11.The best method based on the comparison in Table 11 is 8 th SA2.The method with the fewest number of functional evaluations per point is 8 th SA2 followed closely by 4 th YM.The fastest method is 8 th SA2 followed closely by 8 th YM.The method with highest number of functional evaluation and slowest method is 8 th PNPD.
Table 11.Mean number of iteration (N µ ) and TNFE for each polynomials and each methods.

IF
N µ f or p 1 (z) N µ f or p 2 (z) N µ f or p 3 (z) Average TNFE

Concluding Remarks and Future Work
In this work, we have developed optimal fourth, eighth and sixteenth order iterative methods for solving nonlinear equations using the divided difference approximation.The methods require the computations of three functions evaluations reaching order of convergence is four, four functions evaluations reaching order of convergence is eight and five functions evaluations reaching order of convergence is sixteen.In the sense of convergence analysis and numerical examples, the Kung-Traub's conjecture is satisfied.We have tested some examples using the proposed schemes and some known schemes, which illustrate the superiority of the proposed method 16 th YM.Also, proposed methods and some existing methods have been applied on the Projectile motion problem and Planck's radiation law problem.The results obtained are interesting and encouraging for the new method 16 th YM.The numerical experiments suggests that the new methods would be valuable alternative for solving nonlinear equations.Finally, we have also compared the basins of attraction of various fourth and eighth order methods in the complex plane.
Future work includes: • Now we are investigating the proposed scheme to develop optimal methods of arbitrarily high order with Newton's method, as in [26].• Also, we are investigating to develop derivative free methods to study dynamical behavior and local convergence, as in [27,28].
Author Contributions: The contributions of both the authors have been similar.Both of them have worked together to develop the present manuscript.
Funding: This paper is supported by three project funds:

Figure 1 .
Figure 1.The enveloping parabola with linear impact function.

Table 1 .
Comparison of Efficiency Indices.

Table 2 .
Numerical results for nonlinear equations.

Table 3 .
Numerical results for nonlinear equations.

Table 4 .
Numerical results for nonlinear equations.

Table 5 .
Numerical results for nonlinear equations.

Table 6 .
Results of projectile problem.

Table 7 .
Results of Planck's radiation law problem.

Table 8 .
Results of the polynomials p 1

Table 9 .
Results of the polynomials p 2

Table 10 .
Results of the polynomials p 3 1. National College Students Innovation and entrepreneurship training program of Ministry of Education of the People's Republic of China in 2017: Internet Animation Company in Minority Areas-Research Model of "Building Dream" Animation Company.(Project number: 201710684001).2. Yunnan Provincial Science and Technology Plan Project University Joint Project 2017: Research on Boolean Satisfiability Dividing and Judging Method Based on Clustering and Partitioning (Project number: 2017FH001-056).3. Qujing Normal college scientific research fund special project (Project number: 2018zx003).