# Effective High-Order Iterative Methods via the Asymptotic Form of the Taylor-Lagrange Remainder

Authors: Isaac Fried
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The asymptotic form of the Taylor-Lagrange remainder is used to derive some new, efficient, high-order methods to iteratively locate the root, simple or multiple, of a nonlinear function. Also derived are superquadratic methods that converge contrarily and superlinear and supercubic methods that converge alternatingly, enabling us not only to approach, but also to bracket the root.


---

## Body

## 1. The Asymptotic Form of the Taylor-Lagrange Remainder

The Taylor-Lagrange theorem is a corollary of extended Rolle’s theorem, which for the spare polynomial function(1)f(x)=xn(1-x),0≤x≤1
such that
(2)f(0)=f′(0)=f′′(0)=⋯=f(n-1)(0)=0,f(n)(0)≠0,f(1)=0
is exactly
(3)f(n)(ξ)=0,0<ξ=1n+1<1.Asn, the degree of osculation at point x=0, increases, ξ inexorably edges leftwards toward that point.The Taylor-Lagrange theorem states that if functionf(t) is continuous in the closed interval [a,x] and is n times differentiable in the open interval (a,x) and f(n-1)(a) exists, then
(4)f(x)=pn-1(x)+1n!(x-a)nf(n)(ξ),a<ξ<x
with the osculating polynomial part of the formula being
(5)pn-1(x)=f(a)+f′(a)(x-a)+12!f′′(a)(x-a)2+⋯+1(n-1)!f(n-1)(a)(x-a)n-1.As with extended Rolle’s theorem, asn increases, ξ moves gradually closer to point a. More precisely,Theorem 1.
If in addition to the required conditions on the function in the Taylor-Lagrange formula, alsof(n+1)(a)≠0 and f(n+1)(t) exists in [a,x] and is continuous from the right at t=a, then ξ in (4) is such that
(6)lim⁡x→a+ξ-ax-a=1n+1
or
(7)ξ=nn+1a+nn+1x
nearly, if x is nearly a, implying that
(8)f(x)=pn-1(x)+1n!(x-a)nf(n)(nn+1a+nn+1x)+R(x)
with
(9)R(x)=O((x-a)k),k≥n+2.We will give here an explicit elementary proof to this theorem forn=2 only; for its generalization, see [1, 2].We write the Taylor-Lagrange formula forn=2, a=0(10)f(x)=f(0)+xf′(0)+12x2f′′(ξ),0<ξ<x
so as to have, with ξ=kx,
(11)f(x)=f(0)+xf′(0)+12x2f′′(kx)+R(x)
and seek k so that R(x)=O(x4).Using L’Hopital’s rule, we determine that(12)lim⁡x→0R(x)x3=16(1-3k)f′′′(0).
We take k=1/3=1/(n+1) and establish that the asymptotic error with this k is
(13)R(x~0)=172x4f′′′′(0).For example(14)ln⁡(1+x)=x-92(x3+x)2+R(x),R(x~0)=-112x4.


## 2. From Newton to Halley via the Asymptotic Remainder

We write then=1 version of the Taylor-Lagrange formula
(15)f(x)=f(x0+δx)=f(x0)+δxf′(ξ),kkkkkkkkkkkkkx0<ξ<x0+δx
and obtain from it, by setting f(x)=0,ξ=x0, the classical Newton method
(16)x1=x0+δx=x0-u0,δx=-u0,u0=f0f0′,f0=f(x0),f0′=f′(x0)
which is, as is well known, a second order method
(17)x1-a=12f′′f′(x0-a)2+O((x0-a)3),
provided that f′=f′(a)≠0. In the above equation x0 is the input and x1 is the output of the iterative process.The difficult problem of finding the root of the nonlinear functionf(x) is replaced in Newton’s method by the easy task of repeatedly finding the approximating root of the tangent line g(x)=f0′(x-x0)+f0. This is the essence of all other higher-order methods, to supplant the finding of the root of the original function by the repeated finding of the root of an approximating polynomial.Having computedx1 by (16), it occurs to us to return and replace the initial ξ=x0 by the asymptotic
(18)ξ=x0+x12=x0+12δx,δx=-u0=-f0f0′
to have the two-step, mid-point method
(19)x2=x0-f(x0)f′(x0+(1/2)δx),δx=-u0=-f0f0′
which is cubic or third order
(20)x2-a=1246f′′2-f′f′′′f′2(x0-a)3+O((x0-a)4).
See also Traub [3, page 164, (8)–(12)].The linearization(21)f′(x0+12δx)=f′(x0)+12δxf′′(x0),δx=-u0=-f(x0)f′(x0)=-f0f0′
reproduces out of (19) classical Halley’s method
(22)x1=x0+det⁡[-1f002f0′]det⁡[f0′f0f0′′2f0′]f0=x0-2f0′2f0′2-f0f0′′f0=x0-11-(f0′′/2f0′)u0u0,
which is cubic as well,
(23)x1-a=1123f′′2-2f′f′′′f′2(x0-a)3+O((x0-a)4),
but requires the second derivative f′′(x).Power series expansion modifies rational Halley’s method into the polynomial inu0 form
(24)x1=x0-u0-12f0′′f0′u02,u0=f0f0′,
which is still cubic,
(25)x1-a=163f′′2-f′f′′′f′2(x0-a)3+O((x0-a)4)
provided that f′=f′(a)≠0. See also Traub [3, page 205, (9)–(44)].We write the equation of the osculating parabola(26)g(x)=f0+f0′(x-x0)+12f0′′(x-x0)2
to f(x) at x0 and seek its intersection with the x-axis. The smaller root, x1, of g(x) is given by
(27)x1-x0=-f0′+f0′2-2f0′′f0f0′′=-u0-12f0′′f0′u02+O(u03),kkkkkkkwkkkkkkkkkkkkkkkkkkkkkkkku0=f0f0′
which is Halley’s method of (24).


## 3. Construction of High-Order Iterations by Generalized Undetermined Coefficients

Halley’s method or, for that matter, any other higher-order method such as that in (24) can be derived ab initio by writing δx, x1=x0+δx as a power series of u0=f0/f0′ or merely f0=f(x0), as in
(28)x1=x0+Pf0+Qf02,
and then progressively fixing the undetermined coefficients P and Q, which eventually need not remain constant, so as to achieve the highest possible order of convergence.Thus, at first, we have from (28) that
(29)x1-a=(1+Pf′)(x0-a)+O((x0-a)2).
We substitute variable f′(x0) for the constant f′(a) in (29) and try P=-1/f0′. With this P we have next that
(30)x1-a=2Qf′3+f′′2f′(x0-a)2+O((x0-a)3)
and we set
(31)P=-1f0′,Q=-f0′′2f0′3
with which the polynomial variant of Halley’s method in (24) is regained.Doing the same to the rational method(32)x1=x0+P+Qf0R+Sf0f0,
we verify that cubic convergence is achieved with P=-1/f0′, Q=0, R=1, and S=-f0′′/(2f0′2) as in classical Halley’s method in (22).For other interesting applications of the method of undetermined coefficients, see [4, 5].


## 4. High-Order Iterative Methods Derived from the Weighted Fixed Point Iteration

Consider the fixed point iteration(33)x1=F(x0)
for point a,F(a)=a. We write x1-a=F(x0)-a and have the power series expansion
(34)x1-a=F′(a)(x0-a)+12!F′′(a)(x0-a)2+13!F′′′(a)(x0-a)3+⋯.
Hence, if |F′(x)|<1 around x=a, then the fixed point iteration converges linearly; and if F′(a)=0, then the fixed point iteration converges quadratically, and so on.Suppose now that we are seeking single roota of f(x),f(a)=0,f′(a)≠0. We rewrite f(x)=0 as the equivalent fixed point problem
(35)x=F(x),F(x)=x+w(x)f(x)
for weight function w(x),w(a)≠0, and seek to fix it to our advantage. For a quadratic method, w(x) needs to be such that
(36)F′(x)=1+w′(x)f(x)+w(x)f′(x)=0
for x near a. Since f(a)=0, we choose to ignore w′(x)f(x) in the previous equation and are left with w(x0)=-1/f′(x0) and Newton’s method.From(37)F′(x)=1+f′(x)w(x)+f(x)w′(x)=0,F′′(x)=f′′(x)w(x)+2f′(x)w′(x)+f(x)w′′(x)=0
and ignoring f(x)w′′(x) in the second equation, we obtain the system
(38)[f′ff′′2f′][ww′]=[-10]
which we solve for w(x) as
(39)w(x0)=det⁡[-1f002f0′]det⁡[f0′f0f0′′2f0′]
and arrive at Halley’s method of (22). Higher-order methods are systematically generated in similar fashion. See also [6].


## 5. The Recursive Generation of the High-Order Iteration Function

LetF2(x)(40)F2(a)=a,F2′(a)=0
be the fixed point iteration function of the recursion x1=F2(x0). By dint of F2′(a)=0, the iterative method is quadratic
(41)x1-a=12F2′′(a)(x0-a)2+O((x0-a)3).The recursively constructed iteration function(42)Fn+1(x)=Fn(x)+1nFn′(x)(Fn(x)-x),n=2
assures a third order convergence
(43)x1-a=112(3F2′′(a)2-F2′′′(a))(x0-a)3+O((x0-a)4).For example, taking(44)F2(x)=x-mu(x),u(x)=f(x)f′(x),
we obtain
(45)F3(x)=x-12m(3-m)u(x)-12m2f′′(x)f′(x)u(x)2,KKKKKKKKKKKKKKKu(x)=f(x)f′(x).
See also Traub [3, (7)–(11)].Indeed, takingf(x) as
(46)f(x)=(x-a)mg(x),g(a)≠0,
we have with F3(x) in (45) that
(47)x1-a=12m2(3+m)g′2-mgg′′g2×(x0-a)3+O((x0-a)4),
where g=g(a), g′=g′(a), g′′=g′′(a).For more on such recursive formulas, see Traub [3, Section 8.3] and Petkovic´ et al. [7, Theorem 2 and Remark 1].


## 6. A Finite-Difference Approximation

Wishing to avoid the possibly computationally costly additional derivative in (19), we propose to approximate it by the central difference scheme
(48)f′(x0+12δx)=f(x0+(1/2)δx+h)-f(x0+(1/2)δx-h)2h+O(h2).
Taking h=-δx/2 leaves us with the approximation
(49)f′(x0+12δx)=f(x0)-f(x0-u0)u0=(1-r)f0′,kkkkkkkkδx=-u0=-f0f0′,r=f1f0,
where x1=x0-u0, f1=f(x1), by which (19) becomes the cubic chord or two-step method
(50)x2=x0-11-ru0,u0=f0f0′,x1=x0-u0,f1=f(x1),r=f1f0.
See also Traub [3, page 180, (8)–(55)]. We return to this method in the next section.The second derivative approximation(51)f(x-u)=f(x)-uf′(x)+12u2f′′(x)=12u2f′′(x),f′′(x)=2u2f1,u=ff′,f1=f(x-u)
leads to the same result.


## 7. A Cubic, One-Sided, Two-Step, Secant Method

Having computed by Newton’s method(52)f0=f(x0),f0′=f′(x0),x1=x0-u0,u0=f0f0′,f1=f(x1),
we propose to proceed and predict the next x2 by pseudo-Newton’s method
(53)x2=x1-f1g1′,x1=x0-u0,g1′=f0-f1x0-x1=(1-r)f0′,
skirting the computation of a new f′(x1). In (53)
(54)r=f1f0=f(x0-u0)f0=12f′′f′(x0-a)+O(x0-a)2,kkkkkkkkkkkkkkkkkkkf′=f′(a),f′′=f′′(a).
We write (53) variously as
(55)x2=x0-11-ru0orx2=x0-(1+r)u0orx2=x0-(1+r+r2)u0,r=f1f0
with all three methods being cubic
(56)x2-a=14(f′′f′)2(x0-a)3,x2-a=12(f′′f′)2(x0-a)3,x2-a=14(f′′f′)2(x0-a)3.
Notice the extra r2 in the last method of (55), added to recover the factor 1/4 in the last of error equations (56).Convergence is here one sided: ifx0-a>0, then x2-a>0; and if x0-a<0, then x2-a<0.For example, forf(x)=x+x2, we obtain by the first method of (55) the two oppositely or contrarily converging sequences
(57)x2={-16,-1126,-11,953,126},x2={14,128,11,953,124}
by which root a=0 is bounded as
(58)-11,953,126<a<11,953,124.Method (55) is also obtained from the secant line
(59)g(x)=Ax+B,A=f1-f0x1-x0,B=x1f0-x0f1x1-x0,
passing through the two points, (x0,f0),(x1=x0-u0,f1=f(x1)), and then taking the root of g(x)=0 as the next
(60)x2=-BA=x1-rx01-r=x0-11-ru0,r=f1f0.Includingf0′ in the polynomial interpolation and passing a parabola through the available data (x0,f0,f0′) and (x1,f1) should allow us to obtain a better approximation for f′(x2), and with it a higher-order method, as we will see next.


## 8. Quartic Two-Step Methods

Seeking a possibly higher-order method, we write the slope estimate of (53) as g1′=(1-Pr)u0 for undetermined coefficient P and then advantageously determine it, to have
(61)x2=x1-f1g1′,g1′=(1-2r)f0′,x1=x0-u0,f1=f(x1),r=f1f0
or
(62)x2=x0-1-r1-2ru0,
which is the celebrated quartic
(63)x2-a=124f′′(3f′′2-2f′f′′′)f′3(x0-a)4+O((x0-a)5)
method of Ostrowski [8]. See also Traub [3, page 184, (8)–(78)] and King [9]. Quartic method (62) is also obtained by replacing f1 by 2f1 in the slope estimate g1′ in (53).The method(64)x2=x1-f1g1′,g1′=(1-2r-r2)f0′,x1=x0-u0,f1=f(x1),r=f1f0
is quartic as well but with the simpler error equation
(65)x1-a=112f′′f′′′f′2(x0-a)4+O((x0-a)5).Power series expansion changes rational method (62) into the polynomial in r method
(66)x2=x0-(1+r+2r2)u0,r=f1f0
which is still quartic
(67)x2-a=124f′′f′3(15f′′2-2f′f′′′)(x0-a)4+O((x0-a)5).The quartic method of (66) is also obtained from the parabola passing through the data (x0,f0,f0′) and (x1,f1)(68)g(x)=Ax2+Bx+C,A=rf0′2f0,B=f0′-2Ax0,kkkkkkkkkkkkkkkkkkkkkC=f0-Bx0-Ax02
with the predicted x2 such that g(x2)=0 or by taking g1′=2Ax1+B in pseudo-Newton’s method x2=x1-f1/g1′.Replacingg1′ in (61) by the perturbed slope
(69)g1′=(1-2r-2kr)f0′,k>0
turns the method into the supercubic alternatingly converging method
(70)x2-a=-12k(f′′f′)2(x0-a)3+O((x0-a)4).


## 9. Quintic, Sextic, Septic, and Octic Three-Step Methods

We continue with higher-order multistep methods requiring only the sole derivativef0′=f′(x0) at initial point x0.The pseudo-Newton method(71)x3=x2-f2g2′,g2′=f1-f2x1-x2,x2=x0-11-ru0,f2=f(x2)
or
(72)x3=x2-f2g2′,g2′=(1-r)(1-s)f0′,r=f1f0,s=f2f1,f1=f(x1),f2=f(x2)
is quintic and one sided
(73)x3-a=116(f′′f′)4(x0-a)5+O((x0-a)6).The method(74)x3=x2-f2g2′,x2=x0-(1+r+2r2)u0,g2′=(1-r)(1-s)f0′,r=f1f0,s=f2f1
is sextic
(75)x3-a=196f′′3(15f′′2-2f′f′′′)f′5(x0-a)6+O((x0-a)7).
The method
(76)x3=x2-f2g2′,x2=x0-11-ru0,g2′=(1-r+r2-s(1+r-2r2))f0′,r=f1f0,s=f2f1
is septic
(77)x3-a=1192f′′3(9f′′3-6f′f′′f′′′+f′2f(4))f′6(x-a)7+O((x0-a)8).The method(78)x3=x2-f2g2′,x2=x0-(1+r+2r2)u0,g2′=(1-2r+3r2-s)f0′,
where r=f1/f0 and s=f2/f1, is octic
(79)x3-a=11152×(f′′2(15f′′2-2f′f′′′)kkkk×(3f′′3-6f′f′′f′′′+f′2f(4)))×(f′7)-1(x0-a)8+O((x0-a)9).See also [7, 10–16].Octic method (78) is also obtained from a cubic polynomial g(x) passing through the data (x0,f0,f0′), (x1,f1), (x2,f2), and x3 chosen such that g(x3)=0.However, if the root repeats, that is, iff′(a)=0, then the order of convergence of the method plummets from eighth order to first order.


## 10. Estimates for the Root Multiplicity Index

In this section we derive both first and second order estimates for the root multiplicity indexm. Also derived is an estimate for the relative size of the second term in the Taylor expansion of function f(x) at root point a.Assuming that the power series expansion of functionf(x), whose root a we are seeking, is of the form
(80)f(x)=(x-a)m(A+B(x-a)+C(x-a)2+⋯),fffffffffffffffffffffFFffm≥1
we obtain from it the first order estimate for the root multiplicity index m(81)u′=1m-2m2BA(x-a)+O((x-a)2),u=ff′,(82)-12u′′u′2=BA+O((x-a))
as well as the second order estimate for m(83)u′2-2uu′′=1m2+6m4-B2(1+m)+2ACmA2(x-a)2+O((x-a)3).For example, forf=x2+x3,  m=2, we compute, at x=0.1, the first order and second order approximations
(84)m=2.18,m=2.03,
respectively. For x=0.1 and x=0.01, we compute from (82)
(85)BA=0.78,BA=0.9753,
respectively.We have also that(86)f(x0-ku0)f(x0)=(1-km)m+O((x0-a)),f(x0-mu0)f(x0)=1mm(BA)3(x0-a)3+O((x0-a)4).From the first of the previous equations we have, fork=1,
(87)ln⁡r=mln⁡(1-1m),r=f(x0-ku0)f(x0),
which, with the Padé rational approximation
(88)ln⁡(1+x)=xx+64x+6+R,R(x~0)=-136x4,
becomes the first order estimate for the multiplicity index
(89)m=161+4ln⁡r1+ln⁡r,r=f(x0-u0)f(x0),u0=f0f0′.For example, forf(x)=x3+x4, x={1.0,0.5,0.1}, (89) yields the approximations m={3.72,3.51,3.14} for the exact m=3.Traub [3, Section 7.8] has the pointwise approximation
(90)m=ln⁡|f(x)|ln⁡|u(x)|
which yields for f(x)=x3+x4, x={0.1,0.01}, the less precise estimates m={1.99,2.42}.Suppose thatm>1. Application of Newton’s method x1=x0-f0/f0′ to find a root of multiplicity greater than one reduces the method to first order
(91)x1-a=(1-1m)(x0-a)+1m2BA(x0-a)2+O((x0-a)3).
Ignoring higher-order terms and eliminating root a from the pair of successive approximations
(92)x1-a=(1-1m)(x0-a),x2-a=(1-1m)(x1-a)
leave us with the discrete, first order approximation
(93)μ=x0-x1u0-u1=x0-x1x0-2x1+x2=m+BA2m-1m(x0-a)+O((x0-a)2)to the multiplicity index m of the root, an approximation which is the discrete counterpart to that in (81).Now(94)x3=μx2-(μ-1)x1orx2=x1-μf1f1′,μ=x0-x1u0-u1
is quadratic
(95)x2-a=-BAm-1m2(x0-a)2+O((x0-a)3)
with no need for prior knowledge of m.For example, using method (94), with the updated μ employed in the computation of x1 in each cycle, we generate the alternatingly converging sequences x2={1.0,-6.5·10-2,4.0·10-4,-9.2·10-10} and m={3.73,2.95,3.0004}. Actually, for a pair of successively computed x values in the sequence, |x1|=0.43|x0|2.6.For a high-order method, to realize its full speed of convergence, it is necessary that the estimatedm is appropriately accurate. For example, using the estimate for m from (90) in the modified Newton method x1 = x0-mf0/f0′, we obtain for f(x)=x2+x3 the sequences m={1.27,1.52,1.68,1.78,1.83} and x={0.2,0.082,0.022,0.0037,0.00042,0.000032}, with convergence, that is, barely above the linear: x1=0.32(x0)p, p=1.184.


## 11. Correction for Multiple Roots by Undetermined Coefficients

We rewrite Newton’s method as(96)x1=x-Pu0,u0=f0f0′
for the undetermined coefficient P and have that near a root of multiplicity m≥1(97)x1-a=(1-Pm)(x0-a)+Pm2BA(x0-a)2+O((x0-a)3).
Quadratic convergence is restored, as is well known, with P=m. In the previous equation, A and B are the coefficients in the power series expansion of f(x) in (80).WithP=m(1-k),k<0, the modified Newton method of (96) is reduced to an alternating superlinear method. For example, for f(x)=x2+x3, m=2, k=-1/8, we generate, starting with x0=1,
(98)x1={1,1.0·10-1,-7.6·10-3,9.8·10-4,kk-1.2·10-4,1.5·10-5,-1.9·10-6}.Next, we rewrite the method in (50) as
(99)x2=x0-PQ-ru0,r=f1f0,f1=f(x1),x1=x0-u0
and seek to adjust correction coefficients P and Q so that convergence remains cubic even to root a of multiplicity m>1. By power series expansion we determine that
(100)P=Q=(m-1m)m-1,m>1,P=Q=1ifm=1
upholds cubic convergence
(101)x2-a=mB2-2(m-1)AC2m2A2(x0-a)3+O((x0-a)4),
where A, B, and C are as in (80). Method (99)-(100) is found in [17]. See also [18].The method(102)x2=x0-(P+Qr)u0,P=m(2-m),Q=mm+1(m-1)m-1,m>1
is similarly cubic
(103)x2-a=(m+2)B2-2(m-1)AC2m2A2(x0-a)3+O((x0-a)4).No such correction to account for multiple roots exists for the quartic two-step method of (62).


## 12. Correction of Halley’s Method for Multiple Roots

We parametrize Halley’s method of (22) with the undetermined coefficients P and Q as
(104)x1=x0-Pf0′Qf0′2-f0f0′′f0
and determine by power series expansion that for
(105)P=2,Q=1+1m
convergence remains cubic even to a root of any multiplicity m≥1(106)x1-a=(m+1)B2-2mAC2m2A2(x0-a)3+O((x0-a)4).
Method (104)-(105) is found in Hansen and Patrick [19].Method (24) becomes here for a multiple root
(107)x1=x0+(P+Qf0′′f0′u0)u0,P=12m(m-3),Q=-12m2
with error equation (106).


## 13. From a Cubic to a Quartic Method by Taylor’s Formula

We write the second order,n=2, version of the Taylor-Lagrange formula
(108)f(x)=f(x0+δx)=f(x0)+δxf′(x0)+12δx2f′′(ξ),KKKKKKKwKKKKKKKKKKx0<ξ<x0+δx
and take f(x1=x0+δx)=0,ξ=x0 to obtain the iterative method
(109)x1=x0+δx,0=f(x0)+δxf′(x0)+12δx2f′′(x0).
We propose to approximate the solution of the quadratic increment equation
(110)f0+δxf0′+12δx2f0′′=0
or, for that matter, any such higher-order algebraic equation, by the power series
(111)δx=(P+Qf0+Rf02+Sf03+Tf04+⋯)f0
and have upon substitution and collection
(112)(1+Pf0′)+(Qf0′+12P2f0′′)f0+(Rf0′+PQf0′′)f02+(Sf0′+12Q2f0′′+PRf0′′)f03+⋯=0,from which we deduce, by annulling lower order terms, that
(113)P=-1f0′,Q=-12P2v0,R=-PQv0,S=-(12Q2+PR)v0,T=-(QR+PS)v0,v0=f0′′f0′
and so on.The methods(114)x1=x0+Pf0+Qf02,x2=x0+Pf0+Qf02+Rf03
or
(115)x1=x0-(1+12w0)u0,x2=x0-(1+12w0+12w02)u0,u=ff′,v=f′′f′,w=uv
are both cubic
(116)x1-a=163f′′2-f′f′′′f′2(x0-a)3+O((x0-a)4),x2-a=-13f′′′f′(x-a0)3+O((x0-a)4)provided that f′(a)≠0.As heren=2, we take next
(117)ξ=23x0+13x1,
recalculate f′′(ξ), and verify that the second method in (114) is elevated thereby to quartic
(118)x1-a=172(45f′′3f′3+f(4)f′)(x0-a)4+O((x0-a)5).


## 14. Contrarily Converging Superquadratic Methods

We write(119)x2=x0-(1+Pr)u0,u0=f0f0′,r=f1f0,f1=f(x1),x1=x0-u0
for undetermined coefficient P and have
(120)x2-a=12f′′f′(1-P)(x0-a)2+O((x0-a)3).
We request that
(121)f′′f′(1-P)=2k(f′′f′)2,
for parameter k, or, in view of (51), that
(122)P=1-4kru0,
by which the iterative method in (119) becomes
(123)x2=x0-(1+r)u0+4kr2,u0=f0f0′,x1=x0-u0,f1=f(x1),r=f1f0,x2-a=k(f′′f′)2(x0-a)2+O((x0-a)3).
This superquadratic method converges from above if k>0, and from below if k<0.The interest in the method(124)x2=x0-21-ru0,x1=x0-2u0,f1=f(x1),u0=f0f0′,r=f1f0
is that it ultimately converges oppositely to Newton’s method
(125)x2-a=-12f′′f′(x0-a)2+O((x0-a)3)
as seen by comparing (125) with (17).For example, forf(x)=x+x2 and starting with x0=1/2, we compute from Newton’s method and from method (124), respectively,
(126)x1={12,18,180,16,560,143,046,720},x2={12,-110,-182,-16,562,-143,046,722},12(x1+x2)=12{11,140,13280,121,523,360,1926,510,094,425,920}
and root a=0 is bounded or bracketed thereby as
(127)-143,046,722<a<143,046,720.The average of Newton’s method and the method of (124) is cubic
(128)12(x1+x2)-a=16f′′′f′(x0-a)3+O((x0-a)4).The modified Halley method(129)x1=x0+det⁡[-1f0k2f0′]det⁡[f0′f0f0′′2f0′]f0=x0-2f0′+kf02f0′2-f0f0′′f0
is also superquadratic and one sided
(130)x1-a=-12k(x0-a)2+1123f′′2-2f′f′′′+3kf′f′′f′2×(x0-a)3+O((x0-a)4).According to error equation (130), if k<0, then convergence is at least asymptotically from above; if k>0, convergence is from below.


## 15. Alternating Superlinear and Supercubic Methods

We start by modifying Newton’s method as(131)x1=x0-(1+k)f0f0′
to have the superlinear method
(132)x1-a=-k(x0-a)+O((x0-a)2)
that ultimately converges alternatingly if 0<k<1.For example, forf(x)=x+x2, k=1/8, and x0=1, we compute by method (131) the alternating sequence
(133)x1={1,2.5·10-1,1.6·10-2,-1.7·10-3,2.1·10-4,K-2.7·10-5}
allowing us to bracket root a=0 as
(134)-2.6706·10-5<a<2.1406·10-4.For a higher-order method, we start with the originally quartic(135)x2=x0-(1+r+Qr2)u0,u0=f0f0′,r=f1f0,f1=f(x0-u0)
of (66) and have that
(136)x2-a=-k(f′′f′)2(x0-a)3+O((x0-a)4),KKKKKKKKKKKKKKk=14(Q-2).
This supercubic method converges alternatingly if parameter k>0.For example, forf(x)=x+x2, we generate by methods (135)-(136), with k=1, the alternating sequence
(137)x2={1,-0.012,8.34·10-6,-2.32·10-15}
and root a=0 is bracketed as
(138)-2.32·10-15<a<8.34·10-6.


## 16. Still Higher-Order Taylor Methods

Starting with(139)f(x)=f(x0+δx)=f0+δxf0′+12δx2f0′′+16δx3f′′′(ξ),kkkkkkkkkkkkkkkkkkkkkkkkkkkkx0<ξ<x0+δx,
we obtain the iterative method
(140)x1=x0+δx,f0+δxf0′+12δx2f0′′+16δx3f0′′′=0,
where
(141)δx=(P+Qf0+Rf02+Sf03+⋯)f0
with
(142)P=-1f′,Q=12P3f′′,R=P2(Qf′′+16P2f′′′),S=P(12Q2f′′+PRf′′+12P2Qf′′′).The methods(143)x1=x0+(P+Qf0+Rf02)f0,x2=x0+(P+Qf0+Rf02+Sf03)f0
are both quartic
(144)x2-a=124f(4)f′(x0-a)4+O((x0-a)5)
provided that f′(a)≠0.Recalculatingf′′′(ξ) at
(145)ξ=34x0+14x1
elevates the method to quintic
(146)x2-a=1960f′4(840f′′4-840f′f′′f′′′+80f′f′′′KKKKKK-3ff′3f(5))(x0-a)5+O((x0-a)6).


## 17. Repeated Fourth-Order Method

The repeated Newton method(147)x1=x0-f0f0′,x2=x1-f1f1′
is also quartic
(148)x2-a=18(f′′f′)3(x0-a)4+O((x0-a)5).Similarly, the repeated modified Newton method(149)x1=x0-mf0f0′,x2=x1-mf1f1′
remains quartic even near a root of any multiplicity m≥1(150)x2-a=1m3(BA)3(x0-a)4+O((x0-a)5).The repeated-step method(151)x1=x0-f0′f0′2-f0f0′′f0,x2=x1-f1′f1′2-f1f1′′f1
not requiring prior knowledge of the multiplicity index m of the root is also quartic
(152)x2-a=-1m3(BA)3(x0-a)4+O((x0-a)3).The two single-step methods(153)x1=x0-mf0f0′,x2=x0-f0′f0′2-f0f0′′f0
converge contrarily
(154)x1-a=1mBA(x0-a)2+O((x0-a)3),x2-a=-1mBA(x0-a)2+O((x0-a)3).
Their average is a cubic method
(155)x3-a=B2(m-1)-2ACm2A2m2(x0-a)3+O((x0-a)4),KKKKKKKKKKKKKKKKx3=12(x1+x2).For instance, forf(x)=x2+x3+x4,m=2, we compute by the two methods in (153) the sequences
(156)x1={0.5,1.1·10-1,6.3·10-3,2.0·10-5},x2={0.5,-1.3·10-1,-5.4·10-3,-1.4·10-5}x3=12(x1+x2)={0.5,-1.9·10-2,9.2·10-4,5.6·10-6}.


## 18. Stacked Higher-Order Methods and Simple Root

Higher-order, single-step methods can be written as a built-up power series ofu=f/f′(157)x1=x0-K2u02ndorderx1=x0-K2u0-K3u023rd⁡orderx1=x0-K2u0-K3u02-K4u034thorderx1=x0-K2u0-K3u02-K4u03-K5u045thorderx1=x0-K2u0-K3u02-K4u03-K5u04-K6u5kkkkkkkkkkkkkkkkkkkkkkkkk6thorder,
where
(158)K2=1,K3=12f′′f′,K4=16(3f′′2f′2-f′′′f′),K5=124(15f′′3f′3-10f′′f′′′f′2+f(4)f′),K6=1120(105f′′4f′4-105f′′2f′′′f′3+10f′′′2f′2KKKKKk+15f′′f(4)f′2-f(5)f′),
and so on. All evaluations in (158) are at x=x0.Alternatively, the method may be written as a product, such as the built-up sextic method(159)x1=x0-(1+P3u0)(1+P4u02)(1+P5u03)(1+P6u4)u0,
where
(160)P3=K3,P4=K4,P5=124(9f′′′3f′3-8f′′f′′′f′2+f(4)f′),P6=1240(165f′′4f′4-170f′′2f′′′f′3+20f′′′2f′2kkkkkkkk+25f′′f(4)f′2-2f(5)f′),
with all evaluations in (160) done at x=x0.


## 19. Higher-Order Methods and Multiple Root

The method(161)x1=x0-Pu0-Qu02-Ru03,u0=f0f0′,P=1,Q=12f0′′f0′,R=12(k+1)(f0′′f0′)2-16f0′′′f0′
is supercubic
(162)x1-a=-12k(f′′f′)2(x0-a)3+O((x0-a)4)
and converges alternatingly to a single root if k>0.Similarly, the method(163)x1=x0-Pu0-Qf0′′f0′u02-(R(f0′′f0′)2+Sf0′′′f0′)u03,kkkkkkkkkkkkkkkkkkkkkkkkkkkku0=f0f0′
with
(164)P=16m(11-6m+m2+3k(1-m)2),Q=12m2(2-m+2k(1-m)),R=12m3(1+k),S=-16m3
converges supercubically, and alternatingly if k>0(165)x1-a=-k2m2(BA)2(x0-a)3+O((x0-a)4)
to a root of any multiplicity m≥1.Osada [20] suggests the concise cubic method, compared with method (107),
(166)x1=x0-12m(m+1)f0f0′+12(m-1)2f0′f0′′,m>1x1-a=B2(m+1)2-2ACm(m-1)2A2m2(m-1)(x0-a)3+O((x0-a)4),
where A, B, and C are as in (80). See also Chun et al. [21].The method(167)x1=x0+det⁡[-1f0002f0′f003f0′′3f0′]det⁡[f0′f00f0′′2f0′f0f0′′′3f0′′3f0′]f0=x0-6f0′2-3f0f0′′6f0′3-6f0f0′f0′′+f02f0′′′f0.
is quartic if f′(a)≠0.The method(168)x1=x0-Pf0′2+Qf0f0′′Rf0′3+Sf0f0′f0′′+f02f0′′′f0,P=3+3m,kkkkkkkQ=-3,R=2+3m+1m2,S=-Pis quartic as well(169)x1-a=(g′3(1+3m+2m2)+gg′g′′(-3m-3m2)kkk+g2g′′′m2)(6g3m2)-1(x0-a)4+O((x0-a)5)
even in the event that root a is of multiplicity m≥1. In (169), g=g(a) is as in (46). See also [7].Householder’s [22] concise representation of the single-point iterative method of order p+2(170)x1=x0+(p+1)(f-1)0(p)(f-1)0(p+1)
unfolds, for p=2, into the quartic method of (167).An inexact value for the root multiplicity index plunges method (168) to mere first order. Indeed, replacing m in (168) by m(1+k) results in
(171)x1-a=-k3+2k(x0-a)+O((x0-a)2).To have a cubic method that does not require prior knowledge of the root multiplicity, we replacef by u to have
(172)x1=x0-v0-12u0′′u0′v02,v=uu′,u=ff′
for which
(173)x1-a=(1-m)B2+2mACA2m2(x0-a)3+O((x0-a)4).Still higher-order methods are readily thus generated.


## 20. Numerical Differentiation

All derivatives may, of course, be approximated by sufficiently accurate finite differences. Replacingf′ and f′′ in Halley’s method by the central finite-difference approximations
(174)f(x+h)-f(x-h)2h=f′(x)+16h2f′′′(x)+O(h4),f(x-h)-2f(x)+f(x+h)h2=f′′(x)+112h2f(4)(x)+O(h4)
we have, after some simplification,
(175)x1-a=h2(16f′′′f′(x0-a)+18f(4)f′(x0-a)2)+1123f′′-2f′f′′′f′2(x0-a)3+O((x0-a)4).Following Steffensen’s suggestion (see [23, 24]), we take h=f(x0) and have
(176)x1-a=1123f′′2-2f′f′′′+2f′3f′′′f′2(x0-a)3+O((x0-a)4).


## 21. Pell’s Equation

We consider now the very special case of the iterative solution off(x)=x2-N for accurate rational approximations p/q to N. A pair of natural numbers p,q that satisfies the general integer Pell equation
(177)p2-Nq2=k
is such that
(178)(pq)2=N+kq2orpq=N(1+k2Nq2)nearly, if k/(Nq2)≪1.Forx0=p0/q0 and x1=p1/q1, Newton’s method x1=x0-f0/f0′ assumes the split form
(179)p1=p02+Nq02,q1=2p0q0
and we verify that
(180)p12-Nq12=(p02-Nq02)2=k2ifp02-Nq02=k.Similarly, for Halley’s method(181)p1=(3Nq02+p02)p0,q1=Nq03+3p02q0
we verify that
(182)p12-Nq12=(p02-Nq02)3=k3ifp02-Nq02=k.The iterative method(183)p1=2Np0q0,q1=p02+Nq02
written out from (124) converges quadratically from below to N. Here
(184)p12-Nq12=-Nk2ifp02-Nq02=k.Still higher-order methods may be directly obtained from the expansion of(185)p1+Nq1=(p0+Nq0)n,
where n=2 is taken for Newton’s method, n=3 is taken for Halley’s method, and so on. See also [25, 26].For example, forn=4 and p0=3, q0=2,N=2, and p02-Nq02=1, we obtain from (185)
(186)p1+Nq1=(3+N2)4=577+N408,
and p1 = 577, q1 = 408, p12-Nq12= 1, and (p1/q1)2=2.000006.


## 22. Ratio and Mediant Corrections

In this simplest of all iterative methods, we start withp/q and add 1 to the numerator p if the ratio is an underestimate, or add 1 to the denominator q if the ratio is an overestimate. For example, seeking accurate rational approximations to the root of f(x)=x3-2, we start with p/q=4/4, carry out 625 such corrections, and secure the bounds
(187)286227<23<349277,349277-286227=162,879.The mediant fractionM of the upper and lower bounds in (187) is
(188)M=286+349227+277=635504.
We verify that M3<2 and have the closer bounds
(189)635504<23<349277,349277-635504=1139,608.We recall this well-known approximation theorem.Theorem 2.
Letk=m/n be a positive rational number in lowest terms. For any number ϵ>0 there is only a finite number of rational numbers, p/q≠m/n, such that
(190)|k-pq|<1q1+ϵ,ϵ>0.Proof.
Sincem/n-p/q≠0, the absolute value inequality may be recast as
(191)1≤|mq-np|<nqϵ,nqϵ≥1,ϵ>0.
Only a limited number of integers q satisfy the inequality n>qϵ. Consequently, there is only a limited number of integers p that satisfy the absolute value inequality of (191). End of proof.Unlike the classical method of bisections that may continue indefinitely, the method of mediants that generates rational approximations as in (190) terminates in a finite number of steps if the number it seeks to trap is rational. Conversely, if the method of mediants does not terminate in a finite number of steps, then the number being trapped is surely irrational.



---

*Source: 108976-2014-07-03.xml*