On Homogeneous Linear Recurrence Relations and Approximation of Zeros of Complex Polynomials Bahman Kalantari Abstract. Let p(z) be a complex polynomial of degree n. To each complex number a we associate a sequence called the Basic Sequence {Bm(a) = a − p(a)Dm−2 (a)/Dm−1(a)}, where Dm (a) is defined via a homogeneous linear recurrence relation and depends only on the normalized derivatives p(i) (a)/i!. Each Dm (a) is also representable as a Toeplitz determinant. Except possibly for the locus of points equidistant to two distinct roots, given any input a, the Basic Sequence converges to a root of p. The roots of p partition the Euclidean plane into Voronoi regions. Under some regularity assumption (e.g. simplicity of the roots), for almost all inputs within the Voronoi polygon of a root, the corresponding Basic Sequence converges to that root. The discovery of the Basic Sequence, its error estimates, and several of its properties are consequences of our previous analysis of a fundamental family of iteration functions {Bm(z)}, called the Basic Family. Given any fixed m ≥ 2 and an appropriate input a0 , the fixed-point iterates ak+1 = Bm (ak ) converge to a simple root having order m. The Basic Sequence corresponds to the pointwise evaluation of the Basic Family. We present several fractal images that confirm the theoretical convergence results: as m increases, the basins of attractions to the roots, as computed with respect to the iteration function Bm (z), rapidly converge to the Voronoi regions. Thus the regions with chaotic behavior rapidly shrink to the boundaries of the Voronoi regions. The interplay between the Basic Sequence and the Basic Family result in new algorithms for approximation of polynomial roots. Also as a byproduct of the analysis we give determinantal representation of solution of homogeneous linear recurrence relations, easily computable bounds on the modulus of terms of such recurrence relations, and fast computation of determinant of certain Toeplitz matrices. Finally, we prove a relationship between the Basic Sequence and the iterates of Bernoulli Method for approximation of extreme roots.
1. Introduction. The problem of approximating roots of a real or complex polynomial p(z) is a fundamental and classic problem. A comprehensive bibliography is given by McNamee [18]. Some aspects of the history, applications, and algorithms, are described in Pan [20]. For a complexity result see Smale [22]. New results on this problem continue to be discovered. In [6], Kalantari et al. study 1991 Mathematics Subject Classification. Primary 12D10, 65Y20, 41A20, 11B39; Secondary 28A80, 30C10. Key words and phrases. Zeros of Complex Polynomials, Basic Family of Iteration Functions, Recurrence Relations, Convergence. 1
2
BAHMAN KALANTARI
a fundamental family of iteration function, called the Basic Family, {Bm (z) = z − p(z)Dm−2 (z)/Dm−1 (z)}∞ m=2 , where Dm (z) is a Toeplitz determinant defined with respect to the normalized derivatives of p(z), but is also derivable via a homogeneous linear recurrence relation. Many results on the properties of the members of this family including their close tie with a determinantal generalization of Taylor’s theorem can be found in [5], [6], [7], [9], [10] [8], [11], and [13]. The first member of the sequence, B2 (z), is Newton’s iteration function, and B3 (z) is Halley’s iteration function. For the history of these two iteration functions see Ypma [25] and Traub [23]. In this paper we will prove several new properties of the Basic Family. On the one hand these properties are derived by studying general homogeneous linear recurrence relations that are defined via a special set of initial conditions. On the other hand the properties of the Basic Family reveal new results on such recurrence relations. In tracing a circle of relations, we are led to the study of a classical method for approximation of dominant roots of polynomial, called the Bernoulli Method. As a byproduct we even reveal new properties on this method. Thus, not only in this paper we give new results on approximation of polynomial roots, but on homogeneous linear recurrence relations as well. Let p(z) be a given complex polynomial of degree n. To each point a in the complex plane we associate a sequence {dm (a)} defined via a homogeneous linear recurrence relation, derived from a special representation of the Taylor’s expansion of p(z) at a. The Basic Sequence at a is defined to be {bm (a) = a−p(a)dm−2 (a)/dm−1 (a)}. We will show that dm (a) = Dm (a). Thus the basic sequence corresponds to the pointwise evaluation of the Basic Family. The Basic Sequence has important properties with respect to the approximation of roots of p(z). One of these properties is as follows: Except possibly for the locus of points that are equidistant from two distinct roots (a set of measure zero consisting of the union of a finite set of lines) for any input a, the corresponding Basic Sequence is well-defined and it converges to a root of p(z). Under some regularity assumptions, e.g. the simplicity of all the roots, the Basic Sequence at a converges to the root of p(z) that is nearest to a. By virtue of this, regions of convergence in the complex plane can best be described by Voronoi regions. The roots partition the plane into Voronoi regions. Voronoi region of a root is the set of points that are closer to this root than to any other root. The Voronoi regions consist of convex polygons and their boundaries (see e.g. [21]). Thus, under some regularity assumptions, the Voronoi region of each root is the set of points whose corresponding Basic Sequence converges to that root. As an example of the Basic Sequence consider the complex polynomial p(z) = 1 + p1 z + · · · + pn z n , and let a = 0. The associated sequence becomes dm (0) =
n X (−1)i−1 pidm−i (0), i=1
where d0 (0) = 1, and dj (0) = 0 for j < 0. The evaluation of the corresponding Basic Sequence {bm (0)} in many cases provides a very simple scheme for the approximation of a root of p(z) (often the root nearest to the origin). This approach is very different than the classic iterative method of Newton which requires repeated polynomial evaluations.
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
3
In order to estimate the error in the approximation of a root θ via the individual member of the Basic Sequence, i.e. the quantity |bm (a)−θ|, as well as the discovery of further properties of the Basic Sequence, we will make use a nontrivial determinantal generalization of Taylor’s theorem derived in [9]. We will prove that in fact the Basic Sequence corresponds to the pointwise evaluation of the Basic Family, {Bm (z)} (i.e. Bm (a) = bm (a)). In particular it follows that the Basic Sequence has the following important property: given a fixed m ≥ 2, for any a0 within a certain radius of a simple root, the iterates defined by ak+1 = bm (ak ) converge to that root with order m. The discovery of the Basic Sequence is motivated by our previous analysis of the Basic Family. The results proved in this paper are consequences of the interplay between properties of the Basic Sequence proved here and those of the Basic Family. Not only this interplay gives new algorithms for approximation of polynomial roots, it results in determinantal representation of solution of homogeneous linear recurrence relations defined via a single initial input. This representation allows an easily computable upper bound on the modulus of the terms via Hadamard inequality. Also in some cases it is possible to state lower bounds using a recent determinantal lower bound, as well as the use of Gerschgorin’s theorem. This paper is organized as follows: In Section 2 we derive properties of general homogeneous recurrence relations defined via a single nonzero initial condition. In Section 3, we define the Basic Sequence and state several of its properties. In Section 4, we give the relationship between the Basic Sequence and the Basic Family of iteration functions, thereby state additional properties of the Basic Sequence. In Section 5, we give the connection between the Basic Sequence and the Bernoulli Method. In Section 6, we give some properties on the solution of general homogeneous linear recurrence relations that are defined via a single initial condition. In Section 7, we present some experimental results via a graphics package that produces fractal images. Section 8 consists of our concluding results. 2. Properties of Homogeneous Linear Recurrence Relations Having a Single Nonzero Initial Condition. Let √ numbers. For a complex number c = a + ib √ C be the field of complex (i = −1) its modulus is |c| = a2 + b2 . Consider a general homogeneous linear recurrence relation of degree n, i.e., (2.1)
am = c1 am−1 + c2 am−2 + · · · + cn am−n ,
where ci ’s are given numbers in C, and cn 6= 0. We assume that a0 6= 0 and for all j < 0, aj = 0. Without loss of generality we assume that the initial conditions are (2.2)
a0 = 1, , aj = 0, ∀ j < 0.
We shall refer to (2.1) and (2.2) as homogeneous linear recurrence relation with a single initial condition. The characteristic polynomial of the sequence {am } is the polynomial (2.3)
Q(z) = z n − (c1 z n−1 + c2 z n−2 + · · · + cn−1 z + cn ).
The characteristic equation refers to the equation (2.4)
Q(z) = 0.
4
BAHMAN KALANTARI
It is well known that am can be represented in terms of the roots of Q(z). More precisely, let η1 , . . . , ηt be the set of distinct roots of Q(z) having multiplicities n1 , . . . , nt . Then for each i = 1, . . . , t there exist a polynomial qi(z), such that either qi(z) is identically zero, or it is a polynomial of degree ≤ ni − 1, such that am = q1 (m)η1m + · · · + qt (m)ηtm , ∀ m ≥ −n + 1.
(2.5)
Proposition 2.1. Consider the homogeneous linear recurrence relation with a single initial initial condition defined by (2.1) and (2.2). Assume that |ηi| = 6 |ηj |, if i 6= j. Let r be the index satisfying |ηr | = max{|ηi| : qi(z) 6≡ 0, i = 1, . . . , t}
(i.e. the largest modulus root whose corresponding coefficient polynomial in (2.5) is not the zero polynomial). Then r is well-defined and am lim = ηr . m→∞ am−1 Moreover, if all the roots of Q(z) are simple, then qi(z) is a nonzero constant for all i = 1, . . . , n. In particular, ηr is the root of Q(z) with the largest modulus. Proof. Since a0 = 1 we conclude that not all qi(z)’s are zero. Hence r is welldefined. From the assumption that all distinct roots of Q(z) have different modulus and the definition of ηr , for any i 6= r such that qi (z) 6≡ 0, we have |ηi/ηr | < 1. Since for any such index qi(z) is a polynomial it follows that m ηi (2.6) = 0. lim qi(m) m→∞ ηr First we claim that the sequence {am } cannot have a subsequence {amk } where amk = 0 for all k. Suppose such a subsequence exists. If all but one of the qi(z)’s are identically zero, then from (2.1) we have am = ql (m)ηlm , for some l and ql (m) 6≡ 0. Thus, am 6= 0 for all m ≥ 0. So under the assumption that {amk = 0} there must exist at least two indices for which qi(z) is not identically zero. For any a, b ∈ C, |a − b| ≥ |a| − |b|, we have (2.7) t mk amk X ηi m = qi(mk ) ≥ qr (mk ) − ηr k ηr i=1
i=1, i6=r, qi (z)6≡0
Taking the limit in (2.7) and using (2.6), we get (2.8)
t X
mk qi(mk ) ηi . ηr
0 ≥ lim |qr (mk )|. k→∞
But the limit in (2.8) is either a nonzero number, or infinity. This is a contradiction. Thus am can be zero at most for a finite set of indices m. Since qi(z) is a polynomial, if it is not identically zero then we have (2.9)
lim
m→∞
qi(m) = 1. qi (m − 1)
Thus for m large enough the ratio am /am−1 is well-defined and is given by: m Pt ηrm i=1 qi (m) ηηri am (2.10) = m−1 am−1 ηi m−1 Pt ηr i=1 qi (m − 1) ηr
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
5
Now using (2.6) and (2.9) it follows that the ratio in (2.10) converges to ηr . Now assume that all the roots of Q(z) are simple. In this case t = n, and qi (z) is a constant, say qi ∈ C, for all i = 1, . . . , n. We wish to show that qi 6= 0 for all i = 1, . . . , n. From the initial condition (2.2), while substituted into (2.5), we get the following system of linear equations in qi’s −n+1 η1 η1−n+2 .. . −1 η 1 1
η2−n+1 η2−n+2 .. . η2−1 1
ηn−n+1 q1 0 q2 0 ηn−n+2 .. .. = .. . . . −1 q 0 ηn n−1 q 1 1 n
··· ··· ··· ··· ···
Suppose that one of the qi ’s is zero. Without loss of generality we may assume that q1 = 0. Substituting q1 = 0 in the above system of linear equations gives: −n+1 η2 η2−n+2 .. . η2−1
η3−n+1 η3−n+2 .. .
··· ··· ··· ···
η3−1
ηn−n+1 0 q2 q3 0 ηn−n+2 .. .. = .. . . . ηn−1
qn
0
But the coefficient matrix is a product of a Vandermonde matrix in η2 , . . . , ηn and the diagonal matrix diag(η1−2 , · · · , ηnn−1). Since the ηi’s are distinct the Vandermonde matrix is invertible. Hence if q1 = 0, then q2 = · · · = qn = 0. But this contradicts that a0 6= 0. We now review some properties of am which are important from the computational point of view. From (2.1) and the identity equations ai = ai , i = 1, . . . , m, we get
am c1 am−1 1 .. = .. . . am−n
0
c2 0 .. . 0
··· ··· ··· ···
cn−1 0 .. . 1
am−1 cn 0 am−2 .. .. . . am−n−1 0
Thus, if for each k ≥ 0 we let vk = (ak , · · · , ak−n )T , and let A be the coefficient matrix above, then we have (2.11)
vm = Avm−1 = · · · = Am v0 , v0 = (1, 0, . . . , 0)T .
The matrix A is the companion matrix of the characteristic polynomial Q(z). Using (2.11) it can be shown that vm and hence am can be computed in O(n2 log m) arithmetic operations, see Gries and Levin [1]. It is also possible to compute am in O(n log n log m) arithmetic operations, see Fiduccia [4]. 3. The Basic Sequence: Approximation of Polynomial Roots Using a Homogeneous Linear Recurrence Relation. Let (3.1)
p(z) = p0 + p1 z + · · · + pn z n ,
6
BAHMAN KALANTARI
be a given polynomial of degree n with coefficients pi ∈ C, i = 1, . . . , n. We wish to approximate roots of p(z). For each point a ∈ C we associate a sequence {dm (a)} defined via the following homogeneous linear recurrence relation of degree n, with a single nonzero initial condition: (3.2)
dm (a) =
n X p(i) (a) (−1)i−1 pi−1 (a) dm−i(a), i! i=1
where (3.3)
d0 (a) = 1, dj (a) = 0, ∀ j < 0.
Define the characteristic polynomial of a, denoted by Qa(z), to be the characteristic polynomial corresponding to the above recurrence relation, i.e. n X p(i) (a) n−i n (−1)i−1 pi−1 (a) Qa(z) = z − z . i! i=1
From the relationship between a recurrence relation and its characteristic polynomial it follows that if η1 , . . . , ηt is the set of distinct roots of Qa (z) having multiplicities n1 , . . . , nt , then there exist polynomials q1 (z), . . . , qt (z) such that for each i = 1, . . . , t, qi (z) is either identically zero or a polynomial of degree ≤ ni − 1, such that (3.4)
dm (a) = q1 (m)η1m + · · · + qt (m)ηtm , ∀ m ≥ −n + 1.
Theorem 3.1. If η is a root of Qa (z), then θ = a − p(a)/η is a root of p(z). Conversely, if θ is a root of p(z), then η = −p(a)/(θ − a) is a root of Qa (z). Proof. Suppose that η is a root of Qa (z), i.e. (3.5)
ηn −
n X p(i)(a) n−i (−1)i−1 pi−1 (a) η = 0. i! i=1
Multiplying both sides by p(a)/η n we get (3.6)
p(a) −
n X p(i)(a) n−i η = 0. (−1)i−1 pi (a) i! i=1
Let θ = a − p(a)/η. Then (3.7)
η=−
Substituting (3.7) in (3.6) gives (3.8)
n X p(i) (a) i=0
i!
p(a) . (θ − a) (θ − a)i = 0.
But from Taylor’s theorem the above coincides with p(θ). Hence θ is a root of p(z). The converse follows simply by reversing the steps. Let (3.9)
Rp = {θ1 , . . . , θt }
be the set of distinct roots of p(z). The elements of Rp partition the Euclidean plane into Voronoi regions and their boundaries. The Voronoi region of a root θ is
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
7
a convex polygon defined by the locus of points which are closer to this root than to any other root. More precisely, the Voronoi region of a root θ is (3.10)
V (θ) = {z ∈ C : |z − θ| < |z − θ′ |, θ′ ∈ Rp , θ′ 6= θ}.
Let Sp be the locus of points that are equidistant from two distinct roots, i.e. (3.11)
Sp = {z ∈ C : |z − θ| = |z − θ′ |, where θ, θ′ ∈ Rp , θ 6= θ′ }.
This is a set of measure zero consisting of the union of finite number of lines. Definition 3.2. Given a ∈ C the Basic Sequence at a is defined as bm (a) = a − p(a)
dm−2 (a) , m = 2, 3, . . .. dm−1 (a)
Theorem 3.3. Given p(z), for any input a 6∈ Sp , the Basic Sequence is welldefined satisfying lim bm (a) = θ, m→∞
for some θ ∈ Rp . Moreover, for each root θ of p(z), given any a ∈ V (θ), the dominant root of Qa (z) is η = −p(a)/(θ − a). Thus if η m appears with a nonzero coefficient polynomial in the representation of dm (a) (see (3.4)), then bm (a) converges to θ. In particular, if all roots of p(z) are simple, then for all a ∈ V (θ), limm→∞ bm (a) = θ. Proof. From Theorem 3.1 for any a ∈ C, different than a root of p(z), the set of distinct roots of Qa (z) is (3.12)
Tp (a) = {η1 = −
p(a) p(a) , . . . , ηt = − }. θ1 − a θt − a
If a 6∈ Sp , then |ηi| = 6 |ηj |, for i 6= j. Now from Proposition 2.1 it follows that lim
m→∞
dm (a) = η, dm−1 (a)
for some η ∈ Tp (a). From (3.12) for each a ∈ C the dominant root of Qa (z) is −p(a)/(θ − a), where θ is the closest root to a. From Proposition 2.1 it then follows that if η m appears with a nonzero coefficient in the representation of dm (a), then bm (a) converges to θ. If all roots of p(z) are simple, then so are the roots of Qa (z). Thus from Proposition 2.1 any a ∈ V (θ), bm (a) converges to θ. Remark 3.4. For points that belong to Sp the convergence of bm (a) to a root 2 may not occur. For instance if p(z) √ = z −√2, then the origin is a point that belongs to Sp since it is equidistant to 2 and − 2. It is easy to see that for this choice of a the corresponding sequence {dm (a)} has a subsequence that is zero. Thus, the Basic Sequence is undefined. 4. The Basic Sequence and its Connection to the Basic Family In this and the next sections we prove several important properties of the sequences {dm (a)} and {bm (a)}. In doing so we first need to relate them to some previously known results. We will next describe a family of iteration functions for polynomial root-finding and state some of its known properties. As before let p(z) be a polynomial of degree n ≥ 2 with complex coefficients. Set D0 (z) ≡ 1, and for each natural number m ≥ 1, define
8
BAHMAN KALANTARI
p′′ (z) 2!
p′ (z)
p(z) Dm (z) = det 0 . ..
p′ (z) p(z) .. .
0
0
... .. . .. . .. . ...
p(m−1)(z) (m−1)!
..
.
..
.
..
. p(z)
p(m) (z) (m)!
p(m−1)(z) (m−1)!
.. .
p′′ (z) 2! ′
p (z)
where det denotes the determinant. Also, given any m ≥ 1, for each i = m + 1, . . . , n + m − 1 define p′′ (z) p′′′ (z) p(m) (z) p(i) (z) . . . 3! (m)! i! 2! .. ′ p(i−1) (z) p(m−1) (z) p′′ (z) . p (z) 2! (m−1)! (i−1)! . . . b m,i (z) = det . . . ′ D . . . p(z) p (z) . .. .. p′′ (z) p(i−m+2) (z) .. . . 2! (i−m+2)! p(i−m+1) (z) 0 0 ... p′ (z) (i−m+1)! For each m ≥ 2, define Bm (z) ≡ z − p(z)
Dm−2 (z) . Dm−1 (z)
Theorem 4.1. (Kalantari [13]) The following conditions hold: 1. For all m ≥ 1 we have, Dm (z) =
n X pi−1 (z)p(i) (z) Dm−i (z), Dj = 0, j < 0. (−1)i−1 i! i=1
2. Let θ be a simple root of p(z). Then, Bm (z) = θ +
m+n−2 X
(−1)m
i=m
b m−1,i (z) D (z − θ)i . Dm−1 (z)
3. There exists r > 0 such that given any a0 ∈ Nr (θ) = {z : |z − θ| ≤ r}, the fixed-point iteration ak+1 = Bm (ak ) is well-defined, it converges to θ having order m. Specifically, b b (θ − ak+1 ) m−1 Dm−1,m (θ) m−1 Dm−1,m (θ) = (−1) = (−1) . k→∞ (θ − ak )m Dm−1 (θ) p′ (θ)m−1 lim
4. Let θ be a simple root of p(z). There exists a neighborhood of θ, N ∗ (θ), such that for each a within this neighborhood p′ (a) 6= 0, |a − θ| < 1, and |p′ (a)| − ∗
n X
i=0,i6=1
|p(a)|
i−1 2
|p(i)(a)| 1 ≥ |p′ (a)|. i! 2
For any a ∈ N (θ), if we set 2 1/2 X n (i) |p (a)| , h(a) = i! i=0
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
9
then we have
m 2h(a) 1 |Bm (a) − θ| ≤ ′ (a − θ) . p (a) 1 − |a − θ|
5. In particular, if θ is a simple root of p(z), there exists r∗ ∈ (0, 1) such that given any a ∈ Nr∗ (θ), we have θ = lim Bm (a). m→∞
Theorem 4.2. Given p(z), for any input a ∈ C the associated sequences {dm (a)} and the corresponding Basic Sequence {bm (a)} satisfy: dm (a) = Dm (a), bm (a) = Bm (a), ∀ m ≥ 0.
In particular, the following hold: 1. Dm (a) and hence Bm (a) can be computed in O(n log n log m) operations. 2. dm (a) can be represented as an m × m Toeplitz matrix. 3. Dm (a) (and dm (a)) is the first entry of the vector A(a)m v0 , where v0 = (1, 0, . . ., 0)T ∈ Rn and A(a) is the companion matrix of Qa (z), i.e.
(n)
(−1)n−1 pn−1 (a) p 1 .. . 0
(a) n!
(n−1)
(−1)n−2 pn−2 (a) p 0 .. .
(a) n−1!
0
··· ··· ··· ···
′′
−p(a) p 2!(a) 0 .. . 1
p′ (a) 0 .. . 0
4. For each simple root θ of p(z) there exists an open neighborhood of the root such given any fixed m ≥ 2, and any a0 within this neighborhood, the sequence {ak } defined by ak+1 = bm (ak ),
converges to θ having order of convergence equal to m. Proof. The equality dm (a) = Dm (a) is immediate since they correspond to the same homogeneous recurrence relations satisfying the same initial conditions. The complexity of computing Dm (a) and Bm (a) follows from two known results. Firstly, given any a the evaluation of p(j) (a)/j!, j = 0, . . . , n, called the normalized derivatives, can be established in O(n log2 n) arithmetic operations, see Kung [16]. Secondly, the m-th term of a sequence defined via a homogeneous linear recurrence relation can be computed in O(n log n log m) arithmetic operations, see Fiduccia [4]. Combing these two results the claimed complexity follows. The remaining parts of the theorem follow from Theorems 3.3 and 4.1. 5. The Basic Sequence and the Bernoulli Iteration Here we make a connection between the Basic Sequence and a classic method known as the Bernoulli Iteration or the Bernoulli Method due to Daniel Bernoulli, see Hildebrand [2], Householder [3]. Consider again a complex polynomial p(z) = pn z n + · · · + p1 z + p0 of degree n, with p0 6= 0, having t distinct roots Rp = {θ1 , . . . , θt }, of multiplicity n1 , . . . , nt , respectively. We can consider p(z) as the characteristic polynomial of the recurrence relation (5.1)
pn um + pn−1 um−1 + · · · + p0 um−n = 0.
10
BAHMAN KALANTARI
For any set of initial conditions we get (5.2)
um = β1 (m)θ1 + · · · + βt (m)θt ,
where for each i = 1, . . . , t, βi (z) is either identically zero or a polynomial of degree ni − 1. It thus follows, as in the proof of Proposition 2.1, that if the coefficient polynomial corresponding to the dominant root is nonzero, then the ratio um /um−1 converges to the dominant root of p(z). Different schemes for selecting the initial conditions so that under the assumption of simplicity of the roots convergence of ratio um /um−1 to the dominant root is guaranteed, are discussed in several numerical analysis books, see e.g. Hildebrand [2]. In order to construct sequences convergent to the least modulus root, one can consider the polynomial g(y) = p0 yn + p1 yn−1 + · · · + pn .
(5.3)
Clearly, θ is a root of p(z) if and only if 1/θ is a root of g(y). Thus if the sequence {vm } is defined via the homogeneous recurrence relation (5.4)
p0 vm + p1 vm−1 + · · · + pn vm−n = 0,
then, under some regularity conditions, the ratio vm−1 /vm converges to the least root of p(z). The approximation of the extreme roots via the ratios um /um−1 and vm−1 /vm is what is known as the Bernoulli Method, see Householder [3] and Hildebrand [2]. In what follows we will show that through the use of Taylor’s theorem the notion of least modulus and largest modulus is only relative to the input with respect to which the polynomial p(z) is being represented. Thus we may define the ratios um /um−1 and vm−1 /vm in more generality and produce convergent sequences to any root. Finally, we establish a relationship between the latter generalization and the Basic Sequence. This results in a deeper understanding of the Bernoulli Method as well as the existence of a precise error formula for its iterates. Let a ∈ C be any complex number that is not a root of p(z), and θ a root of p(z). From Taylor’s theorem we have (5.5)
n X p(n−i)(a) i=0
(n − i)!
(θ − a)n−i = 0 =
n X p(i) (a) i=0
i!
(θ − a)i .
Letting z = θ − a, from the left-hand-side of (5.5) we get (5.6)
n X p(n−i)(a) i=0
(n − i)!
z n−i = 0,
which can be viewed as the characteristic equation of the homogeneous recurrence relation n X p(n−i)(a) (5.7) um−i (a) = 0. (n − i)! i=0 Also by dividing the right-hand-side of (5.5) by (θ − a)n and letting w = 1/(θ − a), we get (5.8)
n X p(i)(a) i=0
i!
wn−i = 0.
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
11
The above can be viewed as the characteristic equation of the homogeneous recurrence relation n X p(i) (a) (5.9) vm−i (a) = 0. i! i=0
Definition 5.1. Let the Bernoulli Sequence associated with the farthest root from a be defined as um−1 (a) , αm (a) = a + um−2 (a) and the Bernoulli Sequence associated with the closet root to a be defined as vm−2 (a) βm (a) = a + p(a) . vm−1 (a) In particular, if a = 0 the corresponding Bernoulli Sequences coincides with the usual iterates of the Bernoulli Method. Theorem 5.2. Given p(z), let a be any complex number not in Sp (see (3.11)). Consider the Bernoulli Sequences with the initial conditions u1 (a) = 1, uj (a) = 0, for all j < 0; and v1 (a) = 1, vj (a) = 0, for all j < 0. We have 1. lim αm (a) = θ, lim βm (a) = θ′ , m→∞
m→∞
where θ and θ′ are root of p(z). 2. If all the roots of p(z) are distinct, then θ is the farthest root from a, and θ′ the nearest root to a. 3. dm (a) vm (a) = (−1)m m. p(a) 4. βm (a) = bm (a) = Bm (a), (i.e. the Bernoulli Sequence associated with the nearest root to a coincides with the Basic Sequence associated with a, and hence the pointwise evaluation of the Basic Family at a). Proof. Define (5.10)
Qa (z) =
n X p(i)(a) i=0
i!
b a (z) = z n−i , Q
n X p(n−i)(a) i=0
(n − i)!
z n−i.
b a (z) are, respectively, The set of distinct roots of Qa (z) and Q 1 1 b = {b η1 = (5.11) R = {η 1 = θ1 − a, . . . , ηt = θt − a}, R , . . . , ηbt = }. θ1 − a θt − a If a 6∈ Sp , then |¯ ηi | = 6 |¯ ηj |, for i 6= j. From Proposition 2.1 it follows that αm (a) converges to a root of p(z). Similarly, it follows that βm (a) converges to a root of p(z). From Proposition 2.1 it also follows that if all the roots are simple, then αm (a) converges to the farthest root from a, and βm (a) converges to the nearest root to a. Next we establish the equivalence of βm (a), and bm (a). First we prove by induction that dm (a) vm (a) = (−1)m m. p(a)
12
BAHMAN KALANTARI
The above is true for m = 0. Assume that it is true for all k = 0, . . . , m − 1. k Substituting vk (a) = (−1)k dk (a)/p(a) for k = 0, . . . , m − 1 in (5.9), we get (5.12)
p(a)vm (a) = −
n X p(i) (a) i=1
i!
(−1)m−i
dm−i(a)
m−i .
p(a)
Multiplying the above by p(a)m−1 and since (−1)−i+1 = (−1)i−1 we get (5.13)
p(a)m vm (a) = (−1)m
n X p(i) (a) i−1 (−1)i−1 dm−i (a)p(a) . i! i=1
But from (3.2) the right-hand-side of the above coincides with (−1)m dm (a). From this it easily follows that βm (a) = bm (a). The equivalence of βm (a) and Bm (a) was already established in Theorem 4.2
6. Determinantal Representation of Solution of Homogeneous Recurrence Relations and its Applications In this section we consider the general homogeneous linear recurrence relation with a singe initial condition, i.e., (6.1)
am = c1 am−1 + c2 am−2 + · · · + cn am−n , a0 = 1, , aj = 0, ∀ j < 0.
We first show that am can be represented as the determinant of a certain Toeplitz matrix. Next we show how this representation is useful in obtaining upper bounds on the modulus of am , and in some cases lower bounds. We will demonstrate the application of these results to the generalized Fibonacci numbers of order n. In the following theorems we assume that {am } satisfies the recurrence relation (6.1). Theorem 6.1. For each m ≥ 1 we have c1 −c2 . . . (−1)m−2 cm−1 .. .. 1 . . c1 . . .. .. am = det 0 1 . .. .. .. .. . . . 0 0 ... 1
(−1)m−1 cm
(−1)m−2 cm−1 .. . −c2 c1
where ci = 0, ∀ i > n. Proof. Given any set of numbers c1 , . . . , cn we can always define a polynomial p(z) of degree n and a complex number a such that (6.2)
ci = (−1)i−1 pi−1 (a)
p(i)(a) , i = 1, . . . , n i!
From Theorem 4.1 it follows that am coincides with the quantity Dm (a) corresponding to this polynomial (see definition of Dm (a)). Hence we obtain the claimed determinantal representation of am .
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
13
Theorem 6.2. For each m ≥ 1, we have m/2 n X . |ci | |am | ≤ 1 + i=1
Proof. From Hadamard inequality (see e.g. Marcus and Minc [17]), given an m × m matrix A = (aij ), we have Y 1/2 m X m det(A) ≤ . |aij |2 i=1 i=1
But this implies the claimed bound. The following is an immediate corollary which applies to a sequence known as the the generalized Fibonacci numbers (see Miles [19]): Corollary 6.3. Consider the Fibonacci sequence of order n defined via a single nonzero initial condition, i.e., (6.3)
Fm = Fm−1 + Fm−2 + · · · + Fm−n , F0 = 1, , Fj = 0, ∀ j < 0.
Then, we have Fm ≤ (n + 1)m/2 .
√ For n = 2 it is√well-known that Fm is bound above by [(1 + 5)/2]m. Note m that our bound is 3 and is comparable to the known bound. The problem of computing the determinant of a tridiagonal matrix that corresponds to the determinantal representation of the ordinary Fibonacci numbers (i.e. n = 2) appears as an exercise in Knuth [15] (see page 84). However, other than the determinantal representation given in the above corollary we are not aware of any other determinantal representation of the generalized Fibonacci numbers. Theorem 6.4. Suppose that there exists a positive number δ ∈ (0, 1) such that n X i=0
|ci| = (1 + δ)|c1 |.
Then we have m
m
|am | ≥ (1 − δ)⌈ 2 ⌉ (1 + δ)⌊ 2 ⌋ |c1 |m . The proof of Theorem 6.4 appears in Kalantari [13]. For the sake of completeness we reproduce this proof which in particular makes use of the following result: Theorem 6.5. (Kalantari [10]) Let A be an m × m complex matrix. Assume we are given positive numbers l and u such that if λ is an eigenvalue of A, then l ≤ |λ| ≤ u. Let mu − |trace(A)| κ = κ(l, u) ≡ ⌈ ⌉. u−l If |trace(A)| ≥ ml, then | det(A)| ≥ lκ um−κ .
14
BAHMAN KALANTARI
PROOF OF Theorem 6.4. Let Am be the matrix corresponding to am (see Theorem 6.1). The matrix Am is diagonally dominant. Since δ=
m 1 X |ci|, |c1 | i=0,i6=1
from Gerschgorin’s theorem it follows that if λ is an eigenvalue of Am , we have l = (1 − δ)|c1 | ≤ λ ≤ (1 + δ)|c1 | = u. Since |trace(∆m )| = m|c1 | ≥ m(1 − δ)|c1 |, Theorem 6.4 applies to Am , where κ=⌈ and we get
mu − |trace(∆m )| m(1 + δ)|a1 | − m|a1 | m ⌉=⌈ ⌉ = ⌈ ⌉, u−l (1 + δ)|a1 | − (1 − δ)|a1 | 2 ⌈ m2 ⌉ ⌊ m2 ⌋ |am | ≥ (1 − δ)|c1 | , (1 + δ)|c1 |
which coincides with the claimed lower bound. 7. Experimental Results Via Fractal Images In this section we present some of our experimentation with a graphic package capable of producing fractal images that depict the convergence behavior of the Basic Family Bm (z), as applied to a given polynomial, for a fixed m, as well as the convergence behavior of the Basic Sequence {Bm (a)}∞ m=2 . Consider a polynomial p(z) and a fixed natural number m ≥ 2. The basins of attraction of a root of p(z) with respect to the iteration function Bm (z) are regions in the complex plane such that given an initial point within them the corresponding sequence ak+1 = Bm (ak ) will converge to that root. The fractal nature and the images of the basins of attractions of Newton’s method are now quite familiar for some special polynomials such as p(z) = z n − 1. However, this is not the case as far as the behavior of other members of the Basic Family is concerned. In Figure 1 we consider a polynomial with randomly set of roots, depicted as dots. The figure shows the evolution of the basins of attraction of the roots to the Voronoi regions as m takes the values 2, 4, 10, and 50. In Figure 2-4 we show the basin of attractions for the polynomials p(z) = z 3 − 1, P (z) = z 4 − 1, and p(z) = z 9 − 1; corresponding to different values of m. Note that in all the figures as m increases the regions with chaotic behavior dramatically decrease to merely the boundaries of the Voronoi regions. Note also that, as predicted by the theory, for each root the basins of attractions converge to the Voronoi region of that root. The roots of p(z) = z n − 1 are the roots of unity and hence the Voronoi regions are completely symmetric. In all the figures in the case of m = 2, i.e. Newton’s method, the basins of attractions are chaotic. However, these regions rapidly improve by increasing m. For more vivid and color pictures of the above images as well as many more images we refer the reader to view www.cs.rutgers.edu/ kalantar/fractals. Moreover, in a forthcoming paper [14] we use this graphic package to examine the computational advantages of the higher order members of the Basic Family, as well as the use of the Basic Sequence in computing a single root or all the roots of a given polynomial.
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
Figure 1. Evolution of basins of attraction to Voronoi regions via Bm (z): random points, m = 2, 4, 10, 50 clockwise from upper left corner.
Figure 2. Evolution of basins of attraction to Voronoi regions via Bm (z): p(z) = z 3 − 1, m = 2, 3, 10, 50 clockwise from upper left corner.
15
16
BAHMAN KALANTARI
Figure 3. Evolution of basins of attraction to Voronoi regions via Bm (z): p(z) = z 4 − 1, m = 2, 3, 4, 50 clockwise from upper left corner.
Figure 4. Evolution of basins of attraction to Voronoi regions via Bm (z): p(z) = z 9 − 1, m = 2, 3, 10, 50 clockwise from upper left corner.
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
17
8. Concluding Remarks We have shown that given a complex polynomial p(z), to each complex number a one can associate a Basic Sequence, {Bm (a) = a−p(a)Dm−2 (a)/Dm−1 (a)}, where Dm (a) is dependent only on the normalized derivatives at a, and is also definable via a homogeneous linear recurrence relation of degree n. Except possibly for a set of measure zero (union of finite number of lines) for any other input a the corresponding Basic Sequence converges to a root of p(z). Under some regularity assumption (e.g. simplicity of all the roots) the Basic Sequence will necessarily converge to the nearest root of p(z). The evaluation of the Basic Sequence provides a very simple scheme for the approximation of a root of p(z). As an example given the complex polynomial p(z) = 1 + p1 z + · · · + pn z n , and a = 0, the corresponding Basic Sequence {Bm (0)} in many cases provides a very simple method for the approximation of a root of p(z). Even if the sequence {Dm (0)} has an infinite subsequence of zeros so that the Basic Sequence is not well-defined for all m, the Basic Sequence could still have subsequences that converge to a root of p(z). For instance consider p(z) = (z 2 + 2)(z − θ), say for θ = 1 and θ = 4. For θ = 1 the Basic Sequence {Bm (0)} converges to θ even-though the origin is equidistant to the complex roots. Although the root θ = 4 is not the nearest root to the origin, the Basic Sequence at a = 0 still has subsequences convergent to θ. If we now set p(z) to be p(z) = (z 2 + 2)(z − θ)(z − β)2 , where β is a real number, again at a = 0 we observe the convergence of the corresponding Basic Sequence, or a subsequence of it, to a root of p(z). The Basic Sequence can be efficiently computed. Its m-th term can be evaluated in O(n log n log m) arithmetic operations. For small m it can be computed in O(mn) operations. For each m, Dm (a) can also be represented as a special Toeplitz determinant. When the given input a is close enough to a simple root, for any m ≥ 2 the Basic Sequence can be turned into an iterative method of order m, by replacing the given input a with Bm (a), and repeating. There are many computational implication of these results. As an example one possible algorithm for polynomial root-finding is the following: for a given input a continue producing Bm (a) until the difference between Bm (a) and Bm+1 (a) is small. Then for a desirable m replace a with Bm (a) and repeat in order to obtain higher and higher accuracy. It is also possible to switch back and forth between the two schemes. An approach that is suggested by the properties of the Basic Sequence in order to compute all the roots of a given polynomial is as follows: first we obtain a rectangle that contains all the roots. Then by selecting a sparse number of gird point we can generate the corresponding Basic Sequence and thereby obtain good approximation to a subset of the roots. Then after deflation the same approach can be applied to find additional roots. This method could provide an excellent alternative to the method of Weyl [24] (see also Pan [20]). Experimentation with this method will be reported in [14]. The determinantal generalization of Taylor’s theorem (see Theorem 4.1) gives a mechanism for estimating the error |Bm (a) − θ|. The Basic Sequence can be even defined for functions that are not polynomial. In [12] we have used the pointwise evaluation of the Basic Family (the Basic Sequence in the terminology of this paper) to give new formulas for approximation of π and e. Although the definition of the Basic Sequence was motivated by the properties of the Basic Family of iteration
18
BAHMAN KALANTARI
functions, our new results have implied deeper properties of the Basic Family. Furthermore, the relationship between the Basic Sequence and the Basic Family has motivated new results on general homogeneous linear recurrence relations that are defined via a single nonzero initial condition. Finally, we have shown the equivalence of the Basic Sequence and the Bernoulli Method, once defined at an arbitrary point a. This connection gives a deeper understanding of the Bernoulli Method. In particular, using our determinantal generalization of Taylor’s theorem one can actually estimate the error in approximation of the extreme roots via the Bernoulli Method. In general our results should find many theoretical and practical applications in approximation of roots of polynomials and more general functions. In a forthcoming paper, [14], we will make use of fractal images to study the computational advantages of members of the Basic Family over Newton’s method, as well as the advantages in using the Basic Sequence in computing a single root or all the roots of a given polynomial. Acknowledgments. I would like to thank Mike Fredman for a communication. I specially would like to thank Mr. Zhengfang Zhu for having produced a beautiful graphic package that was used to produce the fractal figures given in this paper as well as those available at www.cs.rutgers.edu/ kalantar/fractals. References [1] D. Gries and G. Levin, Computing Fibonacci numbers (and similarly defined functions) in log time, IPL 11 (1980), 68-69. [2] F.B. Hildebrand, “Introduction to Numerical Analysis,” 2nd, ed., McGraw-Hill, New York, 1974. [3] A.S. Householder, “The Numerical Treatment of a Single Nonlinear Equation,” McGraw-Hill, New York, 1970. [4] C.M. Fiduccia, An efficient formula for linear recurrences, SIAM J. Comput 14 1985, 106-112. [5] B. Kalantari, and I. Kalantari, High order iterative methods for approximating square roots, BIT 36 (1996), 395-399. [6] B. Kalantari, I. Kalantari, and R. Zaare-Nahandi, A basic family of iteration functions for polynomial root finding and its characterizations, J. of Comput. Appl. Math., 80 (1997), 209-226. [7] B. Kalantari, On the order of convergence of a determinantal family of root-finding methods, BIT 39 (1999), 96-109. [8] B. Kalantari and J. Gerlach, Newton’s method and generation of a determinantal family of iteration functions, J. of Comp. and Appl. Math., 116 (2000), 195-200. [9] B. Kalantari, Generalization of Taylor’s theorem and Newton’s method via a new family of determinantal interpolation formulas and its applications, DCS-TR-328. To appear in J. of Comp. and Appl. Math. [10] B. Kalantari, A lower bound on determinants from linear programming, Technical Report DCS-TR 330, Department of Computer Science, Rutgers University, New Brunswick, New Jersey, 1997. [11] B. Kalantari and S. Park, A computational comparison of the first nine members of a determinantal family of rootfinding methods, DCS-TR-375. To appear in J. of Comp. and Appl. Math. [12] B. Kalantari, New formulas for approximations of Pi and other transcendental numbers, Technical Report DCS-TR-389, Department of Computer Science, Rutgers University, New Brunswick, New Jersey, 1999. To appear in Numerical Algorithms. [13] B. Kalantari, Approximation of polynomial root using a single input and the corresponding derivative values, Technical Report DCS-TR 369, Department of Computer Science, Rutgers University, New Brunswick, New Jersey, 1998. [14] B. Kalantari and Z. Zhu, A computational study of the Basic Family in locating all zeros of a complex polynomial using fractal images, forthcoming.
HOMOGENEOUS RECURRENCE RELATIONS AND ZEROS OF POLYNOMIALS
19
[15] D.E. Knuth, “The Art of Computer Programming, Volume I: Fundamental Algorithms,” 2nd Edition, Addison-Wesley, Reading, Mass, 1973. [16] H.T. Kung, A new upper bound on the complexity of derivative evaluation, Info. Proc. Lett., 2 (1974), 146-147. [17] M. Marcus and H. Minc, A Survey of Matrix Theory and Matrix Inequalities. Allyn and Bacon, Boston, Mass., 1964. [18] J.M. McNamee, A bibliography on root of polynomials, J. of Comp. and Appl. Math., 47 (1993), 391-394. [19] E.P. Miles, Generalized Fibonacci numbers and associated matrices, Amer. Math. Monthly, 67 (1960), 745-757. [20] V. Y. Pan, Solving a polynomial equation: some history and recent progress, SIAM Review, 39 (1997),187–220. [21] F.P. Preparata, “Computational Geometry An Introduction,” Springer-Verlag, New York, 1985. [22] S. Smale, The fundamental theorem of algebra and complexity theory, Bull. Amer. Math. Soc. 4 (1981), 1-36. [23] J.F. Traub, “Iterative Methods for the Solution of Equations,” Prentice Hall, Englewood Cliffs, New Jersey, 1964. [24] H. Weyl, Randbemerkungen zu Hauptproblemen der Mathematik. II. Fundamentalsatz der Algebra und Grundlagen der Mathematik, Math. Z. 20 (1924), 142–150. [25] T.J. Ypma, Historical development of Newton-Raphson method, SIAM Rev., 37 (1995), 531551. Department of Computer Science, Rutgers University, New Brunswick, New Jersey 08903 E-mail address:
[email protected]