## 3 The Contour Integral Method for Solving NEP
### Practical Applications  
In what follows, some discussions about the contour integral method are presented from the perspective of practical applications, which involve the numerical implementation of the method itself and computation of dominant eigenvalues of controlled structures with multiple time-delays. First, one typically does not know a priori number \(m\) of the simple eigenvalues inside the prescribed contour \(\partial C\). However, as long as the order number \(q\) of moments in (18) is chosen large enough, such that \(\mathbf{\Sigma}_{0}\) of the SVD of \(\mathbf{B}_{0}\) in (23) contains at least one zero singular value. One can detect \(m\) by the rank test of the matrix \(\mathbf{B}_{0}\), for example, by counting the number of singular values of \(\mathbf{B}_{0}\) that are above a user-defined tolerance, and obtain a truncated SVD of \(\mathbf{B}_{0}\) with rank \(m\) as in (23). Alternatively, counting the eigenvalues can also be achieved by the argument principle,  
\[(2\pi\mathrm{i})^{-1}\oint_{\partial C}\frac{g^{\prime}(z)}{g(z)}\mathrm{d}z= N_{z}-N_{p} \tag{27}\]  
where \(g(z)=\det(\mathbf{T}(z))\), \(N_{z}\) and \(N_{p}\) denote the number of roots (i.e., eigenvalues of \(\mathbf{T}(z)\)) and poles of \(g(z)\) within the contour \(\partial C\) in the complex plane, respectively. By using the following trace relation  
\[\frac{g^{\prime}(z)}{g(z)}=\mathrm{trace}\Big{(}\mathbf{T}^{-1}(z)\mathbf{T}^ {\prime}(z)\Big{)} \tag{28}\]  
The expression (27) can be rewritten as  
\[(2\pi\mathrm{i})^{-1}\oint_{\partial C}\mathrm{trace}\Big{(}\mathbf{T}^{-1}(z )\mathbf{T}^{\prime}(z)\Big{)}\mathrm{d}z=N_{z}-N_{p} \tag{29}\]  
Ordinarily, requiring prior knowledge of \(N_{p}\) in (27) or (29) hinders one from determining \(N_{z}\). Fortunately, for controlled systems involved in this paper the number \(N_{p}\) of all poles of its reduced characteristic function \(f_{\mathrm{m}}(z)=\det(\mathbf{J}_{\mathrm{m}}(z))\) inside a given contour \(\partial C\) is identical to that of all poles of its open-loop receptances \(\mathbf{H}_{0m}(\lambda)\) in (12) or (13) within the same contour (based on the assumption that these poles are simple), and the latter can be determined experimentally. Second, the numerical evaluation of the contour integrals for \(\mathbf{A}_{k}\) in (18) is usually performed by the trapezoid sum for which \(n_{\mathrm{c}}\) quadrature points are used. It is preferred to choose the contour \(\partial C\) to be a circle with \(\partial C=\alpha+\beta\varepsilon^{i\theta}\), where \(\alpha\) and \(\beta\) are the center and radius of \(\partial C\), respectively, because the quadrature error in this case yields the exponential decay with an exponent that depends on the product of the number of quadrature nodes and the smallest distance of the eigenvalues to the contour [22]. To improve numerical conditioning, one also needs to shift and scale the weight function \(z^{k}\) in (18). Thus, a discretized version of (18) is obtained,  
\[\mathbf{A}_{k}\approx\frac{1}{n_{\mathrm{c}}}\sum_{j=0}^{n_{\mathrm{c}}-1} \bigg{(}\frac{\omega_{j}-\alpha}{\beta}\bigg{)}^{k+1}\mathbf{T}^{-1}\Big{(} \omega_{j}\Big{)}\mathbf{R},\ k=0,1,\ldots,2q-1. \tag{30}\]  
where \(\omega_{j}=\alpha+\beta\varepsilon^{2\pi\mathrm{i}(j+1/2)/n_{\mathrm{c}}}\) for \(j=0\), \(1\),..., \(n_{\mathrm{c}}-1\). Additionally, the integrand \(\mathbf{T}^{-1}\Big{(}\omega_{j}\Big{)}\mathbf{R}\) is the solution \(\mathbf{Y}\Big{(}\omega_{j}\Big{)}\) of the linear system  
\[\mathbf{T}\Big{(}\omega_{j}\Big{)}\mathbf{Y}\Big{(}\omega_{j}\Big{)}=\mathbf{R} \tag{31}\]  
evaluated at each quadrature point \(\omega_{j}\). Third, the second-order time delayed systems discussed in this paper are TDSs of retarded type.