- Email: [email protected]
On bounding zeros of analytic functions
Applied Mathematics and Computation 131 (2002) 383–396 www.elsevier.com/locate/amc
On bounding zeros of analytic functions M.A. Wolfe School of Mathematics and Statistics, University of St. Andrews, Fife, UK KY 16 9SS
Abstract An algorithm CZ for bounding simple and multiple zeros z 2 D C of an analytic function f : D C ! C in a compact set z D where D C is an open convex set is given. The algorithm is easy to implement using, for example Fortran 95. Some numerical results are presented. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Interval mathematics; Simple and multiple zeros of analytic functions
1. Introduction Let f : D C ! C be an analytic function and let D be an open convex set. Let z D be defined by z ¼ fz ¼ x1 þ ix2 j x1 6 x1 6 x1 ; x2 6 x1 6 x2 g; pffiffiffiffiffiffiffi where i ¼ 1 and x1 ; x2 2 R1 . An effective algorithm for determining computationally rigorous bounds on the zeros of f in z using variable-precision interval arithmetic is described in [1]. An interval algorithm for determining computationally rigorous bounds on the zeros of algebraic polynomials f in z is described in [2] with Pascal-XSC software. If f ðzÞ ¼ 0 is expressed as a pair of simultaneous equations Fk ðxÞ ¼ 0 T ðk ¼ 1; 2Þ where x ¼ ðx1 ; x2 Þ 2 R2 and z ¼ x1 þ ix2 then an algorithm which is described in [2] with Pascal-XSC software for determining computationally rigorous bounds on (simple) zeros of F in z may be used.
E-mail address: [email protected] (M.A. Wolfe). 0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 1 5 5 - 2
384
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
A large freely available Fortran package for determining computationally rigorous bounds on the zeros of F : Rn ! Rn with F 2 C 2 ðDÞ is described in [3] and may be used with n ¼ 2 to bound zeros of f in z. If, however, it is required to bound simple and multiple zeros of simple analytic functions, without having access to variable-precision interval arithmetic or to a Pascal-XSC compiler or of having the need of a large complicated software package then it is possible to obtain reasonably reliable results using an algorithm which is described in this paper. The results which may be obtained could be made computationally rigorous if one of the computational environments which supports rigorous interval arithmetic was to be used. The paper is arranged as follows. Section 2 contains the notation. Section 3 contains preliminary results on which the algorithm CZ which is described in Section 4 is based. Section 5 contains numerical results for a few simple analytic functions. Section 6 contains some remarks about the algorithm CZ. 2. Notation The fundamental interval mathematics that is required to understand this paper is in [4,5] (see also [6]). The notation that is needed is as follows. A real interval x ¼ ½x; x has infimum x 2 R1 and supremum x 2 R1 with x 6 x. The set of real intervals is denoted by IðR1 Þ. If X R1 then IðX Þ ¼ fx 2 IðR1 Þ j x X g. A real interval vector (a box) x ¼ ðxi Þ 2 IðRn Þ has infimum x ¼ ðxi Þ 2 Rn and supremum x ¼ ðxi Þ 2 Rn . A real interval matrix A ¼ ðai;j Þ 2 IðRmn Þ has infimum A ¼ ðai;j Þ 2 Rmn and supremum A ¼ ðai;j Þ 2 Rmn . The midpoint mðAÞ, the magnitude jAj and the width wðAÞ of A 2 IðRmn Þ are defined by: mðAÞ ¼ ððai;j þ ai;j Þ=2Þ 2 Rmn ; jAj ¼ ðmaxfjai;j j; jai;j jgÞ 2 Rmn ; wðAÞ ¼ ðai;j ai;j Þ 2 Rmn ; respectively. It is sometimes convenient to write A ¼ mðAÞ where A ¼ ðai;j Þ 2 IðRmn Þ. If mn A; B 2 IðR Þ then A B means that B < A 6 A < B, A B means that B < A 6 A 6 B or B 6 A 6 A < B and A B means that B 6 A 6 A 6 B where the usual componentwise partial order on Rmn is assumed. The spectral radius qðAÞ of A 2 Rnn is defined by qðAÞ ¼ maxfjkj j k 2 C ^ 9x 2 C n f0g; Ax ¼ kxg and qðAÞ 6 kAk1 [5]. The interval hull IðX Þ of X Rmn is defined by IðX Þ ¼ ½infðX Þ; supðX Þ: If A 2 IðRnn Þ is regular, so that 9A 1 ð8A 2 AÞ then
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
385
A 1 ¼ IfA 1 j A 2 Ag; and if A 2 IðRnn Þ is regular and b 2 IðRn Þ then AH b ¼ IfA 1 b j A 2 A ^ b 2 bg: 1 1 Thepbox ffiffiffiffiffiffiffi z ¼ x1 þ ix2 2 C has ReðzÞ ¼ x1 2 IðR Þ and ImðzÞ ¼ x2 2 IðR Þ, with i ¼ 1, and
IðCÞ ¼ fz ¼ x1 þ ix2 j x1 ; x2 2 Cg: The midpoint mðzÞ, the magnitude jzj and the width wðzÞ of z 2 IðCÞ are defined by: mðzÞ ¼ mðx1 Þ þ imðx2 Þ; jzj ¼ jx1 j þ jx2 j; wðzÞ ¼ wðx1 Þ þ wðx2 Þ: If F : Rn ! R1 is a differentiable function then oi F ðx1 ; . . . ; xn Þ ¼
@ F ðx1 ; . . . ; xn Þ @xi
ði ¼ 1; . . . ; nÞ:
3. Preliminaries This section contains some fundamental ideas and theorems which are needed in Section 4 and which are based on results from [7,5]. It is shown in [5] that if A 2 IðR22 Þ is regular and b 2 IðR2 Þ then AH b ¼ x where x 2 IðR2 Þ may be computed as follows. Let 8 a2;2 =a1;2 ð0 62 a1;2 Þ; > > < a1;2 =a2;2 ð0 62 a2;2 Þ; ð3:1Þ q¼ a =a ð0 62 a1;1 Þ; > > : 2;1 1;1 a1;1 =a2;1 ð0 62 a2;1 Þ: Then
and
! 8 b1 q b2 > b1 q b2 > > > < I a1;1 q a2;1 [ a1;1 q a2;1 ! x1 ¼ > b b b q 1 b2 q > 1 2 >I > [ : a1;1 a2;1 q a1;1 a2;1 q ! 8 b2 qb1 > b 2 qb1 > > > < I a2;2 qa1;2 [ a2;2 qa1;2 ! x2 ¼ > qb2 b1 qb2 b1 > > > : I qa a [ qa a 2;2 1;2 2;2 1;2
ð0 62 a1;2 Þ; ð3:2Þ ð0 62 a2;2 Þ;
ð0 62 a1;1 Þ; ð3:3Þ ð0 62 a2;1 Þ:
386
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
Furthermore if A 2 IðR22 Þ is regular then [5] AH b A 1 b:
ð3:4Þ
The following three theorems are due essentially to Alefeld [7]. Minor modifications have been made to the theorems from [7] for use with the algorithm which is described in Section 4. Theorem 3.1. Let F : D R2 ! R2 be a given function with F 2 C 1 ðDÞ where D is an open convex set, let x 2 IðDÞ, and let F0 ðÞ be a regular inclusion isotonic interval extension of F 0 ðÞ in x. If ~x 2 x then the following statements hold: (i) if x 2 x and F ðx Þ ¼ 0 H then x 2 ~x F0 ðxÞ F ð~xÞ; H 0 (ii) if f~x F ðxÞ F ð~xÞg \ x ¼ ; then F has no zero in x; (iii) if f~x F0 ðxÞH F ð~xÞg x then F has a unique zero in x. Proof. By (3.4) F0 ðxÞH F ð~xÞ F0 ðxÞ 1 F ð~xÞ. The remainder of the proof follows that of Theorem 1 in [7]. Let the hypotheses of Theorem 3.1 hold and let the sequence ðxðkÞ Þ be generated from xð0Þ ¼ x and for k P 0 ~xðkÞ 2 xðkÞ ; H
NðxðkÞ Þ ¼ ~xðkÞ F0 ðxðkÞ Þ F ð~xðkÞ Þ; ðkþ1Þ
x
ðkÞ
ð3:5Þ
ðkÞ
¼ Nðx Þ \ x :
If Nðxð0Þ Þ xð0Þ then by Theorem 3.1 F has a unique zero x 2 xð0Þ and x 2 xðkÞ ð8k P 0Þ. Theorem 3.2. If the hypotheses of Theorem 3.1 hold, if xð0Þ ¼ x, if Nðxð0Þ Þ xð0Þ , if qðAÞ < 1 where 1
A ¼ wðF0 ðxð0Þ Þ ÞjF0 ðxð0Þ Þj
ð3:6Þ
and if ðxðkÞ Þ is generated from (3.5) then x 2 xðkÞ ð8k P 0Þ and xðkÞ ! x ðk ! 1Þ where x is the unique zero of f in xð0Þ . If ~xðkÞ ¼ mðxðkÞ Þ ð8k P 0Þ then qðAÞ < 1 may be replaced with qðAÞ < 2.
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
387
Proof. The proof follows part of the proof of Theorem 2 in [7]. By (3.4) 1
H
F0 ðxð0Þ Þ F ð~xð0Þ Þ ¼ IfF 0 ðxÞ F ð~xð0Þ Þ j F 0 ðxÞ 2 F0 ðxð0Þ Þg 1
IfF 0 ðxÞ j F 0 ðxÞ 2 F0 ðxð0Þ ÞgF ð~xð0Þ Þ ¼ F0 ðxð0Þ Þ 1 F ð~xð0Þ Þ;
ð3:7Þ
so ð8k P 0Þ H
wðNðxðkÞ ÞÞ ¼ wðF0 ðxðkÞ Þ F ð~xðkÞ ÞÞ 6 wðF0 ðxð0Þ Þ 1 F ð~xðkÞ ÞÞ 1
6 wðF0 ðxð0Þ Þ ÞjF0 ðxð0Þ ÞjwðxðkÞ Þ 6 Akþ1 wðxð0Þ Þ; whence xðkÞ ! 0 ðk ! 1Þ.
Theorem 3.3. If the hypotheses of Theorem 3.2 hold and ð8x xð0Þ Þ there exists a > 0 such that kwðF0 ðxÞ 1 Þk1 6 akwðxÞk1
ð3:8Þ
then there exists c > 0 such that ð8k P 0Þ 2
kwðxðkþ1Þ Þk1 6 ckwðxðkÞ Þk1 :
ð3:9Þ
Proof. The proof is similar to that of Theorem 3 in [7]. 1
wðxðkþ1Þ Þ 6 wðF0 ðxðkÞ Þ ÞjF0 ðxð0Þ ÞjwðxðkÞ Þ; whence 1
kwðxðkþ1Þ Þk1 6 kwðF0 ðxðkÞ Þ Þk1 kF0 ðxð0Þ Þk1 kwðxðkÞ Þk1 whence (3.9) follows with c ¼ akF0 ðxð0Þ Þk1 .
Let f : D C ! C be an analytic function where D is an open convex set and let z 2pD. f satisfies the Cauchy–Riemann relations [8] at z ¼ x1 þ ix2 ffiffiffiffiffiffiThen ffi where i ¼ 1, so that if f ðzÞ ¼ uðx1 ; x2 Þ þ ivðx1 ; x2 Þ ð3:10Þ then o1 uðx1 ; x2 Þ ¼ o2 vðx1 ; x2 Þ
ð3:11Þ
o2 uðx1 ; x2 Þ ¼ o1 vðx1 ; x2 Þ:
ð3:12Þ
and Furthermore f 0 ðzÞ ¼ o1 uðx1 ; x2 Þ þ io1 vðx1 ; x2 Þ:
ð3:13Þ
388
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
It follows from (3.11) and (3.12) that if f 0 : IðDÞ ! IðCÞ is a continuous inclusion isotonic interval extension of f 0 : D ! C and z 2 IðDÞ then ð8z; ~z 2 zÞ f ðzÞ 2 f ð~zÞ þ f 0 ðzÞðz ~zÞ:
ð3:14Þ
This result has been used in [1,9] and is used in the algorithm which is described in Section 4. Let g : IðDÞ ! IðCÞ be defined by gðzÞ ¼ ~z f ð~zÞ=f 0 ðzÞ; where ~z 2 z 2 IðDÞ. Let F : IðR2 Þ ! IðR2 Þ be defined by F1 ðxÞ ¼ uðxÞ;
ð3:15Þ
F2 ðxÞ ¼ vðxÞ;
where u : IðR2 Þ ! IðR2 Þ and v : IðR2 Þ ! IðR2 Þ are continuous inclusion isotonic interval extensions of u : D ! C and v : D ! C, respectively. Then by (3.11) and (3.12) o1 uðxÞ o1 vðxÞ F0 ðxÞ ¼ ð3:16Þ o1 vðxÞ o1 uðxÞ and if F0 ðxÞ is regular, if ~z ¼ ~x1 þ i~x2 and if ~x ¼ ð~x1 ; ~x2 ÞT 2 x then by (3.4), (3.5), (3.11) and (3.12) N1 ðxÞ þ iN2 ðxÞ gðzÞ: ð0Þ
ð0Þ
Therefore if the sequence ðzðkÞ Þ is generated from zð0Þ ¼ x1 þ ix2 kP0
and for
~zðkÞ ¼ mðzðkÞ Þ; gðzðkÞ Þ ¼ ~zðkÞ f ð~zðkÞ Þ=f 0 ðzðkÞ Þ;
ð3:17Þ
zðkþ1Þ ¼ gðzðkÞ Þ \ zðkÞ ; and the sequence ðxðkÞ Þ is generated from (3.5) with ~xðkÞ ¼ mðxðkÞ Þ then ðkÞ ðkÞ x1 þ ix2 zðkÞ . On the other hand, if F0 ðxÞ is regular and GðxÞ ¼ ~x HðxÞF ð~xÞ; where HðxÞ ¼
1 ðfo1 uðxÞg2 þ fo2 uðxÞg2 Þ
then ð8x 2 xÞ F 0 ðxÞ
1
2 F0 ðxÞ
G1 ðxÞ þ iG2 ðxÞ ¼ gðzÞ:
1
o1 uðxÞ o1 vðxÞ
o1 vðxÞ o1 uðxÞ
ð3:18Þ
HðxÞ and by (3.10)–(3.13) ð3:19Þ
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
389
H
Furthermore Theorem 3.1 is valid when HðxÞ replaces F0 ðxÞ , and if ðxðkÞ Þ is generated from xð0Þ ¼ x and for k P 0 ~xðkÞ ¼ mðxðkÞ Þ; GðxðkÞ Þ ¼ ~xðkÞ HðxðkÞ ÞF ð~xðkÞ Þ;
ð3:20Þ
xðkþ1Þ ¼ GðxðkÞ Þ \ xðkÞ and A ¼ wðHðxð0Þ ÞÞjF 0 ð~xð0Þ Þj
ð3:21Þ 0
1
then by Theorem 3.2 with GðxÞ and HðxÞ replacing NðxÞ and F ðxÞ , respectively, it follows that if Gðxð0Þ Þ xð0Þ and qðAÞ < 2 then xðkÞ ! x ðk ! 1Þ where x is the unique zero of F in xð0Þ . Furthermore, by (3.19) if ðzðkÞ Þ is generated from (3.17) and ðxðkÞ Þ is generated from (3.20) then ð8k P 0Þ ðkÞ ðkÞ x1 þ ix2 ¼ zðkÞ and zðkÞ ! z ¼ x1 þ ix2 2 zð0Þ ðk ! 1Þ where z is the unique zero of f in z. Finally, by Theorem 3.3 with GðxÞ and HðxÞ replacing NðxÞ and 1 F0 ðxÞ , respectively, and A defined by (3.21) with qðAÞ < 2 it follows that ðzðkÞ Þ converges quadratically onto z if there exists a > 0 such that ð8z zð0Þ Þ wð1=f 0 ðzÞÞ 6 awðzÞ. The existence and uniqueness of z 2 D may be established by using the following results, even when z is not simple. Let Lþ : IðDÞ R22 R2 ! IðR2 Þ be defined by Lþ ðx; Y ; yÞ ¼ y Y fF ðyÞ þ ð1=2ÞF00 ðxÞðx yÞðx yÞg; 0
þ
ð3:22Þ 00
0
where y ¼ mðxÞ, Y ¼ F ðyÞ is the generalized inverse of F ðyÞ and F : IðDÞ ! IðR222 Þ is a continuous inclusion isotonic interval extension of F 00 : D R2 ! R222 with F00i;j;‘ ðxÞ ¼ oj o‘ Fi ðxÞ
ði; j; ‘ ¼ 1; 2Þ:
Theorem 3.4. If Lþ ðx; Y ; yÞ x, if rankðY Þ ¼ 2 and if wðLþ ðx; Y ; yÞÞ < wðxÞ then 9x 2 Lþ ðx; Y ; yÞ such that F ðx Þ ¼ 0 and x is the unique zero of F in x. Furthermore if Lþ ðx; Y ; yÞ \ x ¼ ; then x contains no zero of F. Proof. The proof is in [10].
Let / : D R2 ! R1 be defined by 2
2
/ðxÞ ¼ fuðxÞg þ fvðxÞg ;
ð3:23Þ
where u and v are as in (3.10). Then ðf ðz Þ ¼ 0Þ () ð/ðx Þ ¼ 0Þ where z ¼ x1 þ ix2 . Let Lþ : IðDÞ R12 R21 ! IðR1 Þ be defined by Lþ ðx; Y ; yÞ ¼ y Y f/ðyÞ þ ð1=2Þ/00 ðxÞðx yÞðx yÞg;
ð3:24Þ þ
where /00 ðxÞ ¼ ðoi oj /ðxÞÞ 2 IðR22 Þ, y ¼ mðxÞ and Y ¼ /0 ðyÞ .
390
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
Theorem 3.5. If Lþ ðx; Y ; yÞ x, if rankðY Þ ¼ 1 and if wðLþ ðx; Y ; yÞÞ < wðxÞ then 9x 2 Lþ ðx; Y ; yÞ such that /ðx Þ ¼ 0 and x is the unique zero of / in x. ~ ¼ Lþ ðx; Y ; yÞ then x 2 x ~ and 9c > 0 such that Furthermore if x 2
kwð~ xÞk1 6 ckwðxÞk1 : Proof. The proof is in [10].
Note that if Y 2 R22 then rankðY Þ ¼ 2 if and only if detðY Þ ¼ 2 and if Y 2 R12 then rankðY Þ ¼ 1 if and only if Y 6¼ 0. 4. The Algorithm CZ An algorithm CZ which determines bounds of at least prescribed width on simple and multiple zeros of analytic functions f : C ! C in a box z 2 IðCÞ is described in this section. The algorithm CZ consists of seven sub-algorithms A1 ; . . . ; A7 . The box z is enqueued into a list L1 (initially empty) which is to contain sub-boxes of z which are to be processed by CZ. Algorithm A1 proceeds as follows. If L1 is not empty then a box z is dequeued from L1 . If 0 2 fðzÞ then mðzÞ, f ðmðzÞÞ and f 0 ðzÞ are computed and an improved bound ~z on ff ðzÞ j z 2 zg is computed from (3.14). If now 0 2 ~z then H F0 ðxÞ F ð~xÞ is computed from (3.1)–(3.3) with x1 þ ix2 ¼ ~z, ~x ¼ mðxÞ and F defined by (3.15). H If F0 ðxÞ is not regular so that F0 ðxÞ F ð~xÞ cannot be computed then x is bisected along its side of greatest length to obtain boxes zð1Þ and zð2Þ . If wðzð1Þ Þ 6 e where e > 0 is given as data then zð1Þ is enqueued into a list L2 of so-called ‘narrow’ boxes each of which might contain a zero of f. If wðzð1Þ Þ > e then zð1Þ is enqueued into L1 . The box zð2Þ is treated similarly. H ^ is computed from If F0 ðxÞ F ð~xÞ can be computed then x ^ ¼ ~x F0 ðxÞH F ð~xÞ: x ^ then x is bisected. If x ^ x then by Theorem 3.1 9x 2 x ^ such that If x x ^1 þ i^ F ðx Þ ¼ 0 and x is the unique zero of F in x, so ^z ¼ x x2 is enqueued into a list L3 of so-called ‘safe’ boxes. Each box ^z in L3 contains a unique zero of f but wð^zÞ could be greater than e. ^ \ x 6¼ ; then ^z ¼ ð^ If x x \ xÞ1 þ ið^ x \ xÞ2 is enqueued into L2 if wð^zÞ 6 e and is enqueued into L1 if wð^zÞ > e. When L1 becomes empty and L3 is not empty then algorithm A2 is used, the purpose of which is to obtain from L3 boxes z each of which contain a unique zero of f and for which wðzÞ 6 e. Algorithm A2 proceeds as follows. If L3 6¼ ; then z is dequeued from L3 . If wðzÞ 6 e then z is enqueued into a list L4 of safe, narrow boxes (i.e., boxes z such that wðzÞ 6 e and which contain a unique zero of f ).
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
391
If wðzÞ > e then the sequence ðxðkÞ Þ is computed from (3.5) and (3.1)–(3.3) ð0Þ ð0Þ with z ¼ x1 þ ix2 and ð8k P 0Þ ~xðkÞ ¼ mðxðkÞ Þ, until for some k^ P 0 wð^zÞ 6 e ðk^Þ ðk^Þ where ^z ¼ x1 þ ix2 . If wðzÞ is sufficiently small then by Theorem 3.2 xðkÞ ! x ðk ! 1Þ so ^z is a safe narrow box which is enqueued into L4 . When L3 becomes empty and L2 6¼ ; then Algorithm A3 is used, the purpose of which is to obtain from L2 boxes z each of which is narrow and safe if such exist. If e > 0 is small then for each box z in L2 wðzÞ might be very small and it is necessary to use epsilon inflation [11] to determine from z a box containing a unique zero of f, if such a zero exists. Algorithm A3 is easier to communicate if it is expressed in (Fortran 95) pseudocode. On exit from A3 the list L5 contains boxes which are narrow, but for which epsilon inflation has been unsuccessful, and the list L4 contains boxes which are narrow and safe. Algorithm A3 Data: kMAX ¼ 3, ‘MAX ¼ 2, eps ¼ 0:1,L2 loop1 : do ‘¼0 if ðL2 ¼¼ ;Þ exit loop1 dequeue z from L2 z ¼ ½mðzÞ; mðzÞ loop2 : do k ¼ 0; zero_exists ¼ false loop3 : do k ¼ k þ 1; y ¼ inflatedðzÞ; y~ ¼ mðyÞ if ð0 2 f 0 ðyÞÞ exit loop2 z ¼ y~ f ð~ y Þ=f 0 ðyÞ if ðz yÞ then zero_exists ¼ true; exit loop2 ! Epsilon inflation has succeeded. end if if ðk ¼¼ kMAX Þ then eps ¼ 5 eps; exit loop3 ! Increase eps. end if end do loop3 ‘¼‘þ1 if ð‘ ¼¼ ‘MAX Þ then ! Epsilon inflation has failed for this z. eps ¼ 0:1; exit loop2 end if end do loop2 if ðzero_existsÞ then enqueue z into L4 else enqueue z into L5
392
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
end if end do loop1 function inflated ðzÞ; result y !g is the smallest positive machine number if ðwðReðzÞÞ= ¼ 0Þ then ReðyÞ ¼ ReðzÞ þ wðReðzÞÞ½ eps; eps else ReðyÞ ¼ ReðzÞ þ ½ g; g end if if ðwðImðzÞÞ= ¼ 0Þ then ImðyÞ ¼ ImðzÞ þ wðImðzÞÞ½ eps; eps else ImðyÞ ¼ ImðzÞ þ ½ g; g end if end function inflated end Algorithm A3 When A3 has terminated the lists L1 , L2 and L3 are empty. However, L5 is either empty or contains narrow and possibly safe boxes, some of which could be almost identical. Therefore algorithm A4 is used to merge L5 by determining the hulls of those boxes which are not disjoint to obtain the list L6 of boxes each of which could contain more than one possibly multiple zero of f. When L6 has been determined, algorithm A5 is used to obtain from L6 boxes ~z which are either narrow and safe, or for which this is uncertain. Algorithm A5 proceeds as follows. A box z is dequeued from L6 . Then Lþ is applied to F defined by (3.15) to obtain ~z such that either 0 62 fð~zÞ (in which case ~z is rejected), or, using Theorem 3.4 such that ~z contains a unique zero of f (in which case ~z is enqueued into L4 Þ or such that it is uncertain whether ~z contains a unique zero of f (in which case ~z is enqueued into the list L7 (originally empty)). When L6 is empty, L4 contains boxes z which are safe and almost certainly narrow even though L5 has been merged. However, L4 could contain a large number of almost identical boxes, so L4 is merged using algorithm A6 to obtain the list L8 , after which L4 is empty. When L4 is empty Lþ together with epsilon inflation is applied to / defined by (3.23) in algorithm A7 to obtain from L8 boxes ~z. If, using Theorem 3.5 it is established that ~z contains a unique zero of f then ~z is enqueued into L4 . Otherwise ~z is enqueued into L7 . On termination of CZ, L4 contains narrow safe boxes and L7 contains boxes which might contain zeros of f. If the width wðzÞ of the initial box z in which it is required to bound zeros of f is large then it has been found to be desirable to subdivide z before beginning
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
393
Sm
execution of CZ to obtain a set fzð1Þ ; . . . ; zðmÞ g such that z ¼ i¼1 zðiÞ where m ¼ m1 m2 according to the following subdivision algorithm. Initially L1 ¼ ;. Subdivision algorithm Data: z ¼ x1 þ ix2 , m1 > 0, m2 > 0, L1 ¼ ; h1 ¼ wðx1 Þ=m1 ; h2 ¼ wðx2 Þ=m2 do j ¼ 0; m1 ðjÞ x 1 ¼ x 1 þ j h1 end do do k ¼ 0; m2 ðkÞ x 2 ¼ x 2 þ k h2 end do do j ¼ 0; m1 1 ðjÞ ðjþ1Þ y1 ¼ ½x1 ; x1 do k ¼ 0; m2 1 ðkÞ ðkþ1Þ y2 ¼ ½x2 ; x2 y1 þ iy2 ¼ ðy1 þ iy2 Þ \ z if ð0 2 fðy1 þ iy2 ÞÞ then enqueue y into L1 end if end do end do end Subdivision algorithm A plan of CZ follows. program CZ use interval_arithmetic ! Interval arithmetic library use lists ! List manipulation procedures use czvp ! Global variables and parameters use czfns ! Test functions use czlib ! Procedures needed in subroutines A1 –A7 call read_the_data call initialize ! Initialize global variables and parameters call subdivision_algorithm (optional) call A1 ! Process L1 to obtain L2 and L3 call A2 ! Process L3 to obtain L4 call A3 ! Process L2 to obtain L5 and possibly part of L4 call A4 ! Merge L5 to obtain L6 call A5 ! Process L6 to obtain L7 and possibly part of L4 call A6 ! Merge L4 to obtain L8 call A7 ! Process L8 to replace L4 and possibly to add to L7 call write_the_output end program CZ
394
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
5. Numerical results The algorithm CZ has been implemented on a PC with dual 400 MHz Pentium II processors and 512 MB RAM using a NAG Fortran 95 compiler release 4.0 (185). In order to demonstrate the value of CZ for users who have no access to state-of-the-art environments such as are provided by Pascal-XSC, a simple interval arithmetic library was used to obtain the numerical results which are reported in this section. The upwardly and downwardly rounded machine numbers " x and # x of the machine number x were obtained as follows. If x > 0 then " x ¼ xð1 þ Þ þ g and # x ¼ xð1 Þ g; if x ¼ 0 then " x ¼# x ¼ 0; if x < 0 then" x ¼ xð1 Þ þ g and # x ¼ xð1 þ Þ g. Here > 0 is such that 1 and 1 þ are consecutive machine numbers and g is the smallest positive machine number [5]. This rounding was applied to the values of the standard functions which were obtained using the Fortran 95 compiler (with allowance for the periodicity of sin and cos) as well as to the basic interval arithmetic. This primitive foundation for interval arithmetic appears to be reasonably reliable for very small-scale problems, and is easily implementable. The following nine examples illustrate the performance of CZ. Example 1. f ðzÞ ¼ expðzÞ þ z; z ¼ ½ 1; 0 þ i½ 1; 1; z " 0:567143 . . . þ i0: Example 2. f ðzÞ ¼ coshðzÞ 1=2; z ¼ ½ 2; 2 þ i½ 2; 2; z ¼ 0 # ip=3: Example 3. f ðzÞ ¼ z3 6z2 þ 11z 6; z ¼ ½ 3; 3 þ i½ 3; 3; z1 ¼ 1 þ i0;
z2 ¼ 2 þ i0;
Example 4. f ðzÞ ¼ sinhðzÞ i; z ¼ ½ 2; 2 þ i½ 2; 2; z ¼ 0 þ ip=2:
z3 ¼ 3 þ i0:
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
395
Example 5. f ðzÞ ¼ sinðzÞ 2; z ¼ ½ 2; 2 þ i½ 2; 2; pffiffiffi z ¼ p=2 # i logð2 þ 3Þ: Example 6. f ðzÞ ¼ z4 1; z ¼ ½ 2; 2 þ i½ 2; 2; z1 ¼ 1 þ i0; z2 ¼ 0 þ i1;
z3 ¼ 1 þ i0;
z4 ¼ 0 i1:
Example 7. f ðzÞ ¼ sinðz2 Þ; z ¼ ½ 1; 1 þ i½ 1; 1; z ¼ 0 þ i0: Example 8. 2 f ðzÞ ¼ ðz 1Þ sinðzÞ; z ¼ ½ 2; 2 þ i½ 2; 2; z1 ¼ 0 þ i0; z2 ¼ 1 þ i0: Example 9. f ðzÞ ¼ ðz 1Þ2 ðexpðzÞ 1Þ; z ¼ ½ 1; 2 þ i½ 1; 1; z1 ¼ 0 þ i0; z2 ¼ 1 þ i0: In Table 1 n1 , n2 and n3 are the numbers of evaluations of f, f 0 and f 00 , respectively, n4 and n5 are the numbers of boxes remaining in the lists L4 and L7 , Table 1 Results from CZ with e ¼ 10 8 Ex.
n1
n2
n3
n4
n5
n6
n7
T (s)
1 2 3 4 5 6 7 8 9
36 164 1469 4923 54 584 2144 2337 1578
18 68 734 3019 22 240 789 840 565
1 2 3 7 2 4 4 3 3
1 1 3 1 2 4 1 2 2
0 0 0 0 0 0 0 0 0
2 2 3 5 2 2 2 2 2
2 2 3 5 2 2 2 2 2
0.00 0.01 0.05 0.34 0.00 0.02 0.10 0.13 0.09
396
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
respectively, on exit from CZ, n6 and n7 are the numbers m1 and m2 needed for the subdivision algorithm and T is the total CPU time in seconds, measured using the Fortran 95 standard procedure CPU_TIME. 6. Remarks about CZ The algorithm CZ has proved to be successful when applied to analytic functions having simple structure such as Examples 1–9 and initial boxes z for which wðzÞ is not too large. Examples 4, 7, 8 and 9 have multiple zeros, with f ðz Þ ¼ f 0 ðz Þ ¼ 0: this accounts for the large numbers of evaluations of f and of f 0 required by A1 and A2 in which several bisections and applications of epsilon inflation are needed. However, even for these examples, CPU times are small and in each case the widths of the enclosures of the zeros obtained are orders of magnitude less than e. The algorithm CZ would probably be much more effective were it to be implemented in Pascal-XSC and would certainly be completely reliable. Furthermore the winding-number idea used by Schaefer [1] could be incorporated into CZ in order to eliminate sub-boxes z for which 0 62 fðzÞ but could prove computationally expensive. References [1] M.J. Schaefer, Precise zeros of analytic functions using interval arithmetic, Interval Computations 4 (1993) 22–39. [2] R. Hammer, M. Hocks, U. Kulisch, D. Ratz, Numerical Toolbox for Verified Computing I, Springer, Berlin, 1993. [3] R.B. Kearfott, Rigorous Global Search: Continuous Problems, Kluwer Academic Publishers, Dordrecht, 1996. [4] G. Alefeld, J. Herzberger, Introduction to Interval Computations, Academic Press, London, 1983. [5] A. Neumaier, Interval Methods for Systems of Equations, Cambridge University Press, Cambridge, 1990. [6] M.A. Wolfe, Interval mathematics, algebraic equations and optimization, Journal of Computational and Applied Mathematics 124 (2000) 263–280. [7] G. Alefeld, Inclusion methods for systems of nonlinear equations – the interval newton method and modifications, in: J. Herzberger (Ed.), Topics in Validated Computations, Elsevier, Amsterdam, 1994. [8] R.V. Churchill, J.W. Brown, R.F. Verhey, Complex Variables and Applications, McGrawHill, New York, 1974. [9] M.A.Wolfe, Interval methods for algebraic equations, in: R.E. Moore (Ed.), International Workshop on the R^ ole of Interval Methods in Scientific Computing, Academic Press, London, 1988. [10] M.A. Wolfe, On Bounding Solutions of Underdetermined Systems, Reliable Computing 7 (2001) 195–207. [11] G. Mayer, Epsilon inflation in verification algorithms, Journal of Computational and Applied Mathematics 60 (1995) 147–169.
On bounding zeros of analytic functions M.A. Wolfe School of Mathematics and Statistics, University of St. Andrews, Fife, UK KY 16 9SS
Abstract An algorithm CZ for bounding simple and multiple zeros z 2 D C of an analytic function f : D C ! C in a compact set z D where D C is an open convex set is given. The algorithm is easy to implement using, for example Fortran 95. Some numerical results are presented. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Interval mathematics; Simple and multiple zeros of analytic functions
1. Introduction Let f : D C ! C be an analytic function and let D be an open convex set. Let z D be defined by z ¼ fz ¼ x1 þ ix2 j x1 6 x1 6 x1 ; x2 6 x1 6 x2 g; pffiffiffiffiffiffiffi where i ¼ 1 and x1 ; x2 2 R1 . An effective algorithm for determining computationally rigorous bounds on the zeros of f in z using variable-precision interval arithmetic is described in [1]. An interval algorithm for determining computationally rigorous bounds on the zeros of algebraic polynomials f in z is described in [2] with Pascal-XSC software. If f ðzÞ ¼ 0 is expressed as a pair of simultaneous equations Fk ðxÞ ¼ 0 T ðk ¼ 1; 2Þ where x ¼ ðx1 ; x2 Þ 2 R2 and z ¼ x1 þ ix2 then an algorithm which is described in [2] with Pascal-XSC software for determining computationally rigorous bounds on (simple) zeros of F in z may be used.
E-mail address: [email protected] (M.A. Wolfe). 0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 1 5 5 - 2
384
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
A large freely available Fortran package for determining computationally rigorous bounds on the zeros of F : Rn ! Rn with F 2 C 2 ðDÞ is described in [3] and may be used with n ¼ 2 to bound zeros of f in z. If, however, it is required to bound simple and multiple zeros of simple analytic functions, without having access to variable-precision interval arithmetic or to a Pascal-XSC compiler or of having the need of a large complicated software package then it is possible to obtain reasonably reliable results using an algorithm which is described in this paper. The results which may be obtained could be made computationally rigorous if one of the computational environments which supports rigorous interval arithmetic was to be used. The paper is arranged as follows. Section 2 contains the notation. Section 3 contains preliminary results on which the algorithm CZ which is described in Section 4 is based. Section 5 contains numerical results for a few simple analytic functions. Section 6 contains some remarks about the algorithm CZ. 2. Notation The fundamental interval mathematics that is required to understand this paper is in [4,5] (see also [6]). The notation that is needed is as follows. A real interval x ¼ ½x; x has infimum x 2 R1 and supremum x 2 R1 with x 6 x. The set of real intervals is denoted by IðR1 Þ. If X R1 then IðX Þ ¼ fx 2 IðR1 Þ j x X g. A real interval vector (a box) x ¼ ðxi Þ 2 IðRn Þ has infimum x ¼ ðxi Þ 2 Rn and supremum x ¼ ðxi Þ 2 Rn . A real interval matrix A ¼ ðai;j Þ 2 IðRmn Þ has infimum A ¼ ðai;j Þ 2 Rmn and supremum A ¼ ðai;j Þ 2 Rmn . The midpoint mðAÞ, the magnitude jAj and the width wðAÞ of A 2 IðRmn Þ are defined by: mðAÞ ¼ ððai;j þ ai;j Þ=2Þ 2 Rmn ; jAj ¼ ðmaxfjai;j j; jai;j jgÞ 2 Rmn ; wðAÞ ¼ ðai;j ai;j Þ 2 Rmn ; respectively. It is sometimes convenient to write A ¼ mðAÞ where A ¼ ðai;j Þ 2 IðRmn Þ. If mn A; B 2 IðR Þ then A B means that B < A 6 A < B, A B means that B < A 6 A 6 B or B 6 A 6 A < B and A B means that B 6 A 6 A 6 B where the usual componentwise partial order on Rmn is assumed. The spectral radius qðAÞ of A 2 Rnn is defined by qðAÞ ¼ maxfjkj j k 2 C ^ 9x 2 C n f0g; Ax ¼ kxg and qðAÞ 6 kAk1 [5]. The interval hull IðX Þ of X Rmn is defined by IðX Þ ¼ ½infðX Þ; supðX Þ: If A 2 IðRnn Þ is regular, so that 9A 1 ð8A 2 AÞ then
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
385
A 1 ¼ IfA 1 j A 2 Ag; and if A 2 IðRnn Þ is regular and b 2 IðRn Þ then AH b ¼ IfA 1 b j A 2 A ^ b 2 bg: 1 1 Thepbox ffiffiffiffiffiffiffi z ¼ x1 þ ix2 2 C has ReðzÞ ¼ x1 2 IðR Þ and ImðzÞ ¼ x2 2 IðR Þ, with i ¼ 1, and
IðCÞ ¼ fz ¼ x1 þ ix2 j x1 ; x2 2 Cg: The midpoint mðzÞ, the magnitude jzj and the width wðzÞ of z 2 IðCÞ are defined by: mðzÞ ¼ mðx1 Þ þ imðx2 Þ; jzj ¼ jx1 j þ jx2 j; wðzÞ ¼ wðx1 Þ þ wðx2 Þ: If F : Rn ! R1 is a differentiable function then oi F ðx1 ; . . . ; xn Þ ¼
@ F ðx1 ; . . . ; xn Þ @xi
ði ¼ 1; . . . ; nÞ:
3. Preliminaries This section contains some fundamental ideas and theorems which are needed in Section 4 and which are based on results from [7,5]. It is shown in [5] that if A 2 IðR22 Þ is regular and b 2 IðR2 Þ then AH b ¼ x where x 2 IðR2 Þ may be computed as follows. Let 8 a2;2 =a1;2 ð0 62 a1;2 Þ; > > < a1;2 =a2;2 ð0 62 a2;2 Þ; ð3:1Þ q¼ a =a ð0 62 a1;1 Þ; > > : 2;1 1;1 a1;1 =a2;1 ð0 62 a2;1 Þ: Then
and
! 8 b1 q b2 > b1 q b2 > > > < I a1;1 q a2;1 [ a1;1 q a2;1 ! x1 ¼ > b b b q 1 b2 q > 1 2 >I > [ : a1;1 a2;1 q a1;1 a2;1 q ! 8 b2 qb1 > b 2 qb1 > > > < I a2;2 qa1;2 [ a2;2 qa1;2 ! x2 ¼ > qb2 b1 qb2 b1 > > > : I qa a [ qa a 2;2 1;2 2;2 1;2
ð0 62 a1;2 Þ; ð3:2Þ ð0 62 a2;2 Þ;
ð0 62 a1;1 Þ; ð3:3Þ ð0 62 a2;1 Þ:
386
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
Furthermore if A 2 IðR22 Þ is regular then [5] AH b A 1 b:
ð3:4Þ
The following three theorems are due essentially to Alefeld [7]. Minor modifications have been made to the theorems from [7] for use with the algorithm which is described in Section 4. Theorem 3.1. Let F : D R2 ! R2 be a given function with F 2 C 1 ðDÞ where D is an open convex set, let x 2 IðDÞ, and let F0 ðÞ be a regular inclusion isotonic interval extension of F 0 ðÞ in x. If ~x 2 x then the following statements hold: (i) if x 2 x and F ðx Þ ¼ 0 H then x 2 ~x F0 ðxÞ F ð~xÞ; H 0 (ii) if f~x F ðxÞ F ð~xÞg \ x ¼ ; then F has no zero in x; (iii) if f~x F0 ðxÞH F ð~xÞg x then F has a unique zero in x. Proof. By (3.4) F0 ðxÞH F ð~xÞ F0 ðxÞ 1 F ð~xÞ. The remainder of the proof follows that of Theorem 1 in [7]. Let the hypotheses of Theorem 3.1 hold and let the sequence ðxðkÞ Þ be generated from xð0Þ ¼ x and for k P 0 ~xðkÞ 2 xðkÞ ; H
NðxðkÞ Þ ¼ ~xðkÞ F0 ðxðkÞ Þ F ð~xðkÞ Þ; ðkþ1Þ
x
ðkÞ
ð3:5Þ
ðkÞ
¼ Nðx Þ \ x :
If Nðxð0Þ Þ xð0Þ then by Theorem 3.1 F has a unique zero x 2 xð0Þ and x 2 xðkÞ ð8k P 0Þ. Theorem 3.2. If the hypotheses of Theorem 3.1 hold, if xð0Þ ¼ x, if Nðxð0Þ Þ xð0Þ , if qðAÞ < 1 where 1
A ¼ wðF0 ðxð0Þ Þ ÞjF0 ðxð0Þ Þj
ð3:6Þ
and if ðxðkÞ Þ is generated from (3.5) then x 2 xðkÞ ð8k P 0Þ and xðkÞ ! x ðk ! 1Þ where x is the unique zero of f in xð0Þ . If ~xðkÞ ¼ mðxðkÞ Þ ð8k P 0Þ then qðAÞ < 1 may be replaced with qðAÞ < 2.
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
387
Proof. The proof follows part of the proof of Theorem 2 in [7]. By (3.4) 1
H
F0 ðxð0Þ Þ F ð~xð0Þ Þ ¼ IfF 0 ðxÞ F ð~xð0Þ Þ j F 0 ðxÞ 2 F0 ðxð0Þ Þg 1
IfF 0 ðxÞ j F 0 ðxÞ 2 F0 ðxð0Þ ÞgF ð~xð0Þ Þ ¼ F0 ðxð0Þ Þ 1 F ð~xð0Þ Þ;
ð3:7Þ
so ð8k P 0Þ H
wðNðxðkÞ ÞÞ ¼ wðF0 ðxðkÞ Þ F ð~xðkÞ ÞÞ 6 wðF0 ðxð0Þ Þ 1 F ð~xðkÞ ÞÞ 1
6 wðF0 ðxð0Þ Þ ÞjF0 ðxð0Þ ÞjwðxðkÞ Þ 6 Akþ1 wðxð0Þ Þ; whence xðkÞ ! 0 ðk ! 1Þ.
Theorem 3.3. If the hypotheses of Theorem 3.2 hold and ð8x xð0Þ Þ there exists a > 0 such that kwðF0 ðxÞ 1 Þk1 6 akwðxÞk1
ð3:8Þ
then there exists c > 0 such that ð8k P 0Þ 2
kwðxðkþ1Þ Þk1 6 ckwðxðkÞ Þk1 :
ð3:9Þ
Proof. The proof is similar to that of Theorem 3 in [7]. 1
wðxðkþ1Þ Þ 6 wðF0 ðxðkÞ Þ ÞjF0 ðxð0Þ ÞjwðxðkÞ Þ; whence 1
kwðxðkþ1Þ Þk1 6 kwðF0 ðxðkÞ Þ Þk1 kF0 ðxð0Þ Þk1 kwðxðkÞ Þk1 whence (3.9) follows with c ¼ akF0 ðxð0Þ Þk1 .
Let f : D C ! C be an analytic function where D is an open convex set and let z 2pD. f satisfies the Cauchy–Riemann relations [8] at z ¼ x1 þ ix2 ffiffiffiffiffiffiThen ffi where i ¼ 1, so that if f ðzÞ ¼ uðx1 ; x2 Þ þ ivðx1 ; x2 Þ ð3:10Þ then o1 uðx1 ; x2 Þ ¼ o2 vðx1 ; x2 Þ
ð3:11Þ
o2 uðx1 ; x2 Þ ¼ o1 vðx1 ; x2 Þ:
ð3:12Þ
and Furthermore f 0 ðzÞ ¼ o1 uðx1 ; x2 Þ þ io1 vðx1 ; x2 Þ:
ð3:13Þ
388
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
It follows from (3.11) and (3.12) that if f 0 : IðDÞ ! IðCÞ is a continuous inclusion isotonic interval extension of f 0 : D ! C and z 2 IðDÞ then ð8z; ~z 2 zÞ f ðzÞ 2 f ð~zÞ þ f 0 ðzÞðz ~zÞ:
ð3:14Þ
This result has been used in [1,9] and is used in the algorithm which is described in Section 4. Let g : IðDÞ ! IðCÞ be defined by gðzÞ ¼ ~z f ð~zÞ=f 0 ðzÞ; where ~z 2 z 2 IðDÞ. Let F : IðR2 Þ ! IðR2 Þ be defined by F1 ðxÞ ¼ uðxÞ;
ð3:15Þ
F2 ðxÞ ¼ vðxÞ;
where u : IðR2 Þ ! IðR2 Þ and v : IðR2 Þ ! IðR2 Þ are continuous inclusion isotonic interval extensions of u : D ! C and v : D ! C, respectively. Then by (3.11) and (3.12) o1 uðxÞ o1 vðxÞ F0 ðxÞ ¼ ð3:16Þ o1 vðxÞ o1 uðxÞ and if F0 ðxÞ is regular, if ~z ¼ ~x1 þ i~x2 and if ~x ¼ ð~x1 ; ~x2 ÞT 2 x then by (3.4), (3.5), (3.11) and (3.12) N1 ðxÞ þ iN2 ðxÞ gðzÞ: ð0Þ
ð0Þ
Therefore if the sequence ðzðkÞ Þ is generated from zð0Þ ¼ x1 þ ix2 kP0
and for
~zðkÞ ¼ mðzðkÞ Þ; gðzðkÞ Þ ¼ ~zðkÞ f ð~zðkÞ Þ=f 0 ðzðkÞ Þ;
ð3:17Þ
zðkþ1Þ ¼ gðzðkÞ Þ \ zðkÞ ; and the sequence ðxðkÞ Þ is generated from (3.5) with ~xðkÞ ¼ mðxðkÞ Þ then ðkÞ ðkÞ x1 þ ix2 zðkÞ . On the other hand, if F0 ðxÞ is regular and GðxÞ ¼ ~x HðxÞF ð~xÞ; where HðxÞ ¼
1 ðfo1 uðxÞg2 þ fo2 uðxÞg2 Þ
then ð8x 2 xÞ F 0 ðxÞ
1
2 F0 ðxÞ
G1 ðxÞ þ iG2 ðxÞ ¼ gðzÞ:
1
o1 uðxÞ o1 vðxÞ
o1 vðxÞ o1 uðxÞ
ð3:18Þ
HðxÞ and by (3.10)–(3.13) ð3:19Þ
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
389
H
Furthermore Theorem 3.1 is valid when HðxÞ replaces F0 ðxÞ , and if ðxðkÞ Þ is generated from xð0Þ ¼ x and for k P 0 ~xðkÞ ¼ mðxðkÞ Þ; GðxðkÞ Þ ¼ ~xðkÞ HðxðkÞ ÞF ð~xðkÞ Þ;
ð3:20Þ
xðkþ1Þ ¼ GðxðkÞ Þ \ xðkÞ and A ¼ wðHðxð0Þ ÞÞjF 0 ð~xð0Þ Þj
ð3:21Þ 0
1
then by Theorem 3.2 with GðxÞ and HðxÞ replacing NðxÞ and F ðxÞ , respectively, it follows that if Gðxð0Þ Þ xð0Þ and qðAÞ < 2 then xðkÞ ! x ðk ! 1Þ where x is the unique zero of F in xð0Þ . Furthermore, by (3.19) if ðzðkÞ Þ is generated from (3.17) and ðxðkÞ Þ is generated from (3.20) then ð8k P 0Þ ðkÞ ðkÞ x1 þ ix2 ¼ zðkÞ and zðkÞ ! z ¼ x1 þ ix2 2 zð0Þ ðk ! 1Þ where z is the unique zero of f in z. Finally, by Theorem 3.3 with GðxÞ and HðxÞ replacing NðxÞ and 1 F0 ðxÞ , respectively, and A defined by (3.21) with qðAÞ < 2 it follows that ðzðkÞ Þ converges quadratically onto z if there exists a > 0 such that ð8z zð0Þ Þ wð1=f 0 ðzÞÞ 6 awðzÞ. The existence and uniqueness of z 2 D may be established by using the following results, even when z is not simple. Let Lþ : IðDÞ R22 R2 ! IðR2 Þ be defined by Lþ ðx; Y ; yÞ ¼ y Y fF ðyÞ þ ð1=2ÞF00 ðxÞðx yÞðx yÞg; 0
þ
ð3:22Þ 00
0
where y ¼ mðxÞ, Y ¼ F ðyÞ is the generalized inverse of F ðyÞ and F : IðDÞ ! IðR222 Þ is a continuous inclusion isotonic interval extension of F 00 : D R2 ! R222 with F00i;j;‘ ðxÞ ¼ oj o‘ Fi ðxÞ
ði; j; ‘ ¼ 1; 2Þ:
Theorem 3.4. If Lþ ðx; Y ; yÞ x, if rankðY Þ ¼ 2 and if wðLþ ðx; Y ; yÞÞ < wðxÞ then 9x 2 Lþ ðx; Y ; yÞ such that F ðx Þ ¼ 0 and x is the unique zero of F in x. Furthermore if Lþ ðx; Y ; yÞ \ x ¼ ; then x contains no zero of F. Proof. The proof is in [10].
Let / : D R2 ! R1 be defined by 2
2
/ðxÞ ¼ fuðxÞg þ fvðxÞg ;
ð3:23Þ
where u and v are as in (3.10). Then ðf ðz Þ ¼ 0Þ () ð/ðx Þ ¼ 0Þ where z ¼ x1 þ ix2 . Let Lþ : IðDÞ R12 R21 ! IðR1 Þ be defined by Lþ ðx; Y ; yÞ ¼ y Y f/ðyÞ þ ð1=2Þ/00 ðxÞðx yÞðx yÞg;
ð3:24Þ þ
where /00 ðxÞ ¼ ðoi oj /ðxÞÞ 2 IðR22 Þ, y ¼ mðxÞ and Y ¼ /0 ðyÞ .
390
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
Theorem 3.5. If Lþ ðx; Y ; yÞ x, if rankðY Þ ¼ 1 and if wðLþ ðx; Y ; yÞÞ < wðxÞ then 9x 2 Lþ ðx; Y ; yÞ such that /ðx Þ ¼ 0 and x is the unique zero of / in x. ~ ¼ Lþ ðx; Y ; yÞ then x 2 x ~ and 9c > 0 such that Furthermore if x 2
kwð~ xÞk1 6 ckwðxÞk1 : Proof. The proof is in [10].
Note that if Y 2 R22 then rankðY Þ ¼ 2 if and only if detðY Þ ¼ 2 and if Y 2 R12 then rankðY Þ ¼ 1 if and only if Y 6¼ 0. 4. The Algorithm CZ An algorithm CZ which determines bounds of at least prescribed width on simple and multiple zeros of analytic functions f : C ! C in a box z 2 IðCÞ is described in this section. The algorithm CZ consists of seven sub-algorithms A1 ; . . . ; A7 . The box z is enqueued into a list L1 (initially empty) which is to contain sub-boxes of z which are to be processed by CZ. Algorithm A1 proceeds as follows. If L1 is not empty then a box z is dequeued from L1 . If 0 2 fðzÞ then mðzÞ, f ðmðzÞÞ and f 0 ðzÞ are computed and an improved bound ~z on ff ðzÞ j z 2 zg is computed from (3.14). If now 0 2 ~z then H F0 ðxÞ F ð~xÞ is computed from (3.1)–(3.3) with x1 þ ix2 ¼ ~z, ~x ¼ mðxÞ and F defined by (3.15). H If F0 ðxÞ is not regular so that F0 ðxÞ F ð~xÞ cannot be computed then x is bisected along its side of greatest length to obtain boxes zð1Þ and zð2Þ . If wðzð1Þ Þ 6 e where e > 0 is given as data then zð1Þ is enqueued into a list L2 of so-called ‘narrow’ boxes each of which might contain a zero of f. If wðzð1Þ Þ > e then zð1Þ is enqueued into L1 . The box zð2Þ is treated similarly. H ^ is computed from If F0 ðxÞ F ð~xÞ can be computed then x ^ ¼ ~x F0 ðxÞH F ð~xÞ: x ^ then x is bisected. If x ^ x then by Theorem 3.1 9x 2 x ^ such that If x x ^1 þ i^ F ðx Þ ¼ 0 and x is the unique zero of F in x, so ^z ¼ x x2 is enqueued into a list L3 of so-called ‘safe’ boxes. Each box ^z in L3 contains a unique zero of f but wð^zÞ could be greater than e. ^ \ x 6¼ ; then ^z ¼ ð^ If x x \ xÞ1 þ ið^ x \ xÞ2 is enqueued into L2 if wð^zÞ 6 e and is enqueued into L1 if wð^zÞ > e. When L1 becomes empty and L3 is not empty then algorithm A2 is used, the purpose of which is to obtain from L3 boxes z each of which contain a unique zero of f and for which wðzÞ 6 e. Algorithm A2 proceeds as follows. If L3 6¼ ; then z is dequeued from L3 . If wðzÞ 6 e then z is enqueued into a list L4 of safe, narrow boxes (i.e., boxes z such that wðzÞ 6 e and which contain a unique zero of f ).
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
391
If wðzÞ > e then the sequence ðxðkÞ Þ is computed from (3.5) and (3.1)–(3.3) ð0Þ ð0Þ with z ¼ x1 þ ix2 and ð8k P 0Þ ~xðkÞ ¼ mðxðkÞ Þ, until for some k^ P 0 wð^zÞ 6 e ðk^Þ ðk^Þ where ^z ¼ x1 þ ix2 . If wðzÞ is sufficiently small then by Theorem 3.2 xðkÞ ! x ðk ! 1Þ so ^z is a safe narrow box which is enqueued into L4 . When L3 becomes empty and L2 6¼ ; then Algorithm A3 is used, the purpose of which is to obtain from L2 boxes z each of which is narrow and safe if such exist. If e > 0 is small then for each box z in L2 wðzÞ might be very small and it is necessary to use epsilon inflation [11] to determine from z a box containing a unique zero of f, if such a zero exists. Algorithm A3 is easier to communicate if it is expressed in (Fortran 95) pseudocode. On exit from A3 the list L5 contains boxes which are narrow, but for which epsilon inflation has been unsuccessful, and the list L4 contains boxes which are narrow and safe. Algorithm A3 Data: kMAX ¼ 3, ‘MAX ¼ 2, eps ¼ 0:1,L2 loop1 : do ‘¼0 if ðL2 ¼¼ ;Þ exit loop1 dequeue z from L2 z ¼ ½mðzÞ; mðzÞ loop2 : do k ¼ 0; zero_exists ¼ false loop3 : do k ¼ k þ 1; y ¼ inflatedðzÞ; y~ ¼ mðyÞ if ð0 2 f 0 ðyÞÞ exit loop2 z ¼ y~ f ð~ y Þ=f 0 ðyÞ if ðz yÞ then zero_exists ¼ true; exit loop2 ! Epsilon inflation has succeeded. end if if ðk ¼¼ kMAX Þ then eps ¼ 5 eps; exit loop3 ! Increase eps. end if end do loop3 ‘¼‘þ1 if ð‘ ¼¼ ‘MAX Þ then ! Epsilon inflation has failed for this z. eps ¼ 0:1; exit loop2 end if end do loop2 if ðzero_existsÞ then enqueue z into L4 else enqueue z into L5
392
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
end if end do loop1 function inflated ðzÞ; result y !g is the smallest positive machine number if ðwðReðzÞÞ= ¼ 0Þ then ReðyÞ ¼ ReðzÞ þ wðReðzÞÞ½ eps; eps else ReðyÞ ¼ ReðzÞ þ ½ g; g end if if ðwðImðzÞÞ= ¼ 0Þ then ImðyÞ ¼ ImðzÞ þ wðImðzÞÞ½ eps; eps else ImðyÞ ¼ ImðzÞ þ ½ g; g end if end function inflated end Algorithm A3 When A3 has terminated the lists L1 , L2 and L3 are empty. However, L5 is either empty or contains narrow and possibly safe boxes, some of which could be almost identical. Therefore algorithm A4 is used to merge L5 by determining the hulls of those boxes which are not disjoint to obtain the list L6 of boxes each of which could contain more than one possibly multiple zero of f. When L6 has been determined, algorithm A5 is used to obtain from L6 boxes ~z which are either narrow and safe, or for which this is uncertain. Algorithm A5 proceeds as follows. A box z is dequeued from L6 . Then Lþ is applied to F defined by (3.15) to obtain ~z such that either 0 62 fð~zÞ (in which case ~z is rejected), or, using Theorem 3.4 such that ~z contains a unique zero of f (in which case ~z is enqueued into L4 Þ or such that it is uncertain whether ~z contains a unique zero of f (in which case ~z is enqueued into the list L7 (originally empty)). When L6 is empty, L4 contains boxes z which are safe and almost certainly narrow even though L5 has been merged. However, L4 could contain a large number of almost identical boxes, so L4 is merged using algorithm A6 to obtain the list L8 , after which L4 is empty. When L4 is empty Lþ together with epsilon inflation is applied to / defined by (3.23) in algorithm A7 to obtain from L8 boxes ~z. If, using Theorem 3.5 it is established that ~z contains a unique zero of f then ~z is enqueued into L4 . Otherwise ~z is enqueued into L7 . On termination of CZ, L4 contains narrow safe boxes and L7 contains boxes which might contain zeros of f. If the width wðzÞ of the initial box z in which it is required to bound zeros of f is large then it has been found to be desirable to subdivide z before beginning
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
393
Sm
execution of CZ to obtain a set fzð1Þ ; . . . ; zðmÞ g such that z ¼ i¼1 zðiÞ where m ¼ m1 m2 according to the following subdivision algorithm. Initially L1 ¼ ;. Subdivision algorithm Data: z ¼ x1 þ ix2 , m1 > 0, m2 > 0, L1 ¼ ; h1 ¼ wðx1 Þ=m1 ; h2 ¼ wðx2 Þ=m2 do j ¼ 0; m1 ðjÞ x 1 ¼ x 1 þ j h1 end do do k ¼ 0; m2 ðkÞ x 2 ¼ x 2 þ k h2 end do do j ¼ 0; m1 1 ðjÞ ðjþ1Þ y1 ¼ ½x1 ; x1 do k ¼ 0; m2 1 ðkÞ ðkþ1Þ y2 ¼ ½x2 ; x2 y1 þ iy2 ¼ ðy1 þ iy2 Þ \ z if ð0 2 fðy1 þ iy2 ÞÞ then enqueue y into L1 end if end do end do end Subdivision algorithm A plan of CZ follows. program CZ use interval_arithmetic ! Interval arithmetic library use lists ! List manipulation procedures use czvp ! Global variables and parameters use czfns ! Test functions use czlib ! Procedures needed in subroutines A1 –A7 call read_the_data call initialize ! Initialize global variables and parameters call subdivision_algorithm (optional) call A1 ! Process L1 to obtain L2 and L3 call A2 ! Process L3 to obtain L4 call A3 ! Process L2 to obtain L5 and possibly part of L4 call A4 ! Merge L5 to obtain L6 call A5 ! Process L6 to obtain L7 and possibly part of L4 call A6 ! Merge L4 to obtain L8 call A7 ! Process L8 to replace L4 and possibly to add to L7 call write_the_output end program CZ
394
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
5. Numerical results The algorithm CZ has been implemented on a PC with dual 400 MHz Pentium II processors and 512 MB RAM using a NAG Fortran 95 compiler release 4.0 (185). In order to demonstrate the value of CZ for users who have no access to state-of-the-art environments such as are provided by Pascal-XSC, a simple interval arithmetic library was used to obtain the numerical results which are reported in this section. The upwardly and downwardly rounded machine numbers " x and # x of the machine number x were obtained as follows. If x > 0 then " x ¼ xð1 þ Þ þ g and # x ¼ xð1 Þ g; if x ¼ 0 then " x ¼# x ¼ 0; if x < 0 then" x ¼ xð1 Þ þ g and # x ¼ xð1 þ Þ g. Here > 0 is such that 1 and 1 þ are consecutive machine numbers and g is the smallest positive machine number [5]. This rounding was applied to the values of the standard functions which were obtained using the Fortran 95 compiler (with allowance for the periodicity of sin and cos) as well as to the basic interval arithmetic. This primitive foundation for interval arithmetic appears to be reasonably reliable for very small-scale problems, and is easily implementable. The following nine examples illustrate the performance of CZ. Example 1. f ðzÞ ¼ expðzÞ þ z; z ¼ ½ 1; 0 þ i½ 1; 1; z " 0:567143 . . . þ i0: Example 2. f ðzÞ ¼ coshðzÞ 1=2; z ¼ ½ 2; 2 þ i½ 2; 2; z ¼ 0 # ip=3: Example 3. f ðzÞ ¼ z3 6z2 þ 11z 6; z ¼ ½ 3; 3 þ i½ 3; 3; z1 ¼ 1 þ i0;
z2 ¼ 2 þ i0;
Example 4. f ðzÞ ¼ sinhðzÞ i; z ¼ ½ 2; 2 þ i½ 2; 2; z ¼ 0 þ ip=2:
z3 ¼ 3 þ i0:
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
395
Example 5. f ðzÞ ¼ sinðzÞ 2; z ¼ ½ 2; 2 þ i½ 2; 2; pffiffiffi z ¼ p=2 # i logð2 þ 3Þ: Example 6. f ðzÞ ¼ z4 1; z ¼ ½ 2; 2 þ i½ 2; 2; z1 ¼ 1 þ i0; z2 ¼ 0 þ i1;
z3 ¼ 1 þ i0;
z4 ¼ 0 i1:
Example 7. f ðzÞ ¼ sinðz2 Þ; z ¼ ½ 1; 1 þ i½ 1; 1; z ¼ 0 þ i0: Example 8. 2 f ðzÞ ¼ ðz 1Þ sinðzÞ; z ¼ ½ 2; 2 þ i½ 2; 2; z1 ¼ 0 þ i0; z2 ¼ 1 þ i0: Example 9. f ðzÞ ¼ ðz 1Þ2 ðexpðzÞ 1Þ; z ¼ ½ 1; 2 þ i½ 1; 1; z1 ¼ 0 þ i0; z2 ¼ 1 þ i0: In Table 1 n1 , n2 and n3 are the numbers of evaluations of f, f 0 and f 00 , respectively, n4 and n5 are the numbers of boxes remaining in the lists L4 and L7 , Table 1 Results from CZ with e ¼ 10 8 Ex.
n1
n2
n3
n4
n5
n6
n7
T (s)
1 2 3 4 5 6 7 8 9
36 164 1469 4923 54 584 2144 2337 1578
18 68 734 3019 22 240 789 840 565
1 2 3 7 2 4 4 3 3
1 1 3 1 2 4 1 2 2
0 0 0 0 0 0 0 0 0
2 2 3 5 2 2 2 2 2
2 2 3 5 2 2 2 2 2
0.00 0.01 0.05 0.34 0.00 0.02 0.10 0.13 0.09
396
M.A. Wolfe / Appl. Math. Comput. 131 (2002) 383–396
respectively, on exit from CZ, n6 and n7 are the numbers m1 and m2 needed for the subdivision algorithm and T is the total CPU time in seconds, measured using the Fortran 95 standard procedure CPU_TIME. 6. Remarks about CZ The algorithm CZ has proved to be successful when applied to analytic functions having simple structure such as Examples 1–9 and initial boxes z for which wðzÞ is not too large. Examples 4, 7, 8 and 9 have multiple zeros, with f ðz Þ ¼ f 0 ðz Þ ¼ 0: this accounts for the large numbers of evaluations of f and of f 0 required by A1 and A2 in which several bisections and applications of epsilon inflation are needed. However, even for these examples, CPU times are small and in each case the widths of the enclosures of the zeros obtained are orders of magnitude less than e. The algorithm CZ would probably be much more effective were it to be implemented in Pascal-XSC and would certainly be completely reliable. Furthermore the winding-number idea used by Schaefer [1] could be incorporated into CZ in order to eliminate sub-boxes z for which 0 62 fðzÞ but could prove computationally expensive. References [1] M.J. Schaefer, Precise zeros of analytic functions using interval arithmetic, Interval Computations 4 (1993) 22–39. [2] R. Hammer, M. Hocks, U. Kulisch, D. Ratz, Numerical Toolbox for Verified Computing I, Springer, Berlin, 1993. [3] R.B. Kearfott, Rigorous Global Search: Continuous Problems, Kluwer Academic Publishers, Dordrecht, 1996. [4] G. Alefeld, J. Herzberger, Introduction to Interval Computations, Academic Press, London, 1983. [5] A. Neumaier, Interval Methods for Systems of Equations, Cambridge University Press, Cambridge, 1990. [6] M.A. Wolfe, Interval mathematics, algebraic equations and optimization, Journal of Computational and Applied Mathematics 124 (2000) 263–280. [7] G. Alefeld, Inclusion methods for systems of nonlinear equations – the interval newton method and modifications, in: J. Herzberger (Ed.), Topics in Validated Computations, Elsevier, Amsterdam, 1994. [8] R.V. Churchill, J.W. Brown, R.F. Verhey, Complex Variables and Applications, McGrawHill, New York, 1974. [9] M.A.Wolfe, Interval methods for algebraic equations, in: R.E. Moore (Ed.), International Workshop on the R^ ole of Interval Methods in Scientific Computing, Academic Press, London, 1988. [10] M.A. Wolfe, On Bounding Solutions of Underdetermined Systems, Reliable Computing 7 (2001) 195–207. [11] G. Mayer, Epsilon inflation in verification algorithms, Journal of Computational and Applied Mathematics 60 (1995) 147–169.