*Numerical Experiments*

The steps outlined in Section 2.1.1 are readily modified using the sequences in (5) to yield spectral enclosures for matrix polynomials.

**Application 5** ([31], Example 3)**.** *We consider the* 50 × 50 *matrix polynomial P*(*λ*) = *<sup>A</sup>*2*λ*<sup>2</sup> + *A*1*λ* + *A*0*, where*

$$A\_2 = l\_{50\prime} \cdot A\_1 = \text{triangle} \{-3, 9, -3\}, \ A\_0 = \text{triangle} \{-5, 15, -5\},$$

*describing a damped mass-spring system [14,32] and set non-negative weights* **w** = {1, 1, <sup>1</sup>}*, measuring perturbations of the coefficient matrices Aj*<sup>2</sup>*j*=<sup>0</sup> *in an absolute sense. We initiate the procedure with 15 equidistributed initial points* #*μj*0\$<sup>15</sup>*j*=<sup>1</sup>*on the semicircle*

$$\left\{ z \in \mathbb{C} : |z| = 15 \left(= \text{median}\_{j=0,1,2} \left( \left\| A\_j \right\|\_1 \right) \right), \text{Im}(\mathbf{z}) \ge 0 \right\} $$

*and proceed to determine eigenvalue approximating sequences* #*μjk*\$*nj k*=0 *(j* = 1, ... 15*) according to (5), so that their final terms all lie in the interior of <sup>σ</sup>*0.01,**w**(*P*)*. This computation requires 722 iterations. As in the constant matrix case, interpolation between the values of jk such that μjk* ∈ *∂σjk*,**<sup>w</sup>**(*P*) *along the trajectories formed by* #*μjk*\$*nj k*=0 *(j* = 1, ... 15*) results in approximations of ∂σ*,**<sup>w</sup>**(*P*) *for* = 0.1, 0.2, 0.3, 0.4, 0.5*, as seen in Figure 5a. Note this yields a sufficiently accurate sketch of σ*,**<sup>w</sup>**(*P*) *and is very competitive when compared to other methods. For instance, Figure 5b is obtained via the procedure in [31] applied to a* 400 × 400 *grid on the relevant region* Ω = [−20, 5] × [−15, <sup>15</sup>]*. This latter approach is far more computationally intensive, requiring 71,575 iterations to visualize σ*,**<sup>w</sup>**(*P*) *for the same tuple of parameters.*

In case a more detailed spectral localization is desired, our method may be adapted, as in Application 1, to identify individual pseudospectrum components. Our next experiment also serves to illustrate the fact that the number of initial trajectories that are attracted by the individual eigenvalues to form the related clusters is intimately connected to eigenvalue sensitivity.

**Application 6** ([13], Example 5.1)**.** *We consider the mass-spring system from ([13], Ex. 5.2) defining the* 3 × 3 *selfadjoint matrix polynomial*

$$P(\lambda) = A\_2 \lambda^2 + A\_1 \lambda + A\_0 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 5 \end{bmatrix} \lambda^2 + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 3 & -1 \\ 0 & -1 & 6 \end{bmatrix} \lambda + \begin{bmatrix} 2 & -1 & 0 \\ -1 & 3 & 0 \\ 0 & 0 & 10 \end{bmatrix}$$

*and set* **w** = {1, 1, <sup>1</sup>}*. As in Application 5, computations are restricted exclusively to the closed upper half-plane. However, the close proximity of the eigenvalues* −0.08 ± **i**1.45, −0.75 ± **i**0.86*,* −0.51 ± **i**1.25 *(indicated by "\*" in Figure 6), as well as the fact that the pair λ* = −0.51 ± **i**1.25 *is less sensitive than the other two, necessitates the use of many initial points. Indeed, as demonstrated in Figure 6a, initiating the procedure with 40 equidistributed initial points on*

#*z* ∈ C : |*z*| = min*Aj*1, Im(z) ≥ 0\$ *results in σ*(*P*) *being under-represented in the resulting clusters. In order to correctly approximate all three elements of the spectrum on the upper half-plane enforces the use of as many as 80 points on the selected semicircle. The length nj of each sequence* #*μjk*\$*nj k*=0 (*j* = 1, ... 80) *is determined, so that all values* #*<sup>s</sup>*min(*P*(*μjnj*))\$<sup>80</sup>*j*=<sup>1</sup> *do not exceed the prefixed parameter value of* 0 = 0.01*; this construction involved 1162 singular value computations. Using the squared Euclidean distance as the metric for computing the cluster evaluation criterion, three distinct groups are correctly identified, each converging to a different eigenvalue in the closed upper half-plane, as in Figure 6b. Note that the least sensitive eigenvalue λ* = −0.51 + **i**1.25 *ends up attracting only one of these sequences; the corresponding group being a singleton. To correctly sketch the boundaries of σ*,**<sup>w</sup>**(*P*) *for the triple of* = 0.24*,* 0.48*,* 0.73 *(>*0 = 0.01)*, we introduce six additional points for each cluster. Indeed, denoting c*1*, c*2*, c*3 *the centroids of the clusters, for each cluster j* = 1*,* 2*,* 3 *we consider the vertices* #*pji*\$<sup>6</sup>*i*=<sup>1</sup> *of a canonical hexagon centered at cj with maximal diameter equal to* min""*cj* − *ci*""*<sup>i</sup>*=*j. These vertices are indicated by circles in Figure 6b and are used as starting points to construct the additional trajectories indicated by the black lines in Figure 6c. Note that all three selected parameters = 0.24, 0.48, 0.73 should be possible to interpolate along these additional lines as well, which explains why most of these trajectories have been extended to the opposite directions as well, modifying the definition of the sequences in (5) in each instance accordingly. The construction of the additional sequences requires 202 singular value computations (leading to a total of 1364 iterations), while the resulting approximations of pseudospectra boundaries for* = 0.24, 0.48, 0.73 *are depicted in Figure 6c.*

**Figure 5.** Pseudospectrum computation for a damped mass-spring system. (**a**) Approximate pseudospectra visualization, interpolating along 15 trajectories of converging sequences. (**b**) Pseudospectra visualization, using the modified grid algorithm in [31].

(**a**) Partial spectral identification, due to close eigenvalue proximity.

(**b**) Complete spectral identification with increased number of initial points.

(**c**) Pseudospectra visualizations for = 0.24, 0.48, 0.73.

**Figure 6.** Pseudospectra computations for a vibrating system.

**Application 7** ([13], Example 5.3)**.** *This experiment tests the behavior of the method on a damped gyroscopic system described by the* 100 × 100 *matrix polynomial*

$$P(\lambda) = M\lambda^2 + (G+D)\lambda + K\_r$$

*with*

$$\begin{aligned} M &= I\_{10} \otimes \frac{4I\_{10} + B + B^T}{6} + 1.30 \frac{4I\_{10} + B + B^T}{6} \otimes I\_{10\prime} \\ G &= 1.35 I\_{10} \otimes (B - B^T) + 1.10 (B - B^T) \otimes I\_{10\prime} \\ D &= \text{triddiag}\{-0.1, 0.3, -0.1\}, \\ K &= I\_{10} \otimes (B + B^T - 2I\_{10}) + 1.20 (B + B^T - 2I\_{10}) \otimes I\_{10} \end{aligned}$$

*and B the* 10 × 10 *nilpotent matrix having ones on its subdiagonal and zeros elsewhere. Note M*, *K are positive and negative definite respectively, G is skew-symmetric, and the tridiagonal D is a damping matrix.*

*Starting with 50 points on*

$$\{z \in \mathbb{C} : |z| = 15 \\ \text{= } \text{median}(\|K\|\_{1'} \|G + D\|\_{1'} \|M\|\_{1}), \text{Im}(\mathbf{z}) \ge 0\}$$

*and then 5 additional points on the perpendicular bisector of the line segment defined by the two centroids of the resulting clusters (indicated by the blue circles), the resulting pseudospectrum approximation required 1212 iterations in total with* 0 = 0.002 *and can be seen in the left part of Figure 7a. The algorithm in [31] applied to a* 300 × 300 *grid on the region* Ω = [−4, 4] × [−3, 3] *required 29,110 iterations for pseudospectra visualization for the same triple* = 0.004, 0.02, 0.1 *to obtain comparable results in Figure 7b.*

**Figure 7.** Comparison of pseudospectra computation for a damped gyroscopic system and = 0.004, 0.02, 0.1. (**a**) Computation using 50 initial points. (**b**) Computation using algorithm in [31].

We conclude this section, examining the behavior of the method on a non-symmetric example. **Application 8** ([31], Example 2)**.** *We consider the* 20 × 20 *gyroscopic system*

$$P(\lambda) = A\_2 \lambda^2 + A\_1 \lambda + A\_0 = I\_{20} \lambda^2 + \mathbf{i} \begin{bmatrix} I\_{10} & 0 \\ 0 & 5I\_{10} \end{bmatrix} \lambda + \begin{bmatrix} 1 & -1 & -1 & \cdots & -1 \\ -1 & 1 & -1 & \cdots & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1 & -1 & -1 & \cdots & 1 \end{bmatrix}$$

*and* **w** = {1, 1, <sup>1</sup>}*. Starting with 21 points on*

$$\{z \in \mathbb{C} : |z| = 25 \\ = \|A\_0\|\_1 + \|A\_1\|\_1), \text{Re}(z) \le 0\}$$

*and* 0 = 0.001*, five clusters are detected (Figure 8a) after 1140 iterations. Then, two additional points are introduced on each of the line segments defined by the centroids of the detected clusters (indicated by the blue circles in Figure 8b), causing the iterations to rise to the total number of 2662 in order to determine the 20 corresponding trajectories (indicated by grey lines in Figure 8c). The corresponding visualizations in Figure 8d), obtained via [31] applied to a* 400 × 400 *grid on the region* Ω = [ −20, 20] × [−15, 10] *required 88,462 iterations.*

**Figure 8.** Comparison of pseudospectra computation for a gyroscopic system and = 0.1, 0.2, 0.4, 0.6, 0.8. (**a**) Cluster detection using 21 initial points. (**b**) Locations of additional points. (**c**) Pseudospectra visualizations interpolating along the trajectories of 21 initial and 20 additional points. (**d**) Pseudospectra visualization, using the modified grid algorithm in [31].

#### **4. Concluding Remarks**

In this note, we have shown how to define sequences which, beginning from arbitrary complex scalars, converge to some element of the spectrum of a matrix. This approach can be applied both to constant matrices and to more general matrix polynomials and can be used as a means to obtain estimates for those eigenvalues that lie in the vicinity of the initial term of the sequence. This construction is also useful when no information on the location of the spectrum is a priori known. In such cases, repeating the construction from arbitrary points, it is possible to detect peripheral eigenvalues, localize the spectrum and even obtain spectral enclosures.

As an application, in this paper we used this construction to compute the pseudospectrum of a matrix or a matrix polynomial. Thus, a useful technique for speeding up

pseudospectra computations emerges, which is essential for applications. An advantage of the proposed approach is that does not make any assumptions on the location of the spectrum. The fact that all computations are performed on some dynamically chosen locations on the complex plane which converge to the eigenvalues, rather than on a large number of predefined points on a rigid grid, can be seen as improvement over conventional grid algorithms.

Parallel implementation of the method can lead to further computational savings when applied to large matrices. Another option would be to apply this method combined with randomization techniques for the selection of the initial points of the sequences. In the large-scale matrix case, this method may be helpful in obtaining a first impression of the shape and size of pseudospectrum and even computing a rough approximation. Then, if desired, this could be used in conjunction with local versions of the grid algorithm and small local meshes about individual areas of interest.

**Author Contributions:** All authors have equally contributed to the conceptualization of this paper, to software implementation and to the original draft preparation. Funding acquisition and project administration: P.P. All authors have read and agreed to the submitted version of the manuscript.

**Funding:** This research is carried out/funded in the context of the project "Approximation algorithms and randomized methods for large-scale problems in computational linear algebra" (MIS 5049095) under the call for proposals "*Researchers' support with an emphasis on young researchers–2nd Cycle'." The project is co-financed by Greece and the European Union (European Social Fund—ESF) by the Operational Programme Human Resources Development, Education and Lifelong Learning 2014–2020.*

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
