Application of a Generalized Secant Method to Nonlinear Equations with Complex Roots

Abstract

The secant method is a very effective numerical procedure used for solving nonlinear equations of the form $$f\left(x\right)=0$$. In a recent work (A. Sidi, Generalization of the secant method for nonlinear equations. Appl. Math. E-Notes, 8:115–123, 2008), we presented a generalization of the secant method that uses only one evaluation of $$f\left(x\right)$$ per iteration, and we provided a local convergence theory for it that concerns real roots. For each integer k, this method generates a sequence $$\left\{x_n\right\}$$ of approximations to a real root of $$f\left(x\right)$$, where, for $$n\ge k$$, $$x_{n+1}=x_n-f\left(x_n\right)/p_{n,k}^{\prime }\left(x_n\right)$$, $$p_{n,k}\left(x\right)$$ being the polynomial of degree k that interpolates $$f\left(x\right)$$ at $$x_n,x_{n-1},\dots ,x_{n-k}$$, the order $$s_k$$ of this method satisfying $$1<s_k<2$$. Clearly, when $$k=1$$, this method reduces to the secant method with $$s_1=(1+\sqrt{5})/2$$. In addition, $$s_1<s_2<s_3<\cdots ,$$ such that $${lim}_{k\to \infty }{s}_k=2$$. In this note, we study the application of this method to simple complex roots of a function $$f\left(z\right)$$. We show that the local convergence theory developed for real roots can be extended almost as is to complex roots, provided suitable assumptions and justifications are made. We illustrate the theory with two numerical examples.

1. Introduction

Let $\alpha$ be the solution to the equation $f\left(x\right)=0.$ An effective iterative method used for solving this equation that makes direct use of $f\left(x\right)$ (but no derivatives of $f\left(x\right)$) is the secant method that is discussed in many books on numerical analysis. See, for example, Atkinson [1], Dahlquist and Björck [2], Henrici [3], Ralston and Rabinowitz [4], and Stoer and Bulirsch [5]. See also the recent note [6] by the author, in which the treatment of the secant method and those of the Newton–Raphson, regula falsi, and Steffensen methods are presented in a unified manner.

Recently, this method was generalized by the author in [7] as follows: Starting with $x_0,x_1,\dots ,x_k$, $k+1$ initial approximations to $\alpha$, we generate a sequence of approximations $\left\{x_n\right\}$, via the recursion

(1)$x_{n+1}=x_n-\frac{f\left(x_n\right)}{p_{n,k}^{\prime }\left(x_n\right)},n=k,k+1,\dots ,$

$p_{n,k}^{\prime }\left(x\right)$ being the derivative of the polynomial $p_{n,k}\left(x\right)$ that interpolates $f\left(x\right)$ at the points $x_n,x_{n-1},\dots ,x_{n-k}$. (Thus, $p_{n,k}\left(x\right)$ is of degree k.) Clearly, the case $k=1$ is simply the secant method. In [7], we also showed that, provided $x_0,x_1,\dots ,x_k$ are sufficiently close to $\alpha$, the method converges with order $s_k$, that is, ${\displaystyle \underset{n\to \infty }{lim}\frac{|x_{n+1}-\alpha |}{|x_n{-\alpha |}^{s_k}}=C\ne 0}$ for some constant C, and that $1<s_k<2$. (We call $s_k$ the order of convergence of the method or the order of the method for short.) Here $s_k$ is the only positive root of the polynomial $s^{k+1}-{\sum }_{i=0}^k{s}^i$. We also have that

$\frac{1+\sqrt{5}}{2}=s_1<s_2<s_3<\cdots <2;\underset{k\to \infty }{lim}{s}_k=2.$

Actually, rounded to four significant figures,

$s_1\dot{=}1.618,s_2\dot{=}1.839,s_3\dot{=}1.928,s_4\dot{=}1.966,s_5\dot{=}1.984,s_6\dot{=}1.992,s_7\dot{=}1.996,\mathrm{etc}.$

Note that to compute $x_{n+1}$ we need knowledge of only $f\left(x_n\right),f\left(x_{n-1}\right),\dots ,f\left(x_{n-k}\right)$, and because $f\left(x_{n-1}\right),\dots ,f\left(x_{n-k}\right)$ have already been computed, $f\left(x_n\right)$ is the only new quantity to be computed. Thus, each step of the method requires $f\left(x\right)$ to be computed only once. From this, it follows that the efficiency index of this method is simply $s_k$ and that this index approaches 2 by increasing k even moderately, as can be concluded from the values of $s_1,\dots ,s_7$ given above.

In this work, we consider the application of this method to simple complex roots of a function $f\left(z\right)$, where z is the complex variable. Let us denote a real or complex root of $f\left(z\right)$ by $\alpha$ again; that is, $f\left(\alpha \right)=0$ and $f^{\prime }\left(\alpha \right)\ne 0$. Thus, starting with $z_0,z_1,\dots ,z_k$, $k+1$ initial approximations to $\alpha$, we generate a sequence of approximations $\left\{z_n\right\}$ via the recursion

(2)$z_{n+1}=z_n-\frac{f\left(z_n\right)}{p_{n,k}^{\prime }\left(z_n\right)},n=k,k+1,\dots ,$

$p_{n,k}^{\prime }\left(z\right)$ being the derivative of the polynomial $p_{n,k}\left(z\right)$ that interpolates $f\left(z\right)$ at the points $z_n,z_{n-1},\dots ,z_{n-k}$. As in [7], we can use Newton’s interpolation formula to generate $p_{n,k}\left(z\right)$ and $p_{n,k}^{\prime }\left(z\right)$. Thus

(3)$p_{n,k}\left(z\right)=f\left(z_n\right)+\sum _{i=1}^{k}f[z_n,z_{n-1},\dots ,z_{n-i}]\prod _{j=0}^{i-1}(z-z_{n-j})$

and

(4)$p_{n,k}^{\prime }\left(z_n\right)=f[z_n,z_{n-1}]+\sum _{i=2}^{k}f[z_n,z_{n-1},\dots ,z_{n-i}]\prod _{j=1}^{i-1}(z_n-z_{n-j}).$

Here, $g[{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m]$ is the divided difference of order m of the function $g\left(z\right)$ over the set of points $\{{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m\}$ and is a symmetric function of these points. For details, we refer the reader to [7].

As proposed in [7], we generate the $k+1$ initial approximations as follows: We choose the approximations $z_0,z_1$ first. We then generate $z_2$ by applying our method with $k=1$ (that is, with the secant method). Next, we apply our method to $z_0,z_1,z_2$ with $k=2$ and obtain $z_3$, and so on, until we have generated all $k+1$ initial approximations, via

(5)$z_{n+1}=z_n-\frac{f\left(z_n\right)}{p_{n,n}^{\prime }\left(z_n\right)},n=1,2,\dots ,k-1.$

Remark 1.

• 1. 

  Instead of choosing $z_1$ arbitrarily, we can generate it as $z_1=z_0+f\left(z_0\right)$ as suggested in Brin [8], which is quite sensible since $f\left(z\right)$ is small near the root α. We can also use the method of Steffensen—which uses only $f\left(z\right)$ and no derivatives of $f\left(z\right)$—to generate $z_1$ from $z_0$; thus,

  $z_1=z_0-\frac{{\left[f\left(z_0\right)\right]}^2}{f(z_0+f\left(z_0\right))-f\left(z_0\right)}.$

• 2. 

  It is clear that, in case $f\left(z\right)$ takes on only real values along the $Rez$ axis and we are looking for nonreal roots of $f\left(z\right)$, at least one of the initial approximations must be chosen to be nonreal.

• 3. 

  We would like to mention that Kogan, Sapir, and Sapir [9] have proposed another generalization of the secant method for simple real roots of nonlinear equations $f\left(x\right)=0$ that resembles our method described in (1). In the notation of (1), this method produces a sequence of approximations $\left\{x_n\right\}$ via

  (6) $x_{n+1}=x_n-\frac{f\left(x_n\right)}{p_{n,n}^{\prime }\left(x_n\right)},n=1,2,\dots ,$

  starting with arbitrary $x_0$ and $x_1$, and it is of order 2. Note that, in (6), $p_{n,n}\left(x\right)$ interpolates $f\left(x\right)$ at the points $x_0,x_1,\dots ,x_n,$ hence is of degree n, which is tending to infinity. In (1), $p_{n,k}\left(x\right)$ is of degree k, which is fixed.

• 4. 

  Yet another generalization of the secant method for finding simple real roots of $f\left(x\right)$ was recently given by Nijmeijer [10]. This method too requires no derivative information, requires one evaluation of $f\left(x\right)$ per iteration, and has the same order of convergence as our method. It follows an idea of applying a convergence acceleration method, such as Aitken’s ${\Delta }^2$-process, to approximations obtained from the secant method, as proposed by Han and Potra [11]. Because Nijmeijer’s method is not based on polynomial interpolation, it is completely different from our method, however. For Aitken’s ${\Delta }^2$-process, see [1,2,3,4,5]. See also [12] (Chapter 15) by the author.

In the next section, we analyze the local convergence properties of the method as it is applied to complex roots. We show that the analysis of [7] can be extended to the complex case following some clever manipulation. We prove that the order $s_k$ of the method is the same as that we discovered in the real case. In Section 3, we provide two numerical examples to confirm the results of our convergence analysis.

2. Local Convergence Analysis

We now turn to the analysis of the sequence ${\left\{z_n\right\}}_{n=0}^{\infty }$ that is generated via (2). Our treatment covers all $k\ge 1$.

In our analysis, we will make use of the Hermite-Genocchi formula that provides an integral representation for divided differences (For a proof of this formula, see Atkinson [1], for example). Even though this formula is usually stated for functions defined on real intervals, it is easy to verify (see Filipsson [13], for example) that it also applies to functions defined in the complex plane under proper assumptions. Thus, provided $g\left(z\right)$ is analytic on E, a bounded closed convex set in the complex plane, and provided ${\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m$ are in E, there holds

(7)$g[{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m]=\underset{S_m}{\int \cdots \int }{g}^{\left(m\right)}(t_0{\zeta }_0+t_1{\zeta }_1+\cdots +t_m{\zeta }_m)d{t}_1\cdots d{t}_m,t_0=1-\sum _{i=1}^m{t}_i.$

Here $S_m$ is the m-dimensional simplex defined as

(8)$S_m=\left\{(t_1,\dots ,t_m)\in {\mathbb{R}}^m:t_i\ge 0,i=1,\dots ,m,\sum _{i=1}^m{t}_i\le 1\right\}.$

We note that (7) holds whether the ${\zeta }_i$ are distinct or not. We also note that $g[{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m]$ is a symmetric and continuous function of its arguments.

By the conditions we have imposed on $g\left(z\right)$, it is easy to see that the integrand $g^{\left(m\right)}\left({\sum }_{i=0}^m{t}_i{\zeta }_i\right)$ in (7) is always defined because ${\sum }_{i=0}^m{t}_i{\zeta }_i$ is in the set E and $g\left(z\right)$ is analytic on E. This is so because, by (7) and (8),

$(t_1,\dots ,t_m)\in S_m\Rightarrow t_i\ge 0,i=0,1,\dots ,m,\mathrm{and}\sum _{i=0}^m{t}_i=1,$

which implies that ${\sum }_{i=0}^m{t}_i{\zeta }_i$ is a convex combination of ${\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m$ hence is in the set $C=\mathrm{conv}\{{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m\}$, the convex hull of the points ${\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m$, and $C\subseteq E$. Consequently, taking moduli on both sides of (7), we obtain, for all ${\zeta }_i$ in E,

(9)$\begin{array}{cc}\hfill |g[{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m]|& \le \underset{S_m}{\int \cdots \int }|g^{\left(m\right)}\left(\sum _{i=0}^m{t}_i{\zeta }_i\right)|d{t}_1\cdots d{t}_m\hfill \\ \hfill & \le \frac{\parallel g^{\left(m\right)}\parallel }{m!},\parallel g^{\left(m\right)}\parallel =\underset{z\in E}{max}|g^{\left(m\right)}\left(z\right)|.\hfill \end{array}$

In addition, since ${\sum }_{i=0}^m{t}_i=1$ in (7), as ${\zeta }_i\to \widehat{\zeta }$ for all $i=0,1,\dots ,m$, there hold ${\sum }_{i=0}^m{t}_i{\zeta }_i\to \widehat{\zeta }$ and $g^{\left(m\right)}\left({\sum }_{i=0}^m{t}_i{\zeta }_i\right)\to g^{\left(m\right)}\left(\widehat{\zeta }\right)$, and hence

(10)$\underset{\begin{array}{c}{\zeta }_i\to \widehat{\zeta }\\ i=0,1,\dots ,m\end{array}}{lim}g[{\zeta }_0,{\zeta }_1,\dots ,{\zeta }_m]=g\underset{m+1\mathrm{times}}{\underbrace{[\widehat{\zeta },\widehat{\zeta },\dots ,\widehat{\zeta }]}}=\frac{g^{\left(m\right)}\left(\widehat{\zeta }\right)}{m!}.$

In (9) and (10), we have also invoked the fact that (see [14] (p. 346), for example)

$\underset{S_m}{\int \cdots \int }d{t}_1\cdots d{t}_m=\frac{1}{m!}.$

We will make use of these in the proof of our main theorem that follows. This theorem and its proof are almost identical to that given in [7] once we take into account, where and when needed, the fact that we are now working in the complex plane. For convenience, we provide all the details of the proof.

Theorem 1.

Let α be a simple root of $f\left(z\right)$, that is, $f\left(\alpha \right)=0$, but $f^{\prime }\left(\alpha \right)\ne 0$. Let $B_r$ be the closed disk of radius r containing α as its center, that is,

(11) $B_r=\{z\in \mathbb{C}:|z-\alpha |\le r\}.$

Let $f\left(z\right)$ be analytic on $B_r$. Choose a positive integer k and let $z_0,z_1,\dots ,z_k$ be distinct initial approximations to α. Generate $z_{k+1},z_{k+2},\dots$ via

(12) $z_{n+1}=z_n-\frac{f\left(z_n\right)}{p_{n,k}^{\prime }\left(z_n\right)},n=k,k+1,\dots ,$

where $p_{n,k}\left(z\right)$ is the polynomial of interpolation to $f\left(z\right)$ at the points $z_n,z_{n-1},\dots ,z_{n-k}$. Then, provided $z_0,z_1,\dots ,z_k$ are in $B_r$ and sufficiently close to α, we have the following cases:

• 1. 

  If $f^{(k+1)}\left(\alpha \right)\ne 0$, the sequence $\left\{z_n\right\}$ converges to α, and

  (13) $\underset{n\to \infty }{lim}\frac{{ϵ}_{n+1}}{{\prod }_{i=0}^k{ϵ}_{n-i}}=\frac{{(-1)}^{k+1}}{(k+1)!}\frac{f^{(k+1)}\left(\alpha \right)}{f^{\prime }\left(\alpha \right)}\equiv L;{ϵ}_n=z_n-\alpha \forall n.$

  The order of convergence is $s_k$, $1<s_k<2$, where $s_k$ is the only positive root of the polynomial $g_k\left(s\right)=s^{k+1}-{\sum }_{i=0}^k{s}^i$ and satisfies

  (14) $2-2^{-k-1}e<s_k<2-2^{-k-1}\mathit{for}k\ge 2;s_k<s_{k+1};\underset{k\to \infty }{lim}{s}_k=2,$

  e being the base of natural logarithms, and

  (15) $\underset{n\to \infty }{lim}\frac{|{ϵ}_{n+1}|}{|{ϵ}_n{|}^{s_k}}={|L|}^{(s_k-1)/k},$

  which also implies that

  (16) $s_k=\underset{n\to \infty }{lim}\frac{log|{ϵ}_{n+1}/{ϵ}_n|}{log|{ϵ}_n/{ϵ}_{n-1}|}.$

• 2. 

  If $f\left(z\right)$ is a polynomial of degree at most k, the sequence $\left\{z_n\right\}$ converges to α, and

  (17)$\underset{n\to \infty }{lim}\frac{{ϵ}_{n+1}}{{ϵ}_n^2}=\frac{f^{″}\left(\alpha \right)}{2f^{\prime }\left(\alpha \right)};{ϵ}_n=z_n-\alpha \forall n.$

  Thus, $\left\{z_n\right\}$ converges with order 2 if $f^{″}\left(\alpha \right)\ne 0$, and with order greater than 2 if $f^{″}\left(\alpha \right)=0$.

Proof. 

We start by deriving a closed-form expression for the error in $z_{n+1}$. Subtracting $\alpha$ from both sides of (12), and noting that

$f\left(z_n\right)=f\left(z_n\right)-f\left(\alpha \right)=f[z_n,\alpha ](z_n-\alpha ),$

we have

(18)$z_{n+1}-\alpha =\left(1-\frac{f[z_n,\alpha ]}{p_{n,k}^{\prime }\left(z_n\right)}\right)(z_n-\alpha )=\frac{p_{n,k}^{\prime }\left(z_n\right)-f[z_n,\alpha ]}{p_{n,k}^{\prime }\left(z_n\right)}(z_n-\alpha ).$

We now note that

(19)$p_{n,k}^{\prime }\left(z_n\right)-f[z_n,\alpha ]=\left\{p_{n,k}^{\prime }\left(z_n\right)-f^{\prime }\left(z_n\right)\right\}+\left\{f^{\prime }\left(z_n\right)-f[z_n,\alpha ]\right\},$

and that

(20)$f^{\prime }\left(z_n\right)-p_{n,k}^{\prime }\left(z_n\right)=f[z_n,z_n,z_{n-1},\dots ,z_{n-k}]\prod _{i=1}^k(z_n-z_{n-i})$

and

(21)$f^{\prime }\left(z_n\right)-f[z_n,\alpha ]=f[z_n,z_n]-f[z_n,\alpha ]=f[z_n,z_n,\alpha ](z_n-\alpha ).$

Note that (20) can be obtained by starting with the divided difference representation of $f\left(z\right)-p_{n,k}\left(z\right)$, namely, $f\left(z\right)-p_{n,k}\left(z\right)=f[z,z_n,z_{n-1},\dots ,z_{n-k}]{\prod }_{i=0}^k(z-z_{n-i}),$ and by computing ${lim}_{z\to z_n}[f\left(z\right)-p_{n,k}\left(z\right)]/{\prod }_{i=0}^k(z-z_{n-i})$ via L’Hôpital’s rule.

For simplicity of notation, let

(22)$-f[z_n,z_n,z_{n-1},\dots ,z_{n-k}]={\widehat{D}}_n\mathrm{and}f[z_n,z_n,\alpha ]={\widehat{E}}_n,$

and rewrite (19) and (20) as

(23)$\begin{array}{cc}& p_{n,k}^{\prime }\left(z_n\right)-f[z_n,\alpha ]={\widehat{D}}_n\prod _{i=1}^k({ϵ}_n-{ϵ}_{n-i})+{\widehat{E}}_n{ϵ}_n,\hfill \end{array}$

(24)$\begin{array}{cc}\hfill & p_{n,k}^{\prime }\left(z_n\right)=f^{\prime }\left(z_n\right)+{\widehat{D}}_n\prod _{i=1}^k({ϵ}_n-{ϵ}_{n-i}).\hfill \end{array}$

Substituting these into (18), we finally obtain

(25)${ϵ}_{n+1}=C_n{ϵ}_n;C_n\equiv \frac{p_{n,k}^{\prime }\left(z_n\right)-f[z_n,\alpha ]}{p_{n,k}^{\prime }\left(z_n\right)}=\frac{{\widehat{D}}_n{\prod }_{i=1}^k({ϵ}_n-{ϵ}_{n-i})+{\widehat{E}}_n{ϵ}_n}{f^{\prime }\left(z_n\right)+{\widehat{D}}_n{\prod }_{i=1}^k({ϵ}_n-{ϵ}_{n-i})}.$

We now prove that convergence takes place. First, let us assume without loss of generality that $f^{\prime }\left(z\right)\ne 0$ for all $z\in B_r$, and set $m_1={min}_{z\in B_r}|f^{\prime }\left(z\right)|>0$. (This is possible since $\alpha \in B_r$ and $f^{\prime }\left(\alpha \right)\ne 0,$ and we can choose r as small as we wish to also guarantee $m_1>0$.) Next, let $M_s={max}_{z\in B_r}|f^{\left(s\right)}\left(z\right)|/s!,$ $s=1,2,\dots .$ Thus, assuming that $\{z_n,z_{n-1},\dots ,z_{n-k}\}\subset B_r$ and noting that $B_r$ is a convex set, we have by (9) that

$|{\widehat{D}}_n|\le M_{k+1},|{\widehat{E}}_n|\le M_2,\mathrm{because}\{\alpha ,z_n,z_{n-1},\dots ,z_{n-k}\}\subset B_r.$

Next, choose the ball $B_{t/2}$ sufficiently small (with $t/2\le r$) to ensure that $m_1>2M_{k+1}{t}^k+M_{2}t/2$. It can now be verified that, provided $z_n,z_{n-1},\dots ,z_{n-k}$ are all in $B_{t/2}$, there holds

$\begin{array}{cc}\hfill |C_n|& \le \frac{M_{k+1}{\prod }_{i=1}^k|{ϵ}_n-{ϵ}_{n-i}|+M_2|{ϵ}_n|}{m_1-M_{k+1}{\prod }_{i=1}^k|{ϵ}_n-{ϵ}_{n-i}|}\hfill \\ \hfill & \le \frac{M_{k+1}{\prod }_{i=1}^k(|{ϵ}_n|+|{ϵ}_{n-i})|+M_2|{ϵ}_n|}{m_1-M_{k+1}{\prod }_{i=1}^k(|{ϵ}_n|+|{ϵ}_{n-i}|)}\le \overline{C},\hfill \end{array}$

where

$\overline{C}\equiv \frac{M_{k+1}{t}^k+M_{2}t/2}{m_1-M_{k+1}{t}^k}<1.$

Consequently, by (25), $|{ϵ}_{n+1}|\le \overline{C}|{ϵ}_n|<|{ϵ}_n|$, which implies that $z_{n+1}\in B_{t/2}$, just like $z_n,z_{n-1},\dots ,z_{n-k}$. Therefore, if $z_0,z_1,\dots ,z_k$ are chosen in $B_{t/2}$, then $|C_n|\le \overline{C}<1$ for all $n\ge k$, hence $\left\{z_n\right\}\subset B_{t/2}$ and ${lim}_{n\to \infty }{z}_n=\alpha$.

As for (13) when $f^{(k+1)}\left(\alpha \right)\ne 0$, we proceed as follows: By the fact that ${lim}_{n\to \infty }{z}_n=\alpha$, we first note that, by (20) and (21),

(26)$\underset{n\to \infty }{lim}{p}_{n,k}^{\prime }\left(z_n\right)=f^{\prime }\left(\alpha \right)=\underset{n\to \infty }{lim}f[z_n,\alpha ],$

and thus ${lim}_{n\to \infty }{C}_n=0$. This means that ${lim}_{n\to \infty }({ϵ}_{n+1}/{ϵ}_n)=0$ and, equivalently, that $\left\{z_n\right\}$ converges with order greater than 1. As a result,

$\underset{n\to \infty }{lim}({ϵ}_n/{ϵ}_{n-i})=0\mathrm{for}\mathrm{all}i\ge 1,$

and

${ϵ}_n/{ϵ}_{n-i}=o({ϵ}_n/{ϵ}_{n-j})\mathrm{as}n\to \infty ,\mathrm{for}j<i.$

Consequently, expanding in (25) the product ${\prod }_{i=1}^k({ϵ}_n-{ϵ}_{n-i})$, we have

(27)$\begin{array}{cc}\hfill \prod _{i=1}^k({ϵ}_n-{ϵ}_{n-i})& =\prod _{i=1}^k\left(-{ϵ}_{n-i}[1-{ϵ}_n/{ϵ}_{n-i}]\right)\hfill \\ \hfill & ={(-1)}^k\left(\prod _{i=1}^k{ϵ}_{n-i}\right)[1+O({ϵ}_n/{ϵ}_{n-1})]\mathrm{as}n\to \infty .\hfill \end{array}$

Substituting (27) into (25), and defining

(28)$D_n=\frac{{\widehat{D}}_n}{p_{n,k}^{\prime }\left(z_n\right)},E_n=\frac{{\widehat{E}}_n}{p_{n,k}^{\prime }\left(z_n\right)},$

we obtain

(29)${ϵ}_{n+1}={(-1)}^k{D}_n\left(\prod _{i=0}^k{ϵ}_{n-i}\right)[1+O({ϵ}_n/{ϵ}_{n-1})]+E_n{ϵ}_n^2\mathrm{as}n\to \infty .$

Dividing both sides of (29) by ${\prod }_{i=0}^k{ϵ}_{n-i}$, and defining

(30)${\sigma }_n=\frac{{ϵ}_{n+1}}{{\prod }_{i=0}^k{ϵ}_{n-i}},$

we have

(31)${\sigma }_n={(-1)}^k{D}_n[1+O({ϵ}_n/{ϵ}_{n-1})]+E_n{\sigma }_{n-1}{ϵ}_{n-k-1}\mathrm{as}n\to \infty .$

Now, by (10), (22), and (26),

(32)$\underset{n\to \infty }{lim}{D}_n=-\frac{1}{(k+1)!}\frac{f^{(k+1)}\left(\alpha \right)}{f^{\prime }\left(\alpha \right)},\underset{n\to \infty }{lim}{E}_n=\frac{f^{\left(2\right)}\left(\alpha \right)}{2f^{\prime }\left(\alpha \right)}.$

Because ${lim}_{n\to \infty }{D}_n$ and ${lim}_{n\to \infty }{E}_n$ are finite, and because ${lim}_{n\to \infty }({ϵ}_n/{ϵ}_{n-1})=0$ and ${lim}_{n\to \infty }{ϵ}_{n-k-1}=0$, it follows that there exist a positive integer N and positive constants $\beta <1$ and D, with $|E_n{ϵ}_{n-k-1}|\le \beta$ when $n>N$, for which (31) gives

(33)$|{\sigma }_n|\le D+\beta |{\sigma }_{n-1}|\mathrm{for}\mathrm{all}n>N.$

Using (33), it is easy to show that

$|{\sigma }_{N+s}|\le D\frac{1-{\beta }^s}{1-\beta }+{\beta }^s|{\sigma }_N|,s=1,2,\dots ,$

which, by the fact that $\beta <1$, implies that $\left\{{\sigma }_n\right\}$ is a bounded sequence. Making use of this fact, we have ${lim}_{n\to \infty }{E}_n{\sigma }_{n-1}{ϵ}_{n-k-1}=0$. Substituting this into (31), and invoking (32), we next obtain ${lim}_{n\to \infty }{\sigma }_n={(-1)}^k{lim}_{n\to \infty }{D}_n=L$, which is precisely (13).

That $s_k$, the order of the method, as defined in the statement of the theorem, satisfies (14) and (15) follows from Traub [15] (Chapter 3). We provide a simplified treatment of this topic in Appendix A.

This completes the proof of part 1 of the theorem.

When $f\left(z\right)$ is a polynomial of degree at most k, we first observe that $f^{(k+1)}\left(z\right)=0$ for all z, which implies that $p_{n,k}\left(z\right)=f\left(z\right)$ for all z, hence also $p_{n,k}^{\prime }\left(z\right)=f^{\prime }\left(z\right)$ for all z. Therefore, we have that $p_{n,k}^{\prime }\left(z_n\right)=f^{\prime }\left(z_n\right)$ in the recursion of (12). Consequently, (12) becomes

$z_{n+1}=z_n-\frac{f\left(z_n\right)}{f^{\prime }\left(z_n\right)},n=k,k+1,\dots ,$

which is the recursion for the Newton–Raphson method. Thus, (17) follows. This completes the proof of part 2 of the theorem. □

3. Numerical Examples

In this section, we present two numerical examples that we treated with our method. Our computations were done in quadruple-precision arithmetic (approximately 35-decimal-digit accuracy). Note that in order to verify the theoretical results concerning iterative methods with order greater than unity, we need to use computer arithmetic of high precision (preferably, of variable precision, if available) because the number of correct significant decimal digits in the $z_n$ increases dramatically from one iteration to the next as we are approaching the solution.

In both examples below, we take $k=2$. We choose $z_0$ and $z_1$ and compute $z_2$ using one step of the secant method, namely,

(34)$z_2=z_1-\frac{f\left(z_1\right)}{f[z_0,z_1]}.$

Following that, we compute $z_3,z_4,\dots ,$ via

(35)$z_{n+1}=z_n-\frac{f\left(z_n\right)}{f[z_n,z_{n-1}]+f[z_n,z_{n-1},z_{n-2}](z_n-z_{n-1})},n=2,3,\dots .$

In our examples, we have carried out our computations for several sets of $z_0,z_1$, and we have observed essentially the same behavior that we observe in Table 1 and Table 2.

Example 1.

Consider $f\left(z\right)=0$, where $f\left(z\right)=z^3-8$, whose solutions are ${\alpha }_r=2e^{\mathrm{i}2\pi r/3}$, $r=0,1,2$. We would like to obtain the root ${\alpha }_1=2e^{\mathrm{i}2\pi /3}=-1+\mathrm{i}\sqrt{3}$. We chose $z_0=2\mathrm{i}$ and $z_1=-2+2\mathrm{i}$. The results of our computations are given in Table 1.

From (13) and (16) in Theorem 1, we should have

$\underset{n\to \infty }{lim}\frac{{ϵ}_{n+1}}{{ϵ}_n{ϵ}_{n-1}{ϵ}_{n-2}}=\frac{{(-1)}^3}{3!}\frac{f^{‴}\left({\alpha }_1\right)}{f^{\prime }\left({\alpha }_1\right)}=\frac{1}{24}(1-\mathrm{i}\sqrt{3})=0.04166\cdots -\mathrm{i}0.07216\cdots$

and

$\underset{n\to \infty }{lim}\frac{log|{ϵ}_{n+1}/{ϵ}_n|}{log|{ϵ}_n/{ϵ}_{n-1}|}=s_2=1.83928\cdots ,$

and these seem to be confirmed in Table 1. Furthermore, in infinite-precision arithmetic, $z_9$ should have close to 60 correct significant figures; we do not see this in Table 1 due to the fact that the arithmetic we have used to generate Table 1 can provide an accuracy of at most 35 digits.

Example 2.

Consider $f\left(z\right)=0$, where $f\left(z\right)=sin\left(\mathrm{i}z\right)-cosz$. $f\left(z\right)$ has infinitely many roots ${\alpha }_r=(1-\mathrm{i})(\pi /4+r\pi )$, $r=0,\pm 1,\pm 2,\dots$. We would like to obtain the root ${\alpha }_0=(1-\mathrm{i})\pi /4$. We chose $z_0=1.5-\mathrm{i}1.3$ and $z_1=0.6-\mathrm{i}0.5$. The results of our computations are given in Table 2.

From (13) and (16) in Theorem 1, we should have

$\underset{n\to \infty }{lim}\frac{{ϵ}_{n+1}}{{ϵ}_n{ϵ}_{n-1}{ϵ}_{n-2}}=\frac{{(-1)}^3}{3!}\frac{f^{‴}\left({\alpha }_1\right)}{f^{\prime }\left({\alpha }_1\right)}=-\frac{\mathrm{i}}{6}=-\mathrm{i}0.16666\cdots$

and

$\underset{n\to \infty }{lim}\frac{log|{ϵ}_{n+1}/{ϵ}_n|}{log|{ϵ}_n/{ϵ}_{n-1}|}=s_2=1.83928\cdots ,$

and these seem to be confirmed in Table 2. Furthermore, in infinite-precision arithmetic, $z_8$ should have close to 50 correct significant figures; we do not see this in Table 2 due to the fact that the arithmetic we have used to generate Table 2 can provide an accuracy of at most 35 digits.

Remark 2.

In relation to the examples we have just presented, we would like to discuss the issue of estimating the relative errors $|{ϵ}_n/\alpha |$ in the $z_n$. This should help the reader when studying the numerical results included in Table 1 and Table 2. Starting with (13) and (15), we first note that, for all large n,

$|{ϵ}_{n+1}{|\approx |L|}^{(s_k-1)/k}{|{ϵ}_n|}^{s_k}.$

Therefore, assuming also that $\alpha \ne 0$, we have

$|{ϵ}_{n+1}/\alpha |\approx D|{ϵ}_n{/\alpha |}^{s_k},D={\left({|L|}^{1/k}|\alpha |\right)}^{s_k-1}.$

Now, if $z_n$ has $q>0$ correct significant figures, we have $|{ϵ}_n/\alpha |=O\left({10}^{-q}\right)$. If, in addition, $D=O\left({10}^r\right)$ for some r, then we will have

$|{ϵ}_{n+1}/\alpha |\approx O\left({10}^{r-q{s}_k}\right).$

For simplicity, let us consider the case $r=0$, which is practically what we have in the two examples we have treated. Then $z_{n+1}$ has approximately $q{s}_k$ correct significant decimal digits. That is, if $z_n$ has q correct significant decimal digits, then, due to the fact that $s_k>1$, $z_{n+1}$ will have $s_k$ times as many correct significant decimal digits as $z_n$.

Funding

This research received no external funding.

Acknowledgments

The author would like to thank Tamara Kogan for drawing his attention to the paper [9] mentioned in the Introduction.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A

Before ending, we would like to provide a brief treatment of the order of convergence of our method stated in (14) and (15) by considering

$\frac{{ϵ}_{n+1}}{{\prod }_{i=0}^k{ϵ}_{n-i}}=L\forall n\iff {ϵ}_{n+1}=L\prod _{i=0}^k{ϵ}_{n-i}\forall n,$

instead of (13). We will show that $|{ϵ}_{n+1}|=Q|{ϵ}_n{|}^{s_k}$ is possible if $s_k$ is a solution to the polynomial equation $s^{k+1}={\sum }_{i=0}^k{s}^i$ and $Q={|L|}^{(s_k-1)/k}$. (For a more detailed treatment, we refer the reader to [15] (Section 3.3)).

We start by expressing all $|{ϵ}_{n-i}|$ in terms of $|{ϵ}_n|$. We have

$|{ϵ}_{n-i}|=\frac{|{ϵ}_n{|}^{1/s_k^i}}{Q^{m_i}},m_i=\sum _{j=1}^i\frac{1}{s_k^j},i=1,2,\dots .$

Substituting this into $|{ϵ}_{n+1}|=|L|{\prod }_{i=0}^k|{ϵ}_{n-i}|$, we obtain

${Q|{ϵ}_n|}^{s_k}=|L||{ϵ}_n|\prod _{i=1}^k\frac{|{ϵ}_n{|}^{1/s_k^i}}{Q^{m_i}}=\frac{|L|}{Q^M}{|{ϵ}_n|}^{\rho };\rho =\sum _{i=0}^k\frac{1}{s_k^i},M=\sum _{i=1}^k{m}_i.$

Of course, this is possible when $s_k=\rho$ and $Q^{M+1}=|L|$.

Now, the requirement that $s_k=\rho$ is the same as $s_k^{k+1}={\sum }_{i=0}^k{s}_k^i$, which implies that the order $s_k$ should be a root of the polynomial

$g_k\left(s\right)=s^{k+1}-\sum _{i=0}^k{s}^i=\frac{s^{k+2}-2s^{k+1}+1}{s-1}.$

By Descartes’ rule of signs, $g_k\left(s\right)$ has only one positive root, which we denote by $\tilde{s}$. Since $g_k\left(1\right)=-k<0$ and $g_k\left(2\right)=1>0$, we have that $1<\tilde{s}<2$. The remaining k roots of $g_k\left(s\right)$ are the zeroes of the polynomial $\tilde{g}\left(s\right)=g_k\left(s\right)/(s-\tilde{s})={\sum }_{j=0}^k{c}_j{s}^j$, the $c_j$ satisfying $\tilde{s}{c}_0=1$ and $\tilde{s}{c}_j-c_{j-1}=1$, $j=1,\dots ,k,$ hence

$c_j=\frac{1}{\tilde{s}}\sum _{i=0}^j\frac{1}{{\tilde{s}}^i},j=0,1,\dots ,k\Rightarrow 0<c_0<c_1\cdots <c_{k-1}<c_k=1.$

Therefore, by the Eneström–Kakeya theorem, all k roots of $\tilde{g}\left(s\right)$ are in the unit disk. We thus conclude that $\tilde{s}=s_k$ since we already know that the order of our method is greater than 1. (For Descartes’ rule of signs and the Eneström–Kakeya theorem, see, for example, Henrici [16] (pp. 442, 462)).

Next, we note that $g_k\left(s\right)=s{g}_{k-1}\left(s\right)-1$. Therefore, $g_{k-1}\left(s_{k-1}\right)=0$ implies $g_k\left(s_{k-1}\right)=-1<0$, which, along with $g_k\left(2\right)=1>0$, implies that $s_{k-1}<s_k<2$. Therefore, the sequence ${\left\{s_k\right\}}_{k=1}^{\infty }$ is monotonically increasing and is bounded from above by 2. Consequently, ${lim}_{k\to \infty }{s}_k=\widehat{s}$ exists and $\widehat{s}\le 2$. Now,

$g_k\left(s_k\right)=0\Rightarrow s_k^{k+2}-2s_k^{k+1}+1=0\Rightarrow s_k^2-2s_k=-\frac{1}{s_k^k}.$

Upon letting $k\to \infty$ on both sides, we obtain ${\widehat{s}}^2-2\widehat{s}=0$, which gives $\widehat{s}=2$.

The expression given for M can be simplified considerably as we show next. First, it is easy to verify that

$M=\sum _{i=1}^k\frac{k-i+1}{s_k^i}=\frac{1}{s_k^k}\sum _{i=1}^{k}i{s}_k^{i-1}.$

Next,

$s_k^{k}M=\left(\frac{d}{ds}\sum _{i=0}^k{s}^i\right){|}_{s=s_k}=\left(\frac{d}{ds}\frac{s^{k+1}-1}{s-1}\right){|}_{s=s_k}=\frac{(k+1)s^k(s-1)-(s^{k+1}-1)}{{(s-1)}^2}{|}_{s=s_k}.$

By $s^{k+1}-1=(s-1){\sum }_{i=0}^k{s}^i$, this becomes

$s_k^{k}M=\frac{(k+1)s^k-{\sum }_{i=0}^k{s}^i}{s-1}{|}_{s=s_k}=\frac{k{s}_k^k-{\sum }_{i=0}^{k-1}{s}_k^i}{s_k-1}.$

Now, by the fact that $g_k\left(s_k\right)=0$, we have ${\sum }_{i=0}^{k-1}{s}_k^i=s_k^{k+1}-s_k^k$. Consequently,

$M=\frac{k-(s_k-1)}{s_k-1}\Rightarrow M+1=\frac{k}{s_k-1},$

which is the required result.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Tables

Table 1

Results obtained by applying the generalized secant method with $k=2$, as shown in (34) and (35), to the equation $z^3-8=0$, to compute the root ${\alpha }_1=-1+\mathrm{i}\sqrt{3}$. The entries denoted “**” mean that the limit of the extended-precision arithmetic has been reached.

n	|${\mathit{ϵ}}_{\mathit{n}}$|	${\displaystyle \frac{{\mathit{ϵ}}_{\mathit{n}+1}}{{\mathit{ϵ}}_{\mathit{n}}{\mathit{ϵ}}_{\mathit{n}-1}{\mathit{ϵ}}_{\mathit{n}-2}}}$	${\displaystyle \frac{log|{\mathit{ϵ}}_{\mathit{n}+1}/{\mathit{ϵ}}_{\mathit{n}}|}{log|{\mathit{ϵ}}_{\mathit{n}}/{\mathit{ϵ}}_{\mathit{n}-1}|}}$
0	$1.035D+00$	-	-
1	$1.035D+00$	-	-
2	$4.808D-01$	$-8.972D-02$ $+\mathrm{i}1.015D-01$	2.516
3	$6.979D-02$	$1.224D-01$ $-\mathrm{i}2.727D-02$	1.437
4	$4.355D-03$	$1.009D-01$ $-\mathrm{i}4.079D-02$	2.023
5	$1.591D-05$	$4.561D-02$ $-\mathrm{i}9.794D-02$	1.839
6	$5.223D-10$	$3.793D-02$ $-\mathrm{i}7.268D-02$	1.839
7	$2.967D-18$	$3.741D-02$ $-\mathrm{i}7.579D-02$	1.838
8	$2.083D-33$	**	**
9	$0.000D+00$	**	**

Table 2

Results obtained by applying the generalized secant method with $k=2$, as shown in (34) and (35), to the equation $sin\left(\mathrm{i}z\right)-cosz=0$, to compute the root ${\alpha }_0=(1-\mathrm{i})\pi /4$. The entries denoted “**” mean that the limit of the extended-precision arithmetic has been reached.

n	|${\mathit{ϵ}}_{\mathit{n}}$|	${\displaystyle \frac{{\mathit{ϵ}}_{\mathit{n}+1}}{{\mathit{ϵ}}_{\mathit{n}}{\mathit{ϵ}}_{\mathit{n}-1}{\mathit{ϵ}}_{\mathit{n}-2}}}$	${\displaystyle \frac{log|{\mathit{ϵ}}_{\mathit{n}+1}/{\mathit{ϵ}}_{\mathit{n}}|}{log|{\mathit{ϵ}}_{\mathit{n}}/{\mathit{ϵ}}_{\mathit{n}-1}|}}$
0	$6.608D-01$	-	-
1	$3.403D-01$	-	-
2	$1.341D-01$	$3.163D-01$ $+\mathrm{i}1.397D-01$	2.743
3	$1.043D-02$	$1.466D-01$ $-\mathrm{i}1.846D-01$	1.774
4	$1.122D-04$	$-2.943D-03$ $-\mathrm{i}1.117D-01$	1.934
5	$1.755D-08$	$9.223D-03$ $-\mathrm{i}1.614D-01$	1.766
6	$3.320D-15$	$-7.686D-04$ $-\mathrm{i}1.658D-01$	1.857
7	$1.084D-27$	**	**
8	$9.630D-35$	**	**