Finding Zeros of Analytic Functions and Local Eigenvalue Analysis Using Contour Integral Method in Examples

A numerical method for computing zeros of analytic complex functions is presented. It relies on Cauchy’s residue theorem and the method of Newton’s identities, which translates the problem to finding zeros of a polynomial. In order to stabilize the numerical algorithm, formal orthogonal polynomials are employed. At the end the method is adapted to finding eigenvalues of a matrix pencil in a bounded domain in the complex plane. This work is based on a series of papers of Professor Sakurai and collaborators. Our aim is to make their work available by means of a systematic study of properly chosen examples.


Introduction
We present a numerical method for computing zeros of an analytic complex function.The method is further extended to finding eigenvalues of a matrix pencil.This is an essential problem in many areas of engineering such as analysis of mechanical vibrations, electrical networks, optical waveguides, or in quantum chemistry to name a few.Apart from traditional algorithms, cf.[10], which often compute all eigenvalues and rely on costly diagonalization of the system matrix, the presented contour integral method: • computes the roots/eigenvalues only in the region of interest, • and requires solutions to forward perturbed (possibly nonlinear) systems, which can often be achieved at a cost proportional to the problem size.
The problem to find a root of a complex function f is usually solved by fixed-point iterations, where the related mapping is contractive in a neighbourhood of the root/fixed-point.If f is smooth the Newton method is in a sense the best choice.Up to some extent the root can be eliminated from the function and the process repeats.Perhaps the main bottleneck is that the convergence of the Newton method requires to start close enough to the root.
In this paper we follow a conceptually different approach.We assume f to be analytical and we calculate moments of F (z) := f (z) • f (z) −1 along a given curve in the complex plane.By Cauchy's residue theorem the moments give us a complete information about the poles of F , i.e., the roots of f inside the curve.This method dates back to the pioneering work of Delves and Lyness [2].In the method, which is also referred to as the method of Newton's identities, one searches for roots of a polynomial the coefficients of which are unstable.A remedy is proposed by Kravanja, Sakurai and van Barel [5].They stabilize the method by formal orthogonal polynomials.
The contour integral method was further extended by Sakurai and Sugiura [11] towards computing local eigenvalues of matrix pencil (A, B) by finding poles of matrix-valued function F (z) := V T (zB − A) −1 V. Later it was reformulated by Ikegami, Sakurai, and Nagashima [4] using the resolvent theory, see [6], [7] and a more accurate algorithm was proposed.One of the algorithmic parameters is the number of distinct eigenvalues, for which a good estimate is given in [12].The method was also generalized to nonlinear eigenvalue problems [1].This paper relies on Master thesis of the first author [13].Essentialy, it is a compilation of a series of papers on contour integral method towards the local eigenvalue analysis.Our aim is to make the method available by documenting its functionality and the necessity of particular ingrediences on properly chosen examples.
The paper is organized as follows: In Sec. 2. we describe the method of Newton's identities and give examples when it succeeds and when it fails.In Sec. 3. we present the concept of formal orthogonal polynomials by which the problem is decomposed into a separate search for the distinct zeros and a subsequent search for their multiplicities.Again, we give both kinds of examples.In Sec. 4. we add final ingrediences by which the method of formal orthogonal polynomials becomes accurate.In Sec. 5. we sketch an extension of the method towards generalized local eigenvalue analysis of a matrix pencil and give an example of eigenvalue analysis of the eddy current case of Maxwell's equations.We give conclusions in Sec. 5.
In the paper we use mathematical terminology which might be behind the scope of the journal.We recommend readers interested in a deeper understanding to consult the monographs [9] and [3].

Newton's Identities
Let Ω ⊂ C be a simply connected domain, f : C → C be holomorfic in Ω, and γ be a positively oriented curve such that f is nonvanishing along γ.We consider the problem of locating zeros of f in the interior of γ.
According to Cauchy's residue theorem, we obtain: where z 1 , . . ., z n are the mutually distinct zeros of f inside γ and α 1 , . . ., α n are their respective multiplicities.We denote by N := s 0 the total number of the zeros and construct the polynomial: having the same zeros (including their multiplicities) as function f in the interior of γ.The coefficients of P N (z) are given by following Newton's identities: We obtain s k from Eq. ( 1) by numerical integration.However, the method of Newton's identities is usually ill-conditioned due to bad conditioning of the polynomial, i.e., small changes of σ k generate larger changes in the zeros of Eq. ( 2).Therefore, the contour integrals have to be approximated with high accuracy.

Adaptive Numerical Integration
For the sake of simplicity we assume the circular curve γ with the parametrization z(t) := c+ρe 2πit , t ∈ 0, 1 .We employ the composite trapezoidal rule: where As g k is periodic with period 1, then and This allows us to double n successively until In all examples throuhout the paper we underline the accurate digits.
The poor results of the last example are caused by the instability of the mapping of the polynomial coefficients to the roots.In the next section we introduce the concept of formal orthogonal polynomials that will allow to separate the problem into two subtasks: • first determine all mutually distinct zeros of f , • and then determine their multiplicities.
By this approach the last example gets well-posed.

Formal Orthogonal Polynomials
Let P be the linear space of polynomials with complex coefficients and •, • : P × P → C be the following symmetric bilinear form: Note that Eq. ( 1) and Eq. ( 9) are related, By definition the coefficients of an FOP solve the following linear system with a Hankel matrix: Thus, FOP ϕ t is unique iff H t is nonsingular.In such a case ϕ t is called a regular FOP and the index t is a regular index.

We further introduce Hankel matrix
Theorem 1. [5] • Rank H n+p = n for all p ∈ N ∪ {0}, • for a regular index t ≥ 1 the zeros of FOP ϕ t are the eigenvalues of matrix pencil H t − λH t , t −λH t .We have no information about the remaining t − n eigenvalues.
Theorem 1 suggests to replace the computation of the zeros of ϕ n by determining the eigenvalues of matrix pencil H (1) n − λH n .The multiplicities α 1 , . . ., α n solve the following Vandermonde system, for which an efficient algorithm (relying on Newton polynomial interpolation formula) exists [8], The algorithm starts with computing N := s 0 by numerical integration.Then it continues to compute s 1 , . . ., s 2N −2 .Number n of mutually distinct zeros is equal to rank of H N .Yet, the method has two bottlenecks: • In practice, it may be difficult to determine rank H N , since the difference between the zero singular values and the least nonzero singular value is often small.
• The approximation of the eigenvalues z 1 , ..., z n can be inaccurate, since the matrix pencil H The method of FOP works fine for small numbers of distinct zeros regardless their multiplicity.However, the next example documents that the larger numbers of distinct zeros are troublesome.Their respective multiplicities are as follows:

Accurate FOP Method
The problem with the method of FOP is that matrix pencil H n − λH n is ill-conditioned.Therefore, we shall represent ϕ n in a monic, but generally nonmonomial basis ψ j , This translates Eq. ( 10) to: or equivalently where we introduced Gramm matrix We further define Theorem 2.
• Then λ * is an eigenvalue of H t − λG t .
• Let r be the largest regular index less or equal to t.
Then the eigenvalues of G r − λG r are eigenvalues of G (1) t − λG t .We have no information about the remaining t − r eigenvalues.
The theorem suggests to replace the pencil of Hankel matrices by the related pencil of Gramm matrices.For a regular index t ≥ 1 the zeros t − λG t , where µ = s1 s0 .In particular, the eigenvalues of G (1) When the zeros have positive real parts the condition number of the matrix pencils improves as follows: The multiplicities can be computed as follows: c

ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING
A good choice for the basis turns out to be FOPs ϕ j .From now on let ψ j := ϕ j .In case when H n is strongly regular, meaning that all the leading principal submatrices are regular, then all FOPs ϕ 0 , ϕ 1 , ..., ϕ n are regular, G n is diagonal, and G (1) n is tridiagonal.In the other case, H n is not strongly regular, we establish a set of regular indices {i k }, k = 0, ..., K, where K is the number of regular blocks in H n .If n ≥ 1, then i 0 = 0, i 1 = 1, and i K = n.We define sequence {ϕ t } ∞ t=0 of monic polynomials as follows: If t is a regular index, then ϕ t is the regular FOP.Otherwise, t is not regular, i.e. ϕ t := z t−r ϕ r (z), where r is the largest regular index less than t.In this case ϕ t is called an inner polynomial.The polynomials are grouped into blocks such that every block starts with a regular polynomial and the remaining polynomials in the block are inner polynomials, . . .
for p, q = 0, ..., K − 1, where δ p ∈ C lp×lp is nonsingular and symmetric.The entries of δ p vanish above the main antidiagonal and they are equal to z ip+lp−1 , ϕ ip along the main antidiagonal.
The following theorem states that G n is nonsingular, symmetric and block tridiagonal.
for p, q = 0, ..., K − 1, where δ (1) p ∈ C lp×lp is symmetric and lower anti-Hessenberg.The entries of δ (1) p are equal to z ip+lp−1 , ϕ ip along the first antidiagonal and they vanish above it.Matrix p vanishes up to the entry in the south-east corner, which is equal to For instance, assume n = 10 and the following blocks Then the structure of G n and G n is as follows: G The entries marked are nonzero and in each block they are all equal.
We call a regular FOP to be well-conditioned if the corresponding system Eq.( 11) is well-conditioned, otherwise, the regular FOP is referred to as ill-conditioned.To obtain a numerically stable algorithm, it is crucial to generate only well-conditioned regular FOPs and replace the ill-conditioned regular FOPs by inner polynomials.
The algorithm takes three entries on the input: the bilinear form •, • , ε cond , and ε stop with ε stop < ε cond .The value ε cond determines the length of block l p in Eq. ( 23) and ε stop is a stopping criterion.
If | ϕ r , ϕ r | ≥ ε cond for some regular index r, then ϕ r+1 is generated as an FOP.Otherwise, we search for the smallest t such that t ≤ N −1−r and | z t ϕ r , ϕ r | > ε cond .If we succeed, t + 1 is the length of the block and t is the number of inner polynomials in the block.If we fail to find such t we shall check | z t ϕ r , ϕ r | < ε stop for all t ∈ {0, • • • , N − 1 − r}.In case that the latter condition holds true, then n := r and we shall compute the zeros of ϕ r .Otherwise, if | z t ϕ r , ϕ r | ≥ ε stop for some t, then the length of the block is t := arg max c 2017 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING Algorithm 1 Accurate FOP method [5].
generate ϕ r+1 as a regular FOP  If we set ε cond := 100 and ε stop := 10 −12 , then the algorithm finds that n = 10.Polynomials ϕ 0 , ϕ 1 , ϕ 4 , ϕ 6 , ϕ 8 , ϕ 10 are regular FOPs, while ϕ 2 , ϕ 3 , ϕ 5 , ϕ 7 , ϕ 9 are inner polynomials.We obtain the following approximations of zeros: This approximation of the zeros is less accurate.The higher ε cond gives rise to a higher number of inner polynomials.However, this does not improve the accuracy since the polynomials ϕ 0 , • • • , ϕ 10 are well-conditioned.This approximation of the zeros is more accurate.Here the higher ε cond leads to a better accuracy as we are replacing ill-conditioned regular FOPs by the inner polynomials.We do not observe a remarkable rise of computational time.

Local Eigenvalue Analysis
The theory and algorithm presented in this section can be found in [1] and [4].
Theorem 5. (Weierstrass Canonical Form) Let (zB − A) ∈ C n×n be a regular pencil.Then there exist nonsingular matrices P, Q ∈ C n×n such that in Eq. ( 32), where J i , N i ∈ C ki×ki are Jordan blocks, N i is nilpotent and I ki ∈ C ki×ki is the identity matrix.
According to the structure of the Jordan blocks we divide P and Q into block rows P i ∈ C ki×n and block columns Q i ∈ C n×ki , respectively, for i = 1, 2, • • • , r.By Eq. (32) we obtain Let α i be an eigenvalue of matrix J i .Then and The regular pencil (zB − A) −1 can be decomposed into: Let γ be a positively oriented closed Jordan curve and G interior of γ.Define: and We can extract Jordan blocks whose eigenvalues are contained in G.The following collective notations are constructed: The Jordan blocks J i ; α i ∈ G are collected to form the k γ × k γ Jordan matrix J γ , where k γ = i;αi∈G k i .Similarly, the corresponding Q i and P i are collected to form Q γ ∈ C n×kγ and P γ ∈ C kγ ×n , respectively.
By using the Residue theorem, the matrix in Eq. (37) and Eq. ( 38) can be written as: and Theorem 6.Let V be arbitrary n × k γ matrix.Define a size-reduced moment matrices The number of eigenvalues inside γ is given by The trace of ((zB − A) −1 B) in Eq. ( 41) is approximated by with some integer L 0 , where the elements of the sample vectors v i ∈ R N are taken as −1 or 1 with equal probability.
Require: matrices A, B, curve γ := r+ρ e 2πit , number of integration points N Ensure: eigenvalues of A − λB inside γ 10: compute the eigenvalues of the matrix pencil zM 0 − M 1 Example 7. We consider a matrix pencil (A, B) ∈ R 800×800 , where both matrices are real valued and symmetric positive definite.The matrices arise from the coupled FEM-BEM discretization of the eddy current case of Maxwell's equations, which is a parabolicelliptic system, for simulations of electromagnetic forming of metalic sheets.The matrices determine transient electric field e(t) as follows: A e (t) + B e(t) = b(t), e(0) = e 0 , for which the knowledge of the generalized eigenvalues and eigenvectors is essential.For a more detailed description we refer to Fig. 1.
We search for the generalized eigenvalues inside the curve z(t) := 3000+1500 e 2πit .The radius of the curve is chosen to locate exactly ten eigenvalues.The number of integration points of the trapezoidal rule is 1000, L 0 := 30.
The following eigenvalues are obtained:

Conclusion
We presented contour integral method for localization of roots of analytic functions and its extensions towards local generalized eigenvalue analysis Ae = λBe.
For the latter the linear system with perturbed matrix A − z j B is solved at each integration point z j .Thus, our further research shall address two essential problems: • To find a suitable quadrature method to minimize the number of evaluations.
• To find a reasonable way to update the solver, when perturbing the matrix.The mechanical displacements of the metalic sheet is depicted by the solid lines, while the shape of the form against which the sheet is pressed is depicted with the dashed lines.Note that after the Lorentz forces vanish the deformation of the sheet continues due to the inertial forces.

Fig. 1 :
Fig.1: Snapshots of a simulation of electromagnetic forming of a cylindrical metalic sheet: The electromagnetic field and the mechanical displacements are sketched in the cylindrical coordinates (the axis is on the left hand side) at four times (topdown).Positions of line circular turns of the excitation coil are marked by the blue crosses.The colour map corresponds to the magnitude of the eddy current distribution, the black circles and crosses correspond to the outward and inward orientation, respectively, of the eddy currents.The Lorentz forces are depicted with arrows.The mechanical displacements of the metalic sheet is depicted by the solid lines, while the shape of the form against which the sheet is pressed is depicted with the dashed lines.Note that after the Lorentz forces vanish the deformation of the sheet continues due to the inertial forces.