Computing zeros of analytic functions in the complex plane without using derivatives☆
Introduction
In several fields of computational science one is required to obtain numerical solutions to the problem where is an analytic function defined only on a restricted domain of the complex plane , and is bounded by a simply connected closed curve, Γ. We always assume that no root of actually lies on Γ; this is not a severe restriction because the contour can generally be chosen such that its does not pass through a root.
A complicating factor in such kind of problems is usually that, although is analytic on , either the function cannot be easily differentiated or it is such that the numerical calculation of its derivatives, in addition to repeated evaluation of the function, is computationally prohibitive. Such a situation is typical for the electromagnetic problems solved by moment methods where is represented as the determinant of a matrix, each element of which is expressed by integrals or series [1], [2], [3], [4], [5], [6], [7]. Similar problems arise in the study of intermediate energy electron scattering by atoms and molecules [8] where z represents the energy of the incident electron. Computation of the resonance parameters for the electron scattering requires the solution of where is defined with respect to a matrix. For each value of z the matrix, but not its derivative, may be calculated using large computer program suites [9]. In these circumstances, it is therefore necessary to employ solution methods which do not require derivative evaluations and this paper presents a computer program developed for that purpose. In order to rigorously test and validate our program, we demonstrate the verification results for functions whose roots are already known or can be obtained analytically.
Iterative methods [10] are often employed for the refinement of the numerical solution of (1), i.e. if we have, an approximation for a solution, then we may obtain its better estimate, by applying the iterative procedure [10] where is an interpolation function. This process is repeated until we satisfy the condition where ϵ is a user defined threshold. The distinction between various iteration methods essentially arises in the different forms of the interpolation function H used.
We have found Halley's method [10], [11] to be particularly computationally efficient and adopted it in this paper. Although the standard formulation of Halley's method involves derivatives, the use of finite differences, instead of the derivatives, maintains the same convergence rate [10], and we have implemented this finite difference approach in our work.
The set of values generated by (2) is an iterated map of the complex plane. In order to ensure convergence, the initial estimate must lie in the basin of attraction of the root with respect to H, that is where is the filled-in Julia Set for H and its Julia set [12]. However a priori, it is essentially impossible to know whether a given satisfies (4) or not. Furthermore this can only be determined after executing a sufficiently large number of iterations [13]. Alternatively, one needs a priori knowledge of the function properties and root location when picking . Yau and Ben-Israel [14] have shown that, at least in some instances, the Julia set for Halley iterations has a less chaotic shape that the equivalent for Newton iterations.
In applying an iteration technique directly to the stated problem, no account of the function analyticity is taken at the onset. It is only checked that . Clearly, if after a few iterations it appears that , one has no choice but to halt the iterations, as there may be no way of converging to a zero of located inside . Then the process should be restarted from the beginning with a different initial value of . The occurrence of this situation is a reflection of the fact that in general, the basic iterative methods do not guarantee convergence within the finite domain, especially when possesses multiple roots. It is therefore important that the root isolation step produces suitable values.
Indeed, making use of the Cauchy theorem, function can be expressed, at any internal point , by an integral of around the contour Γ which encloses . Previous work in this journal by Botten et al. [15] extended the technique of Delves and Lyness [16] which highlighted the power of using contour integrals to compute the points at which . This work of Botten et al. [15] employed circular contours which are traditionally used in complex analysis. However, when has branch cuts or has complicated shape which should be spanned without gaps and overlaps by several sub-domains, circular contours are hardly applicable. Therefore in this work we use rectangular contours. In particular, this obviates the need to perform co-ordinate transformations to circles, a process which is generally computationally expensive [17].
There is one situation however in which the contour integral method may encounter problems. This happens when the root or, indeed, multiple roots lie very close to the contour, a situation in which the evaluation of the contour integrals is extremely difficult. In some problems, there are additional roots lying nearby but outside the contour, Γ, so the simple idea of perturbing the contour slightly may be inappropriate.
From a practical computational perspective, the contour integral method will only produce estimates, though reasonably good approximations in general, to the exact roots. Limited accuracy of this approach is associated with the inherent errors of numerical quadratures and by the finite dynamic range of the floating-point arithmetic available on the particular computer processors used. Thus, if a computed estimate is such that then the iterative methods may be applied to refine the root. Splitting the numerical solution of (1) into two stages is of fundamental significance for reliable evaluation of multiple roots of . Thus, roots are separated or bracketed at first. Then each root estimate is refined independently inside the specified sub-domain of . We have implemented this two-stage approach in a Fortran 90 program where the Halley iterative algorithm has been adapted for the root refinement in the restricted domain containing a single root.
In the next section, we present a situation in computational electromagnetics which requires use of our algorithm. Next we present the analytical and computational details of our algorithm. The following section describes the format for inputing data to our Fortran computer program; section four presents some details of the implementation of our code and the following section describes the output from some test runs. The final section of the paper outlines our conclusions and some directions for future work.
Section snippets
The dispersion equation for a semiconductor layer
Silicon based integrated micro systems combine multiple active and passive components densely packed in layers of substrates and when operated in the terahertz and gigahertz range, such assemblies cause delay, cross-talk and distortion of transmitted signals. In order to minimize the losses, it is necessary to model and to compute electromagnetic wave propagation in multilayered structures [18].
One model, which is an extension of previous work [18], is illustrated in Fig. 1 where a layer of
The method
In our work we take to be a rectangular region in the complex plane and Γ to be its perimeter. For any analytical function defined on the rectangle, we implement a multi-step process to solve Eq. (1).
- (1)
Count the number of roots enclosed by the rectangle.
- (2)
If is greater than a user defined threshold then sub-divide the region into smaller rectangles, and repeat, until each sub-rectangle has at most roots within. This is similar to the approach taken by Dellnitz et al. [19].
- (3)
Implementation details
We have implemented the code in Fortran 90 in a traditional procedural fashion. and have not therefore used any of the object oriented facilities provided by the language We have made extensive use of the dynamic memory allocation facilities in the language to reduce to a minimum the number of fixed dimensions in the code and thereby making the code more adaptable to problem size and this is controlled by the user input. We have found this to be particularly useful when running multiple
Test runs
We have executed test runs on a on Dell Optiplex GX270 PC which has a 3.2 GHz Intel Pentium processor, 1 Gbyte of RAM and 80 Gbytes of hard disk. The operating system, executing on the PC, was the Fedora Core 2 version of Linux. Specifically, this means that the 2.6.5 kernel was used. The code was compiled using version 8.1 of the Intel Fortran 90 for Linux compiler. We have used the following compiler options in all test runs: Taken together, we believe
Concluding remarks
We have presented a new package, written in Fortran 90, which computes estimated solutions to the problem where is defined over a finite rectangular domain in the complex plane. Computation of is not performed but it is a fundamental requirement that be analytic within and on the boundary of . Our package is targeted at computing problems where the evaluation of the derivative is computationally expensive. Our package computes estimated solutions and the derivative
Acknowledgements
The Institute for Electronics, Communications and Information Technology (ECIT) is jointly funded by The Queen's University of Belfast, Invest Northern Ireland and The Department of Employment and Learning for Northern Ireland.
References (27)
- et al.
Comput. Phys. Comm.
(1983) - et al.
J. Comput. Appl. Math.
(2002) - et al.
J. Comput. Phys.
(1982) - et al.
Comput. Phys. Comm.
(2000) - et al.
Int. J. Numer. Model. Electron. Netw. Device Fields (UK)
(1995) - et al.
IEEE Trans. on Microwave Theory and Techniques
(1998) - et al.
ACES J.
(2001) - et al.
IEEE Trans. on Circuits and Systems I: Fundamental Theory and Applications
(2002) - et al.
IEEE Trans. on Microwave Theory and Techniques
(2002) - A.G. Schuchinsky, E. Brennan, V.F. Fusco, in: Fifth IEE International Conference on Computation in Electromagnetics—CEM...
Int. J. of RF and Computer-Aided Engineering
J. Phys. B: At. Mol. Opt. Phys.
Cited by (24)
CCOMP: An efficient algorithm for complex roots computation of determinantal equations
2018, Computer Physics CommunicationsCitation Excerpt :Outside core, other techniques are used to deliver the good initial point estimation. The argument principle method (APM) is quite popular [2,3,33–37], however, falls short when the region of investigation includes poles. The Abd-ellal–Delves–Reid (ADR) method takes into account the latter restriction [4,38].
Robust location of optical fiber modes via the argument principle method
2017, Computer Physics CommunicationsCitation Excerpt :This procedure is adopted in several implementations [6,8,9], but suffers from several drawbacks. Firstly, it can lead to failure if the contour passes through the vicinity of several roots [9]. However, the more pertinent disadvantage for our purposes is that the enlarged contour may cross a branch point or branch cut, which is fatal to the method.
Noise reduction by the application of an air-bubble curtain in offshore pile driving
2016, Journal of Sound and VibrationA frequency-independent and parallel algorithm for computing the zeros of strictly proper rational transfer functions
2016, Applied Mathematics and ComputationA self-adaptive complex root tracing algorithm for the analysis of propagation and radiation problem
2021, IEEE Transactions on Antennas and Propagation
- ☆
This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).