Approximation of PDE eigenvalue problems involving parameter dependent matrices

We discuss the solution of eigenvalue problems associated with partial differential equations that can be written in the generalized form $\mathsf{A}x=\lambda\mathsf{B}x$, where the matrices $\mathsf{A}$ and/or $\mathsf{B}$ may depend on a scalar parameter. We consider in particular the approximation of Poisson eigenvalue problem using the Virtual Element Method (VEM) and show that the presence of one (or both) parameters can produce unexpected results.


Introduction
Several schemes for the approximation of eigenvalue problems arising from partial differential equations lead to the algebraic form: find λ ∈ R and x ∈ R n with x = 0 such that where A and B are matrices in R n×n . We consider the case when the matrices A and B are symmetric and positive semidefinite and may depend on a parameter. This is a typical situation found in applications where elliptic partial differential equations are approximated by schemes that require suitable parameters to be tuned (for consistency and/or stability reasons). In this paper we discuss in particular applications arising from the use of the Virtual Element Method (VEM), see [9,4,7,10,11,6,12], where suitable parameters have to be chosen for the correct approximation.
In general, it may be not immediate to describe how the matrices A and B depend on the given parameters. For simplicity, we consider the case when the dependence is linear: under suitable assumptions it is easy to discuss how the computed spectrum varies with respect to the parameters.
We assume that the matrices A and B satisfy the following condition for C = A, B.
Assumption 1. The matrix C can be split into the sum where γ is a non negative real number and C 1 and C 2 are symmetric. The matrices C 1 and C 2 satisfy the following properties: a) C 1 is positive semidefinite with kernel K C1 ; b) C 2 is positive semidefinite and positive definite on K C1 ; c) C 2 vanishes on K ⊥ C1 , the orthogonal complement of K C1 in R n .
In Section 2 we describe the spectrum of (1) as a function of the parameters, in various situations that mimic the behavior of matrices A and B originating from several discretization schemes.
Section 3, which is the core of this paper, discusses the influence of the parameters on the VEM approximation of eigenvalue problems. Several numerical examples complete the papers, showing that the parameters have to be carefully tuned and that wrong choices can produce useless results.

Parametric algebraic eigenvalue problem
Given two symmetric and positive semidefinite matrices A and B that can be written as (3) A = A 1 + αA 2 and (4) with nonnegative parameters α and β, we consider the eigensolutions to the generalized problem (1). We assume that the splitting of the matrices A and B is obtained with symmetric matrices and satisfies Assumption 1 for C 1 = A 1 , B 1 and C 2 = A 2 , B 2 . Moreover we denote by n A1 and n B1 the dimension of K A1 and K B1 , respectively. Remark 1. Problem (1) has n eigenvalues if and only if rankB = n, see [8]. If B is singular the spectrum can be finite, empty, or infinite (if A is singular too). If A is non singular, usually one can circumvent this difficulty by computing the eigenvalues of Bx = µAx and setting λ = 1/µ. The kernel of B is the eigenspace associated with the vanishing eigenvalue with multiplicity m, and the original problem has exactly m eigenvalues conventionally set to ∞.
We want to study the behavior of the eigenvalues as the parameters α and β vary. We consider three cases.
2.1. Case 1. We fix β > 0 so that B is positive definite. This implies that the eigenvalues of (1) are all non negative. Let us consider first α = 0 so that (1) reduces to (5) A 1 x = λBx.
Since A 1 is positive semidefinite, λ = 0 is an eigenvalue of (5) with multiplicity equal to n A1 = dim(K A1 ) and K A1 is the associated eigenspace. In addition, we have m A = n−n A1 positive eigenvalues {µ 1 ≤ · · · ≤ µ m A } counted with their multiplicity (since we are dealing with a symmetric problem, we do not distinguish between geometric and algebraic multiplicity). We denote by v j ∈ K ⊥ A1 the eigenvector associated with µ j , that is A 1 v j = µ j Bv j . Thanks to property c) of Assumption 1 when C = A, we observe that Therefore (µ j , v j ), for j = 1, . . . , m A , are eigensolutions of the original system (1).
On the other hand, the eigensolutions of are characterized by the fact that n A1 eigenvalues ν i (i = 1, . . . , n A1 ) are strictly positive with corresponding eigenvectors w i belonging to K A1 , while the remaining m A eigenvalues vanish and have K ⊥ A1 as eigenspace. Thus, property a) of Assumption 1, for C = A, yields which means that (αν i , w i ), for i = 1, . . . , n A1 , are eigensolutions of (1).
Summarizing the eigenvalues of (1) are: The left panel in Figure 1 shows the eigenvalues of a simple example where A ∈ R 6×6 is obtained by the combination of diagonal matrices with entries (7) diag(A 1 ) = [3, 4, 5, 6, 0, 0], diag(A 2 ) = [0, 0, 0, 0, 1, 2]. and B = I 6 is the identity matrix. Along the vertical lines we see the eigenvalues corresponding to a fixed value of α. The eigenvalues 3, 4, 5, 6 are associated with eigenvectors in K ⊥ A1 and do not depend on α. The solid lines starting at the origin display the eigenvalues 1, 2 multiplied by α.
Remark 2. We observe that if A 2 is not positive definite on K A1 , its kernel has a nonempty intersection with K A1 . Let n 12 be the dimension of this intersection, then problem (1) admits n 12 vanishing eigenvalues which appear in the first case of (6).

Case 2.
Let us now fix α > 0, so that A is positive definite. We have that all the eigenvalues are positive. We observe that when β = 0, the matrix B = B 1 may be singular, therefore it is convenient to consider the following problem: where χ = 1 λ . If χ = 0, we conventionally set λ = ∞. Problem (8) reproduces the same situation we had in Case 1, with the matrices A and B switched. Repeating the same arguments as before, we obtain that problem (8) has two families of eigenvalues . Going back to the original problem (1), we can conclude that the eigensolutions of (1) are the following ones: , r k−n B 1 for k = n B1 + 1, . . . , n.
In the right panel of Figure 1, we report the eigenvalues of a simple example where A = I 6 and B is obtained by combining B 1 = A 1 and B 2 = A 2 defined in (7). We see that the eigenvalues 1 3 , 1 4 , 1 5 , 1 6 are independent of β and that the remaining two eigenvalues lie along the hyperbolas 1 β and 1 2β , plotted with solid line. 2.3. Case 3. We consider now the case when α and β can vary independently from each other. We have different situations corresponding to the relation between K A1 and K B1 . To ease the reading, let us introduce the following notation: In this case the space R n can be decomposed into four mutually orthogonal subspaces ). Let us denote by n A1∩B1 the dimension of K A1 ∩ K B1 . If K A1 ∩ K B1 = ∅, for x ∈ K A1 ∩ K B1 the eigenproblem to be solved is αA 2 x = λβB 2 x, hence the eigenvalues are given by α β η i i = 1, . . . , n A1∩B1 , see (10d). Next, if x ∈ K A1 ∩ K ⊥ B1 we have to solve αA 2 x = λB 1 x, which admits (αχ i , y i ) i = 1, . . . , n A1 −n A1∩B1 as eigensolutions where (χ i , y i ) are defined in (10c). Similarly, if x ∈ K ⊥ A1 ∩ K B1 , we find that the eigensolutions are 1 β ν i , w, i i = 1, . . . , n B1 − n A1∩B1 with (ν i , w i ) given by (10b). In the last case, x ∈ K ⊥ A1 ∩K ⊥ B1 , the matrices A and B are non singular and thanks to property c) in Assumption 1, for C = A and C = B, we obtain that the eigenvalues are positive and independent of α and β and correspond to those of (10a). In conclusion, we have We report in Figure 2 the eigenvalues illustrating this last case when we have diagonal matrices given by diag( The surface contains the eigenvalues depending on both α and β, the hyperbolas those depending only on β and the straight lines those depending only on α. If we cut the three dimensional picture with a plane at β > 0 fixed we recognize the behavior analyzed in Section 2.1 and shown in Figure 1 left. Analogously, taking a plane with α > 0 fixed, we recover Case 2 (see Section 2.2). If K A1 ∩ K B1 = ∅, we set n A1∩B1 = 0, hence the eigenvalues are In order to illustrate the case K A1 ∩K B1 = ∅, we report in Figure 3 the eigenvalues computed using the following diagonal matrices with entries For a fixed α, we can see in solid line the hyperbolas νj β , j = 1, 2 while when β is fixed we can see the straight lines αχ j , j = 1, 2. The remaining two eigenvalues are independent of α and β.

Virtual element method for eigenvalue problems
In this section we recall how algebraic eigenvalue problems similar to the ones discussed in the previous section can be obtained withing the framework of the Virtual Element Method (VEM) for the discretization of elliptic eigenvalue problems, see [7,6]. in Ω u = 0 on ∂Ω.
In view of the application of VEM, we consider the weak form: find λ ∈ R and u ∈ H 1 0 (Ω) with u = 0 such that , and (·, ·) is the scalar product in L 2 (Ω).
It is well-known that problem (11) admits an infinite sequence of positive eigenvalues 0 < λ 1 ≤ · · · ≤ λ i ≤ · · · repeated according to their multiplicity, each one associated with an eiegenfunction u i with the following properties Let us briefly recall the definition of the virtual element spaces and of the discrete bilinear forms which we are going to use in this section, see [2,1]. We present only the two dimensional spaces, the three dimensional ones are obtained using the 2D virtual elements on the faces. We decompose Ω into polygons P , with diameter h P and area |P |. Similarly, if e is an edge of an element P , we denote by h e = |e| its length. Depending on the context ∂P refers to either the boundary of P or the set of the edges of P . The notation T h and E h stands for the set of the elements and the edges, respectively. As usual, h = max P ∈T h h P . We assume the following mesh regularity condition (see [2]): there exists a positive constant γ, independent of h, such that each element P ∈ T h is star-shaped with respect to a ball of radius greater than γh P ; moreover, for every element P and for every edge e ⊂ ∂P , it holds h e ≥ γh P .
For k ≥ 1 and P ∈ T h we definẽ , v| e ∈ P k (e) ∀e ⊂ ∂P, ∆v ∈ P k (P )}. We consider the following linear forms on the spaceṼ k h (P ) D1 : the values v(V i ) at the vertices V i of P , D2 : the scaled edge moments up to order k − 2 1 |e| e vm ds ∀m ∈ M k−2 (e), ∀e ⊂ ∂P, D3 : the scaled element moments up to order k − 2 with x ω the barycenter of ω, and with the convention that M −1 (ω) = ∅.
From the values of the linear operators D1-D3, on each element P we can compute a projection operator Π ∇ k :Ṽ k h (P ) → P k (P ) defined as the unique solution of the following problem: where a P (u, v) = (∇u, ∇v) P and (·, ·) P denotes the L 2 (P )-scalar product. The local virtual space is defined as where (P k \ P k−2 )(P ) contains the polynomials in P k (P ) L 2 -orthogonal to P k−2 (P ). We recall that by construction P k (P ) ⊂ V k h (P ), so that the optimal rate of convergence is ensured. Moreover, the linear operators D1-D3 provide a unisolvent set of degrees of freedom (DoFs) for V k h (P ), which allows us to define and compute The global virtual space is In order to discretize problem (11), we introduce the discrete counterparts a h and b h of the bilinear forms a and b, respectively. Both discrete forms are obtained as sum of the following local contributions: where b P (u, v) = (u, v) P , and S P a and S P b are symmetric positive definite bilinear The virtual element counterpart of (11) reads: find λ h and u h ∈ V k h with u h = 0 such that and the corresponding eigenfunctions u ih , for i = 1, . . . , N h , enjoy the discrete counterpart of properties in (12).
The following convergence result has been proved in [7].
Theorem 1. Let λ be an eigenvalue of (11) of multiplicity m and E λ the corresponding eigenspace. Then there are exactly m discrete eigenvalues of (18) λ j(i)h (i = 1, . . . , m) tending to λ. Moreover, assuming that u ∈ H 1+r (Ω), for all u ∈ E λ , the following inequalities hold true: where t = min(k, r),δ(E, F) represents the gap between the spaces E and F, and E h is the eigenspace spanned by u h .
Remark 3. It is also possible to consider on the right hand side of (18) the bilinear . This leads to the following discrete eigenvalue problem: The analogue of Theorem 1 holds true for this partially non stabilized discretization as well.
3.1. Computational aspects and numerical results. In order to compute the solution of problems (18) and (19), we need to describe how to obtain the matrices associated to our bilinear forms. By construction the matrix A 1 (respectively, We observe that the local contributions of the bilinear forms displayed in (16) mimic the following exact relations Let us denote by A 1 , A 2 , B 1 and B 2 the matrices whose entries are given by . Even if the global matrices A and B do not satisfy the properties stated in Assumption 1 it turns out that Assumption 1 is fulfilled by C = B 1 +βB 2 ; moreover, C = A 1 + αA 2 is characterized by the situation described in Remark 2.
We start with the pair A 1 and A 2 . The kernel K A 1 , with abuse of notation, is characterized by , and thus the first term vanishes, while for the second term it is enough to observe that Π ∇ k (Π ∇ k w) = Π ∇ k w. Thus property c) of Assumption 1 is verified for C = A.
Concerning property b) of Assumption 1, we have by construction, that a P (( (20). On the other hand, if v is constant on P , then Π ∇ k v = v is constant, therefore v ∈ K A 1 and (I − Π ∇ k )v = 0 so that v belongs also to the kernel of A 2 . Hence the pair A 1 and A 2 does not satisfy property b), but it is in the situation described in Remark 2.
Let us now consider the pair B 1 and B 2 . We observe that the kernel of B 1 is characterized by Π 0 k v = 0. The analysis performed for the pair A 1 and A 2 can be repeated and gives that in this case Assumption 1 is verified for C = B.
As a consequence of the assembling of the local matrices, the global matrices A 1 and A 2 (B 1 and B 2 , respectively) do not satisfy anymore the properties listed in Assumption 1. In particular, for k = 1 we shall see that the matrices A 1 and B 1 are not singular. Nevertheless, we are going to show that the numerical results look pretty much similar to the ones reported in Section 2.
Moreover, in practice the matrices A 2 and B 2 are not available and they are replaced by using the local bilinear forms S P a and S P b given in (16) as follows. Let us denote by u h , v h ∈ R N P the vectors containing the values of the N P local DoFs associated to u h , v h ∈ V k h (P ). Then, we define the local stabilized forms as the stability parameters σ P and τ P are positive constants which might depend on P but are independent of h. We point out that this choice implies the stability requirements in (17). In the applications, the parameter σ P is usually chosen depending on the mean value of the eigenvalues of the matrix stemming from the term a P (Π ∇ k ·, Π ∇ k ·), and τ P as the mean value of the eigenvalues of the matrix resulting from 1 The choice of the stabilized form S P a is discussed in some papers concerning the source problem, see, e.g., [3] and the references therein. One can find an analysis of the stabilization parameters σ P in [5].
If σ P and τ P vary in a small range, it is reasonable to take σ P = α and τ P = β for all P and this is the situation which we discuss further. Therefore, the structure of the matrices is A = A 1 +αA 2 and B = B 1 +βB 2 where A 2 and B 2 are the matrices with local contribution given by u h v h and h 2 P u h v h , respectively. We study the behavior of the eigenvalues as α and β vary in given ranges.
In the following tests Ω is the unit square partitioned using a sequence of Voronoi meshes with a given number of elements. In Figure 4 we report the coarsest mesh with 50 elements (h = 0.2350, 151 edges, 102 vertices). We recall that the exact eigenvalues are given by (i 2 + j 2 )π 2 for i, j ∈ N \ {0} with eigenfunctions sin(iπx) sin(jπy). The following numerical results have been obtained using Matlab and, in particular, the routine eig for the computation of the eigenvalues. In the following figures, we shall always report the computed eigenvalues divided by π 2 .  Table 1. Dimension of K A1 with respect to k and the number of elements we see that for k = 1 the matrix A 1 is nonsingular. We have computed the lowest eigenvalue of A 1 x = λB 1 x, which gives an estimate of the inf-sup constant of the discrete problem (18). The results, presented in Table 3, show that the first eigenvalue is decreasing, and this behavior corresponds to the fact that the bilinear form P a P (Π ∇ k ·, Π ∇ k ·) is not stable. We now discuss some tests, where we present the behavior of the eigenvalues as the parameters α and β vary, for the mesh with N = 200 and different degree k of the polynomials in the space V k h . The rows of Figure 5 contain the results for fixed k and the values β = 0, 1, 5, while, in the columns, β is fixed and k varies. In each picture, we plot in red the exact eigenvalues and with different colors those corresponding to α = 10 r with r = −3, . . . , 1.
These plots clearly confirm that the choice of the parameters for optimal performance is not so immediate. Consider, in particular, that we are solving the Laplace eigenvalue problem (isotropic diffusion) on a domain as simple as a square. For an arbitrary elliptic problem and more general domains the situation could be much more complicated. For β = 0, the first 30 eigenvalues are well approximated with higher degree of polynomials whenever α ≥ 0.1. The value α = 0.1 seems to be the best choice in the case k = 1. Increasing β does not produce much improvement. All the pictures seem to indicate that higher values of α might give better results. In particular, for k = 2, 3 the first 30 eigenvalues are approximated with a reasonable accuracy for α = 10 and β = 1. Increasing β and keeping α = 10, we see that a smaller number of eigenvalues are captured. Figure 6 shows the behavior of the eigenvalues as α varies from 0 to 10. At a first glance the pictures remind of Figure 1 (left) even if, as it has been explained before, the situation is not exactly matching what we discussed in Section 2.
Each subplot reports all computed eigenvalues between 0 and 40; the dotted horizontal lines represent the exact solutions. The first 30 computed eigenvalues are connected together with lines of different colors in an automated way. An "ideal" good approximation would correspond to a series of colored lines matching the dotted lines of the exact eigenvalues. It is interesting to look at the differences between various degrees (k from 1 to 3 moving from the top to the bottom) and values of β (equal to 0, 1, and 5 from left to right). 1.92654e+00 1.74193e+00 1.06691e+00 6.81927e-01 5.54346e-01 Table 3. First eigenvalues of A 1 x = λB 1 x for different meshes More reliable results seem to be obtained for large k and small β. Actually, the limit case of β = 0 appears to be the safest choice. This is in agreement with the claim of [4] where the authors remark that "even the value σ E = 0 yields very accurate results, in spite of the fact that for such a value of the parameter the stability estimate and hence most of the proofs of the theoretical results do not hold" (note that σ E = 0 in [4] has the same meaning as β in our paper). It is interesting to observe that the analysis of [7], summarized in Theorem 1, covers the case β = 0 as well. On the other hand β = 0 may produce a singular matrix B and this could be not convenient from the computational point of view.
In order to better understand the behavior of the eigenvalues reported in Figure 6(h), we highlight in Figure 7 four eigenvalues that are apparently aligned along an oblique line. The corresponding eigenfunctions are reported in Figure 8. The four eigenfunctions look similar, so that the analogy with Figure 1 (left) is even more evident. We conclude this discussion with an example where, for a given value of α, a good eigenvalue (i.e., an eigenvalue corresponding to a correct approximation) is crossing a spurious one (i.e., an eigenvalue belonging to an oblique line). In this case it may happen that the two eigenfunctions mix together, thus yielding to an even more complicated situation. This behavior is reported in Figure 9, where a region of the plot shown in Figure 6(h) is blown-up close to an intersection point: actually three eigenvalues (a spurious one and two corresponding to good ones) are clustered at the marked intersection points. Figure 10 shows the computed eigenvalues smaller that 40 when β varies from 0 to 5 and for a fixed value of α. As in Figure 5 and in analogy with Figure 6, the rows correspond to the degree k of polynomials, while the columns refer to different values of α. The dotted horizontal lines represent the exact eigenvalues. The lines with different colors in each picture follow the n-th eigenvalue for n = 1, . . . , 30. It turns out that all lines are originating from curves that look like hyperbolas when β is large. Following each of these hyperbolas from β = +∞ backwards, it happens that when the hyperbola meets a correct approximation of an eigenvalue of the continuous problem, it deviates from its trajectory and becomes a (almost   Figure 7 horizontal) straight line. In the case k = 1, we see that the higher eigenvalues are computed with decreasing accuracy as β approaches 0.
We recognize in these pictures the situation presented in Subsection 2.2, corresponding to the behavior of the eigenvalues when the parameter β in matrix B Figure 9. Intersections of good and spurious eigenvalues varies. In this test, the kernel of matrix B 1 is not empty only for k = 3. Nevertheless, we can see that when β approaches 0, there are several eigenvalues going to ∞. On the other side, for greater values of β we obtain several spurious eigenvalues. The range of β, which gives eigenvalues close to the exact ones, clearly depends on k and α. Figure 11 displays, in separate pictures, the first four eigenvalues, with k = 1, α = 10, different values of h, and 0 ≤ β ≤ 400. Taking into account that the routine eig sorts the eigenvalues in ascending order, the four pictures display, in lexicographical order, the first, second, third and fourth computed eigenvalues. In each subplot, each line refers to a particular mesh. We can see that the eigenvalues computed with the finest mesh seem to be insensitive with respect to the value of β. On the opposite side the coarsest mesh gives approximations of the correct values only when β is very small and, furthermore, the accuracy is rather low. For each eigenvalue and each fixed mesh we recognize a critical value of the parameter such that greater values of β produce spurious eigenvalues. The behavior of these eigenvalues clearly reproduces that of the eigenvalues in Figure 1 (right) referring to Case 2. The results are plotted with a different perspective depending on the fact that the results now depend also on the computational mesh. The right bottom plot of Figure 11 highlights a phenomenon which already appears in Figure 10(i). Indeed, we see that the red line corresponding to the fourth computed eigenvalue for N = 400 lies along an hyperbola until β = 65 where it reaches the value 5 associated with second and third exact eigenvalues. Between β = 65 and β = 55 the red line remains close to 5, then decreasing β it follows a different hyperbola until it reaches the expected value for β = 35.

Conclusions
In this paper we have discussed how numerically computed eigenvalues can depend on discretization parameters. Section 2 shows the dependence on α and β of the eigenvalues of (1) when A and B have the forms (3) and (4), respectively. In Section 3 we have studied the behavior of the eigenvalues of the Laplace operator computed with the Virtual Element Method. The presence of two parameters resembles the abstract setting of Section 2; even if assumptions satisfied by the VEM matrices are more complicated than the ones previously discussed, the numerical results are pretty much in agreement. The present work opens the question of a viable choice of the parameters for eigenvalue computations when the discretization scheme depends on a suitable tuning of them (such as in the case of VEM).