In this paper, we study the problem of how to locate the zeros and poles of a meromorphic function in a given region , by using the generalized argument principle. This class of problems can be traced back to the classic work of Harry Nyquist in 1932 cit:feb:09:16:37 , which determines the stability of a dynamical system by searching zeros in the upper-half plane, now known as Nyquist stability criterion and widely used in many fields. However, it can only tell us the number of zeros, without any information about their precise locations. Delves and Lyness suggested a more modern approach to locate the zeros of an analytic function cit:feb:09:18:14 . They construct a polynomial whose zeros coincide with the zeros of
, with its coefficients calculated via the momentsin Eq.(1). This work was extended to the problem of finding the zeros and poles of a meromorphic function by Abd-Elall et al. cit:feb:09:18:47 . The stability is hard to be guaranteed in this approach. Kravanja et al. cit:feb:09:19:28 generalized the approach of Delves and Lyness by using the theory of formal orthogonal polynomials and Prony’s methodcit:jan:27:11:02 . Their approach involves a complicated procedure of constructing the so-called Gram matrix, and refers the sensitivity to the condition of a Yule-Walker system.
The major drawback of these earlier works on this problem is a lack of rigorous theoretical analysis of the sensitivity and error estimation. From this perspective, a stable, easy-to-implement method is still in needed. The main motivation of this work is to find a simple and practical answer to the problem, based on an active restriction on the condition of the Prony system.
The remainder of this paper is organized as follows. In Section II we present the mathematical background of finding zeros and poles by argument principle, and review the algorithm referred to as Prony’s method. In Section III we analyse the sensitivity of this algorithm, and give an explicit expression for the condition number, which enables us to control the algorithm stability proactively. In Section IV we prove a theorem for our singular pencil corrupted by noise, then apply it to estimate the error level. Section V provides three numerical examples to support our results. We state our conclusions in Section VI.
Ii Argument principle method for locating the zeros and poles
As with the earlier works, we consider a meromorphic function in a closed complex domain , bounded by a Jordan curve . Assume has zeros and poles in but no zeros or poles on , from the well-known generalized argument principle in complex analysis, we have
where the summation is over all the zeros and poles of counted with their multiplicities , and is an analytic function in , which will be referred to as a probe function hereafter. Given the possibility to evaluate along , the argument principle method wishes to recover from a series of .
Clearly, the key point is to choose suitable probe functions . Although many sophisticated schemes have been suggested as mentioned above, the fundamental elements are of the form
from which general polynomial probe functions can be constructedcit:jan:24:18:24 . In the present work, for the purpose of practical applications, we restrict our attention to the scaled probe functions
here , and is the scaled variable,
with a shift and scaling radius will be determined later. By using the scaled variable, we redefine the moments as
Once the moments have been calculated, following the approach in cit:jan:24:16:28 , we can construct two Hankel matrices ,
Since an arbitrary Hankel matrix of finite rank admits a Vandermonde decomposition, these matrices can be factorized ascit:feb:10:12:17 :
is a Vandermonde matrix,
We then have the following theorems ( Theorem 3.2.1 and 3.2.2 of cit:jan:24:16:28 ):
Let be the number of zeros and poles, then for every .
Theorem 1 gives us the number of the zeros and poles .
Theorem 2 is actually a reformulated Prony’s methodcit:jan:27:11:02 , which has been extensively applied into many fields. The generalized eigenvalue problem, Eq.(10), can be solved by a QZ algorithm with operationscit:jan:27:15:00 , and this will give us the zeros and poles . Once are known, the multiplicities can be obtained through a Vandermonde system Eq.(5). As the multiplicities must be integers, this step is relatively robust. With such a procedure, we can recover all the information of zeros and poles in , at least in theory.
Iii Sensitivity Analysis
In this section, we discuss the sensitivity of the algorithm presented in the previous section. Since both Vandermonde and Hankel matrices can be very ill-conditioned, it’s impossible to assert the algorithm stability in any universal sense. However, we can focus on the questions that what factors will affect the sensitivity significantly in our case, and how to find a stable parameter interval in the practical application.
To simplify the problem and facilitate theoretical analysis, we assume (without loss of generality) that the complex domain is bounded by a circle with the center and radius , and there are zeros and poles in it. We take a similar approach to cit:jan:24:18:24 , and perturb the generalized eigenvalue problem, Eq.(10), as
where is a small number, are normalized perturbations with Hankel structures. Then to the first order, we have
for specific eigenvalue and its corresponding eigenvector , in consideration of the symmetry of Hankel matrices, we have
pre-multiply Eq.(12) by and take the norm, we get
will give the error estimation. In view of , we have
thus the sensitivity can be estimated by
By assuming , , and noticing that ,
and the well-known matrix norm of inverse of Vandermonde matrix cit:feb:03:10:43 ,
it’s easy to tell that small can improve the condition of the problem, therefore we should take . Now Eq.(16) becomes
with the condition number . Furthermore, if we define
here we have adopted the polar form . Therefore the condition number becomes
So we find small results in better conditioned problems, and set .
Beyond that, Eq.(25) indicates that the sensitivity of the algorithm presented in previous mainly depends on the number and positions of zeros and poles in . Although large and the clusters may make the condition dramatically worse, in the case of zeros and poles distributed uniformly around the unit circle, the algorithm can be stable, even for relatively large .
This inspires us to adopt a subdivision-transformation scheme to guarantee the algorithm stability. Such a scheme asks a basic numerical question: Once a set of zeros and poles is preliminarily calculated, how large is the condition number of the Prony system? If we encounter a condition number larger than the preassigned, we should subdivide the region into smaller subregions and do the transformation Eq.(4) in each subregion. and are chosen according to preliminary , such that the zeros and poles are unifromly distributed around the unit circle in new varible if possible. In more concrete terms, a recommended subdivision strategy is to construt subregions bounded by two concentric circles or squares, and a reasonable should make and close to , should be taken to avoid the appearance of clusters. Then the question is asked repeatedly, untill we get an acceptable condition number. Therefore, Eq.(25) gives us the possibility to decide when and how to subdivide the region and transform to new variable without any need for human supervision. Since the required computation time majorily depends on the number of subregions, and such a subdivision-transformation scheme can bear comparatively large with the premise of ensuring the algorithm stability, it can save time considerably.
Iv Error Analysis
As the generalized eigenvalue problem, Eq.(10), can be solved in well-condition, the algorithm still suffers from an inherent error source, i.e., the numerical errors arised from the contour integration Eq.(5).
Note that the numerical differentiation is time-consuming and error-prone, we perform an integration by parts in Eq.(5) and exploit the logarithmic derivative to avoid evaluating . In doing this, the main difficulty is the multivalueness of the complex logarithm. So in order to recognize the same branch of numerically, we rewrite
and keep track of the numerical calculated principal value and select appropriate s to ensure the extended argument is continuous with . Then Eq.(5) becomes
with and artificial starting-point and end-point of the contour integration. Introducing the complex logarithm in this way provides several benefits. First, it avoids the evaluating of . Second, since is known to be an integer, it can be specified very accurately as with the Nyquist stability criterion. Third, the complex logarithm is more robust against the overflow and underflow issues.
When the algorithm yields the final result , a very relevant question is: What’s the error level? The next theorem and the following discussion allows us to tackle this problem easily.
For small , the eigenvalues of the singular pencil which is corrupted by noise fall into two categories:
, i.e., the eigenvalues of the pencil , which are independent of ;
, which depend on .
First, we redefine Eq.(3) as
and denote the minor matrices of as and , respectively.
Then we obtain the decompositions:
which is just the result of removing -th row from a Vandermonde matrix . Analogously, we can get by removing -th column from a matrix .
Now we construct a new generalized eigenvalue problem . By noting that
we find that the pencil have the same eigenvalues with , i.e., are also solutions to . Furthermore, by reference to the previous section, we see their condition numbers are in the same order for small .
Then the cofactor expansion along a row gives
where represents the number of minor matrices . Since we have restricted our calculation to near the unit circle, we find around for small . It follows that are also eigenvalues of the corrupted pencil . Also notice the QZ-algorithm does not involve the problems of rank determination and matrix inversion, indeed we can get from the ill-conditioned pencil by a QZ-algorithm in well-condition.
This proves the theorem. ∎
Theorem.(3) provides a method to estimate the order of numerical errors at little cost. Specifically, let and be the same eigenvalue calculated, respectively, from and , then the error can be estimated by
Owing to cofactors involved and , implies a weighted average over the results of equivalent pencils, Eq.(IV).
V Numerical examples
In this section, we summarize the algorithm and apply it to three numerical examples to illustrate the main points of this paper.
The algorithm discussed above can be executed as follows:
Set a critical condition number and an error tolerance of the contour integration ;
Choose a region of interest and transform to the scaled variable using Eq.(4);
Construct Hankel matrices, determine their ranks and solve the pencil ;
Calculate the condition number using Eq.(25), if , subdivide into smaller subregions and go back to Step 2, otherwise, continue;
Solve to give the error estimation;
Obtain the multiplicities through the Vandermonde system, Eq.(5).
Now we present three practical examples to show the effectiveness of the algorithm. In particular, the first two examples give the zeros and poles of two specific meromorphic functions, while in the third example, we apply the algorithm to the nonlinear eigenvalue problem. All these examples demonstrate that Theorem.(3) estimates the order of error reasonably.
v.0.1 Example 1:
v.0.2 Example 2:
In this example, we consider the zeros of plasma dispersion functioncit:feb:27:11:02 , which occurs repeatedly in the kinetic theory of hot plasmas. This function is defined as
and as the analytic continuation of this for , with a symmetry property . By choosing the same as before and focusing on the zeros in the range , we get the simple zeros listed in Table.(2). Moreover, from Figure.(1) we see that, as grows, the zeros become denser and closer to , which is a typical characteristic of the error function. In fact, an alternative representation of is just
v.0.3 Example 3:
Nearly all of the eigenvalue problems in plasma physics are nonlinear. However, no general effective solution method is yet available. From the analysis given above, we observe that the algorithm described in this paper has the advantages of being stable and easy-to-parallel, together with the fast-increasing computational power, it might be used to handle these awkward problems. As an example, we solve a transcendental eigenvalue problemcit:feb:08:13:49 : , with
We list in Table.(3) a summary of the simple zeros in the range and their associated errors, with their distribution is depicted in Figure.(2). More plasma physics related nonlinear eigenvalue problems solved by the proposed algorithm will be reported in future publications.
In this paper we have presented a numerical method for locating the zeros and poles of a meromorphic function in a given region, with easy implementation. Sensitivity is analyzed in detail, and an explicit formulation for the condition number is found, which is the core of the algorithm. Such an explicit formulation is utilized to limit the condition number and maintain the algorithm stability. In addition, the proposed algorithm can estimate the numerical errors automatically, supported by Theorem.(3).
In summary, the proposed algorithm has the following advantages:
: It can locate the zeros and poles in a given region stably. Bad convergence may be encountered when handling complicated functions by Newton-type methods, due to the Newton fractal phenomenon;
- Suitability for parallelization
: Intricate problems can be solved accurately by massively-parallel calculations;
- Automatic result verification
: It can estimate the error level;
- Wide applicability
: Meromorphic functions are very common in physics.
Given these advantages, this algorithm may have many applications in physical problems, particularly those difficult for other efficient methods.
One of the authors (H. T. Chen ) would like to thank G. Y. Fu, M. Y. Yu, Z. Y. Qiu, Y. Xiao and H. S. Xie for useful conversations.
This work was supported by National Natural Science Foundation of China under Grant No. 11905097.
- (1) H. Nyquist, Bell Syst. Tech., 11, 126 (1932).
- (2) L. M. Delves and J. N. Lyness, Math. Comp., 21, 543 (1967).
- (3) L. F. Abd-Elall, L. M. Delves and J. K. Reid, in Numerical Methods for Nonlinear Algebraic Equations, edited by P. Rabinowitz (Gordon and Breach, 1970), pp 47-59.
- (4) P. Kravanja, T. Sakurai and M. Van Barel, BIT Numer. Math., 39, 646 (1999).
- (5) F. B. Hildebrand, Introduction to Numerical Analysis, second edn, (Dover, 1987), pp 457-462.
- (6) G. H. Golub, P. Milanfar and J. Varah, SIAM J. Sci. Comput., 21, 1222 (1999).
- (7) P. Kravanja and V. V. Barel, Computing the Zeros of Analytic Functions (Springer, 2000), pp 83-89.
- (8) D. L. Boley, F. T. Luk and D. Vandevoorde, Vandermonde Factorization of a Hankel Matrix, in Scientific Computing, (Springer, 1997), pp 27-39.
- (9) G. H. Golub and C. F. Van Loan, Matrix Computations, third edn, (Johns Hopkins University Press, 1996), pp 375-390.
- (10) W. Gautschi, Numer. Math., 4, 119 (1962).
- (11) B. D. Fried and S. D. Conte, The Plasma Dispersion Function, (Academic Press, 1961).
- (12) T. Betcke, N. J. Higham, V. Mehrmann, C. Schroder and F. Tisseur, NLEVP: A Collection of Nonlinear Eigenvalue Problems, (MIMS EPrint, 2011).