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

The secant method is a very effective numerical procedure used for solving nonlinear equations of the form $f(x)=0$. In a recent work [A. Sidi, Generalization of the secant method for nonlinear equations. {\em Appl. Math. E-Notes}, 8:115--123, 2008] we presented a generalization of the secant method that uses only one evaluation of $f(x)$ per iteration, and we provided a local convergence theory for it that concerns real roots. For each integer $k$, this method generates a sequence $\{x_n\}$ of approximations to a real root of $f(x)$, where, for $n\geq k$, $x_{n+1}=x_n-f(x_n)/p'_{n,k}(x_n)$, $p_{n,k}(x)$ being the polynomial of degree $k$ that interpolates $f(x)$ at $x_n,x_{n-1},\ldots,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 and $\lim_{k\to\infty}s_k=2$. In this note, we study the application of this method to simple complex roots of a real or complex function $f(z)$. 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.


Introduction
Let α be the solution to the equation f (x) = 0. An effective iterative method used for solving (1.1) that makes direct use of f (x) [but no derivatives of f (x)] is the secant method that is discussed in many books on numerical analysis. See, for example, Atkinson [1], Henrici [5], Ralston and Rabinowitz [6], and Stoer and Bulirsch [9]. See also the recent note [7] 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 [8] as follows: Starting with x 0 , x 1 , . . . , x k , k + 1 initial approximations to α, we generate a sequence of approximations {x n }, via the recursion , n = k, k + 1, . . . , p ′ n,k (x) being the derivative of the polynomial p n,k (x) that interpolates f (x) at the points x n , x n−1 , . . . , x n−k . Clearly, the case k = 1 is simply the secant method. In [8], we also showed that provided x 0 , x 1 , . . . , x k are sufficiently close to α, the method converges of order s k , where s k ∈ (1, 2) and is the only positive root of the polynomial s k+1 − k i=0 s i . We also have that 1 + √ 5 2 = s 1 < s 2 < s 3 < · · · < 2; lim k→∞ s k = 2.
Actually, rounded to four significant figures, Note that to compute x n+1 we need knowledge of only f (x n ), f (x n−1 ), . . . , f (x n−k ), and because f (x n−1 ), . . . , f (x n−k ) have already been computed, f (x n ) is the only new quantity to be computed. Thus, each step of the method requires f (x) 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.
In this work, we consider the application of this method to simple complex roots of a function f (z), where z is the complex variable. Let us denote a real or complex root of f (z) by α again; that is, f (α) = 0 and f ′ (α) = 0. Thus, starting with z 0 , z 1 , . . . , z k , k + 1 initial approximations to α, we generate a sequence of approximations {z n } via the recursion , n = k, k + 1, . . . , (1.1) p ′ n,k (z) being the derivative of the polynomial p n,k (z) that interpolates f (z) at the points z n , z n−1 , . . . , z n−k . As in [8], we can use Newton's interpolation formula to generate p n,k (z) and p ′ n,k (z). Thus Here g[ζ 0 , ζ 1 , . . . , ζ m ] is the divided difference of order m of the function g(z) over the set of points {ζ 0 , ζ 1 , . . . , ζ m } and is a symmetric function of these points. For details, we refer the reader to [8].
As proposed in [8], 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 initial approximations, via Remarks.
1. Instead of choosing z 1 arbitrarily, we can generate it as z 1 = z 0 + f (z 0 ) as suggested in Brin [2], which is quite sensible since f (z) is small near the root α. We can also use the method of Steffensen-which uses only f (z) and no derivatives of f (z)-to generate z 1 from z 0 ; thus, 2. It is clear that, in case f (z) takes on only real values along the Re z axis and we are looking for nonreal roots of f (z), at least one of the initial approximations must be chosen to be nonreal.
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 [8] 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.

Local convergence analysis
We now turn to the analysis of the sequence {z n } ∞ n=0 that is generated via (1.1). We treat the case k ≥ 2. The case k = 1 is similar but much simpler.
In our analysis, we will make use of the Hermite-Genocchi formula that provides an integral representation for divided differences. 1 Even though this formula is usually stated for functions defined on real intervals, it is easy to verify (see Filipsson [4], for example) that it also applies to functions defined in the complex plane under proper assumptions. Thus, provided g ∈ C m (E), E being a bounded closed convex set in the complex plane, and provided ζ 0 , ζ 1 , . . . , ζ m are in E, there holds Here S m is the m-dimensional simplex defined as and We also note that (2.1) holds whether the ζ i are distinct or not. By the conditions we have imposed on g(z) it is easy to see that the integrand which implies that m i=0 t i ζ i is a convex combination of ζ 0 , ζ 1 , . . . , ζ m hence is in the set C = conv{ζ 0 , ζ 1 , . . . , ζ m }, the convex hull of the points ζ 0 , ζ 1 , . . . , ζ m , and C ⊂ E. Consequently, taking moduli on both sides of (2.1), we obtain, for all ζ i in E, In (2.4) and (2.5), we have also invoked the fact that (see [3, p. 346], for example) 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 [8], but taking 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 2.1 Let α be a simple root of f (z), that is, f (α) = 0, but f ′ (α) = 0. Assume that f ∈ C k+1 (B r ), where B r is a closed ball of radius r containing α as its center, that is, Let z 0 , z 1 , . . . , z k be distinct initial approximations to α, and generate z n , n ≥ k + 1, via , n = k, k + 1, . . . , where p n,k (z) is the polynomial of interpolation to f (z) at the points z n , z n−1 , . . . , z n−k . Then, provided z 0 , z 1 , . . . , z k are in B r and sufficiently close to α, we have the following cases: 1. If f (k+1) (α) = 0, the sequence {z n } converges to α, and The order of convergence is s k , 1 < s k < 2, where s k is the only positive root of the equation s k+1 = k i=0 s i and satisfies e being the base of natural logarithms, and lim n→∞ |ǫ n+1 | |ǫ n | s k = |L| (s k −1)/k ⇒ s k = lim n→∞ log |ǫ n+1 /ǫ n | log |ǫ n /ǫ n−1 | . (2.10) 2. If f (z) is a polynomial of degree at most k, the sequence {z n } converges to α, and lim n→∞ ǫ n+1 ǫ 2 n = f ′′ (α) 2f ′ (α) ; ǫ n = z n − α ∀ n. Proof. We start by deriving a closed-form expression for the error in z n+1 . Subtracting α from both sides of (2.7), and noting that we have We now note that and that and and rewrite (2.13) and (2.14) as Substituting these in (2.12), we finally obtain . (
From (2.8) and (2.10) in Theorem 2.1, we should have and lim n→∞ log |ǫ n+1 /ǫ n | log |ǫ n /ǫ n−1 | = s 2 = 1.83928 · · · , and these seem to be confirmed in Table 1. Also, in infinite-precision arithmetic, x 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.
and lim n→∞ log |ǫ n+1 /ǫ n | log |ǫ n /ǫ n−1 | = s 2 = 1.83928 · · · , and these seem to be confirmed in Table 2. Also, in infinite-precision arithmetic, x 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. Before closing, we would like to discuss the issue of estimating the relative errors |ǫ n /α| in the z n . This should help the reader when studying the numerical results included in Tables 1 and 2. Starting with (2.8) and (2.10), we first note that, for all large n, |ǫ n+1 | ≈ |L| (s k −1)/k |ǫ n | s k .
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 qs 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 .