{'fraction_non_alphanumeric': 0.02845480591859963, 'fraction_numerical': 0.009615762000078496, 'mean_word_length': 4.746504285069914, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chebyshev spectral collocation method approximations of the Stokes eigenvalue problem based on penalty techniques

The numerical solution strategies for the Stokes eigenvalue problem based on the use of penalty formulations are investigated in this study. It is shown that the penalty method approach can successfully be adapted for the eigenproblem to rectify the associated problems such as the existence of zero diagonal entries in the resulting algebraic system. Two different schemes, namely, the standard penalisation with a small penalty parameter, and the iterative penalisation that enables relatively large parameters, are implemented. The employment of the latter leads to a so-called inhomogeneous generalised eigenvalue problem which requires a special attention. A feasible solution strategy is presented which is adapted from a procedure based on Newton’s method proposed for the corresponding standard (inhomogeneous) eigenvalue problems. Concerning the spatial discretisation, among other possible options, the Chebyshev spectral collocation method based on expanding the unknown fields in tensor product of Chebyshev polynomials is employed. It is shown that the method constitutes a novel way of efficiently examining the approximate eigensolutions of the Stokes operator with the use of Chebyshev spectral collocation method directly, without a decoupling of velocity and pressure.


Introduction
The Stokes eigenvalue problem is a subject of an extensive research due to its significance on both fundamental and practical grounds. For instance, the eigenmodes of the Stokes operator form a natural basis in the analysis of homogeneous component of any flow. They are used in describing the fluctuating part of a turbulent flow. Moreover, the Stokes eigenvalue problem is used as benchmark for analysing convergence and accuracy of the numerical algorithms designed in fluid dynamics (see, e.g, [9,15] and the references therein). There are numerous works devoted to approximating the Stokes eigenproblem based on different methodologies. Among them are finite element methods [1,8,12,18], mesh free methods based on radial basis functions [6], spectral Chebyshev methods based on decoupling the velocity and pressure operators [10,11], and spectral Lagrange method using a staggered grid system [4].
There are several difficulties in approximating accurately the eigenmodes of the Stokes problem both in two and three dimensional domains, such as the lack of theoretical knowledge of the corresponding spectrum, and algorithmic limitations due to the size of the resulting algebraic system. These difficulties are emphasised when the problem is formulated in primitive variables, velocity and pressure, with the solenoidal condition, as the discretisation of the operator leads to a generalised eigenvalue problem with zero diagonal entries in the resulting algebraic system. It is well known that these and similar issues arising in solving the corresponding source problems, can be resolved with the use of penalty methods in which a perturbed form of the problem is considered. Moreover, using a penalty method allows one to avoid the problem in the choice of the degrees of freedom for the pressure (which is determined up to an additive constant) for which, a zero mean value constraint is introduced.
In this paper, we investigate numerical solution strategies for the Stokes eigenvalue problem based on the use of penalty formulations. The motivation is to extend the widely used application of penalisation techniques to the Stokes eigenvalue problem. We show that the penalty method approach can successfully be adapted for the eigenproblem to deal with associated issues. Two different schemes, namely, the standard penalisation with a small penalty parameter, and the iterative penalisation that allows the use of relatively large parameters, are implemented. The employment of the iterative method leads to a so-called inhomogeneous generalised eigenvalue problem which cannot be treated in the classical framework of eigenproblems. We adapt an efficient solution strategy from a procedure based on Newton's method proposed for the corresponding standard (inhomogeneous) eigenvalue problems in [13]. We employ the Chebyshev spectral collocation method (CSCM) based on expanding the unknown fields in tensor product of Chebyshev polynomials for the spatial discretisation. As a matter of fact, the idea introduced here can be applied with any type of discretisation, however, the simplicity and efficiency of the CSCM make it a viable option for the demonstration of the idea. More importantly, we show that the proposed idea constitutes an efficient way of numerically examining the eigensolutions of the Stokes operator with the use of Chebyshev collocation approximation directly, that is, without a decoupling of velocity and pressure.

Problem statement
The Stokes eigenvalue problem consists of finding [u, p, λ] on a bounded and polyhedral domain Ω ⊂ R d , d = 2, 3, where u : Ω → R d , with u = 0, is the displacement or velocity field, p : Ω → R is the pressure, and λ ∈ R, such that in Ω, where ν > 0 is a physical parameter. The eigenproblem that seeks the pressure eigenfunctions (a mixed eigenvalue problem of second type, see [2]) will not be considered in this work. It is well known that the eigenvalues satisfy with the associated eigenfunctions which are assumed to satisfy where (·, ·) denotes the standard L 2 inner product.
As we aim to focus on the demonstration of the idea that the Stokes eigensolutions can be approximated by the use of different penalisation techniques, we consider (1) on both two-and three-dimensional configurations, namely, a square domain and a cube domain. For the former case, the velocity field components are denoted by u and v, whereas, the third dimensional component is denoted by w in the latter case. We note that the accurate determination of the Stokes eigensolutions even on simple domains is a challenging problem. On the other hand, an extension to more general domains can be made by using a domain decomposition technique or a coordinate transformation. Nevertheless, the numerical strategies devised in this study are independent of the choice of domains.

Numerical approximation
We explain the spatial discretisation of the Stokes eigenproblem (1) which is carried out by a collocation approach which is explained in Section 3.1. The discretisation will lead to an algebraic eigenvalue problem that, in general, can be solved in numerous ways. For the reason to be clear later, we will present a procedure based on Newton's method for approximating eigenvalues in Section 3.2.

The Chebyshev spectral collocation method
The method we consider is based on requiring the numerical approximation of each unknown to be exactly satisfied on the abscissae of the extreme points of the Chebyshev polynomials. In this approach, each function spans the whole domain under consideration and thus, the derivatives of the function depend on the entire discretisation.
The collocation points are given as  The n−th derivative of Φ(x) is then approximated by The first derivatives at the collocation points satisfy C Now, the discrete values of the first derivative of the function Φ N can be obtained as This equation can be written in a matrix-vector form as where D (1) is called the first order Chebyshev spectral differentiation matrix, which is of size (N + 1) × (N + 1). Φ N in (4) is an array of size (N + 1) whose components are Φ N (x i ), and d dx (Φ N ) is the array that contains Φ In order to minimise the round-off errors for the calculation of the first derivatives, the diagonal entries d ii are computed as [16] The n−th order derivative of the function Φ(x) is now approximated as where D (n) = [D (1) ] n , that is, n-times matrix multiplication of D (1) . In general, D (n) is referred as the n-th order Chebyshev spectral differentiation matrix. The use of matrix multiplication for higher order derivatives, and the use of Equation (5) for obtaining diagonal entries, lead to a significantly greater accuracy in the computation of second and higher order derivatives for a wide range of functions. The Chebyshev spectral differentiation matrix for functions defined on an arbitrary interval [a, b] can be constructed by a linear transformation The utilisation of the Chebyshev spectral differentiation matrices to construct the approximate discrete operators for solving eigenvalue problems in several space dimensions is described in the subsequent subsections.

Two-dimensional formulation
In a two-dimensional configuration, we first set a tensor product grid based on CGL points assuming that the domain is a square, and using the same polynomial degree N in each direction for all the unknowns. Let D (i) , i = 1, 2, denote the Chebyshev spectral differentiation matrices in the first and second order differentiation in each spatial dimension. Next, we substitute the approximations u N , v N and p N to u, v and p, respectively, into the system of equations describing the Stokes eigenproblem (1), and approximate the differential operators defining this problem. Then, the discretised equations are written as Here, the (N + 1) 2 × (N + 1) 2 matrices K, G x and G y are defined with the use of the Kronecker product as where I N denotes the identity matrix of order (N +1). The approximation vectors are of order (N + 1) 2 , and are computed in the following pattern The previous equations form the following generalised matrix eigenvalue problem where Here we note that the matrix K associated to the Laplace operator is not symmetric, and thus the block matrix L is not symmetric.
Let us write the reduced system corresponding to (11) after imposing the boundary conditions as

Three-dimensional formulation
The problem within a cubical domain is tackled in an analogous way as in the two-dimensional case. We assume without loss of generality that the same polynomial degree N is used in all directions; thus, the first and second order Chebyshev spectral differentiation matrices D N , respectively, are computed and are used to approximate the derivatives in all coordinates. If the approximations u N , v N , w N and p N to u, v, w and p, respectively, are substituted into the system of equations describing the Stokes eigenproblem (1), then the discretised equations are written in matrix-vector form as Here, the matrix associated to the Laplace operator in three dimensions is of order (N + 1) 3 , and it is defined as The matrices G ′ x , G ′ y , and G ′ z are defined as Next, we can write the following generalised matrix eigenvalue problem in the three-dimensional setting as where As before, let us write the reduced system corresponding to (13) with the boundary conditions imposed as

Eigenvalue approximation
As the eigenvalue approximation is established in the same manner for both two-and three-dimensional problems, it is sufficient to consider the former case. Thus, let us consider the reduced algebraic eigenproblem (12). This problem is a system of non-linear equations written as for the unknowns λ and φ. Here, the matrix R, in the term φ T Rφ, has been introduced for the normalisation of the eigenvectors corresponding to the velocity field associated to (2).
For convenience, we can reformulate this system as If we introduce the partitioned vector Λ defined by then the problem can be viewed as a non-linear equation in the form F (Λ) = 0. The Jacobian of this system is We can then use Newton's method for the solution of the non-linear problem, with a suitable initial estimate Λ 0 , written in the following way: Naturally, the inefficient way of inverting the Jacobian matrix is avoided, and a linear system is solved for the increment vector that is added to the previous iterate at each step.

Penalty methods for the Stokes eigenproblem
The penalty method is a widely used approach in incompressible Stokes and Navier-Stokes models for relaxing the solenoidal condition. It has some advantages, such as the possibility of condensing discontinuous pressures and writing the problem in terms of the velocity only (see, e.g., [5,14]). Below, we present the extension of the application of the penalisation idea to approximate the spectrum of the Stokes operator.

The classical penalty problem
The idea is to approximate the solution of the Stokes eigenproblem by the solution to the penalised problem where ε is a penalty parameter chosen such that 0 < ε/ν < 1.
The CSCM discretisation steps are followed as before yielding where the modified matrix L ε is now given as The reduced system after the imposition of the boundary conditions is written in the form The effect of removing zero diagonal entries of the penalisation can readily be seen in Figure 2, where the structures of L and L ε which are calculated on a square domain with N = 16, are visualised using the MATLAB function spy. It is well known that the solutions [u, p] and [u ε , p ε ] to the source problems that correspond to (1) and (17), respectively, satisfy where C 1 is a positive constant independent of ε. That is, the convergence of both u ε and p ε (for the source problem) is O(ε). Our calculations have revealed the fact that the convergence order is the same for the eigenvalue problem. Furthermore, we have observed that the approximations obtained from the classical penalty method satisfy for a positive constant C 2 . We do not provide a proof for this estimate, however, it can be done following the perturbation analysis given in [13].

The iterative penalty method
The iterative penalisation idea we consider here, is proposed in [5] for the Stokes source problem. In this approach, the penalised equations are solved in each iteration with the addition of the residual of the incompressibility equation of the previous iteration. This method allows the use of large penalty parameters leading to a system with better conditioning.
The resulting problem reads as: given initially p 0 ε ; find u ε i , p i ε , and λ i ε such that for i = 1, 2, . . .. The discretised system of (18) can be written as at each iteration i, where ϕ i−1 ε is the inhomogeneity vector partitioned as Clearly, the first iterate, that is, the case i = 1, with p 0 ε = 0, and thus ϕ 0 ε = 0, corresponds to the classical penalty problem and can be solved as described in the previous section. On the other hand, the next iterations where ϕ i−1 ε = 0 are in the form of an inhomogeneous generalised eigenvalue problem which cannot be solved as a standard eigenvalue problem. A natural contender for solving such a problem is a procedure based on Newton's method proposed for the corresponding inhomogeneous but standard eigenvalue problems in [13].
To begin with, we assume that a normalisation condition in the form (φ i ε ) T Rφ i ε = 1, accompanies Equation (19), since the eigenproblem we consider corresponds to finding velocity eigenfunctions. Consequently, we can write the iterative problem in the form This problem can be viewed as a non-linear equation in the form H(Λ i ε ) = 0, where the partitioned vector Λ i ε is defined by We can then introduce Newton's method for the solution of the non-linear problem at each iteration i = 1, 2, . . . , yielding the nested algorithm in the following way: where J H is the Jacobian of the system H(Λ i ε ) = 0. A convergence criteria is set at each inner iteration (with respect to j), i.e. Newton's method, and then the outer cycle (with respect to i), i.e. the iterative penalisation, is carried out until a corresponding condition for the convergence is met.

The coupled iterations
Instead of following the nested iterative scheme (20), we may implement a combined iterative scheme defined as starting from an initial ϕ 0 ε formed by a given p 0 . In this way, a computationally cheaper scheme is obtained that is also capable of approximating the eigenspectrum of the Stokes operator.
It is shown in [5] that for the corresponding source problem, the error of the iterative penalty method ν 1/2 ∇u − ∇u ε i + ν −1/2 p − p i ε is of order O(ε i ), so it tends to zero both if ε → 0 or i → ∞. For the Navier-Stokes source problem, it is possible to couple the iterations due to penalisation and to non-linearity as explained for the eigenvalue problem. It is shown in [5] that convergence is driven by the slower of the two errors. We have experimentally observed the same behaviour for the eigenvalue problem.

Numerical examples
In this section we present the numerical results for approximating the Stokes eigensolutions on the square domain [0, 1] 2 and the cube domain [−1, 1] 3 . The case ν = 1 is considered, and the results are obtained with different penalisation techniques described in Section 4. The computations are carried out by a computer program created by us, using MATLAB. In all the simulations presented below, we have taken 17 collocation points in each spatial dimension, that is, the interpolation polynomials are of degree N = 16, for all cases. For the iterative penalty method, the convergence criteria has been set to 10 −12 for both Newton's method and the iterative penalisation method.
The eigensolutions to the Stokes problem for all cases we consider are not known analytically, and therefore, we take as reference values existing in the literature to compare our results.

The square domain 5.1.1. The classical penalty method
Firstly, we consider approximating the Stokes eigenproblem on the square using the classical penalty method described in Subsection 4.1. The first 10 eigenvalues are listed in Table 1, where the penalty parameter is taken as ε = 10 −6 .
The tabulated values show that the approximations obtained from the classical penalisation agree reasonably well with the reference values. Further, the plot of the eigenfunction associated to the minimum eigenvalue is given in Figure 3, revealing a good agreement with the existing ones in the literature (see, [17]). In order to investigate the convergence behaviour of the approximated eigenvalue with respect to the penalty parameter, the minimum eigenvalue (denoted by λ 1 ) is calculated for ε values varying from 10 −1 to 10 −10 . For this test, we consider the corresponding reference value as λ ref = 52.34469138411319. This value is obtained by manually removing the spurious pressure modes from the non-penalised problem (for the case N = 16). The variation of |λ 1 − λ ref |/λ ref with respect to ε is presented in Figure 4, from which it is inferred that the convergence is linear.

The iterative penalty method
Here we present numerical results obtained from the techniques described in Subsection 4.2, for the approximation of the iterative penalty problem (18).  A plot of the relative change in the approximated first eigenvalue, defined by |λ i − λ i−1 |/λ i against the number of the accumulated iterations i is provided in Figure 5. Two different penalty parameters ε = 10 −1 and ε = 10 −3 are tested. In both cases, the quadratic dependence on the iteration number at initial steps shows that Newton's method governs the convergence behaviour of the procedure. The jumps in the error correspond to the new iterations of the nested loops. The previous experiment for ε = 10 −1 and ε = 10 −3 is now repeated by implementing the iterative penalisation procedure in a combined loop defined in (21). The convergence behaviours are illustrated in Figure 6. In progressive iterations, the convergence is driven by the penalisation scheme where an almost linear profile is observed. This phenomena is more pronounced for the larger value of the penalty parameter where the relative errors are larger compared to the smaller value for all iterates as expected. For the larger ε case, the profile changes to linear after the sixth iteration. A similar tendency can be observed for the smaller ε case, however, with faster decrease in the residual.

The cube domain
In the previous subsection, we have focused on the square domain to show that the penalty formulations are successfully applied to solve the Stokes eigenvalue problem with the use of CSCM. Now we proceed to present the corresponding results for the cube domain.

The classical penalty method
The classical penalisation approach has been tested to approximate the eigensolutions of the Stokes operator on the cube domain. The approximations to the first 33 eigenvalues are listed in Table 2, together with their multiplicities. We have obtained these results using ε = 10 −6 . Comparing them with the reference values published in [11], one sees that they are reasonably well approximated, with exactly matching multiplicities. The minimum eigenvalue (having multiplicity three, and is denoted by λ 1 ), is calculated for ε values varying from 10 −1 to 10 −10 . In the same way as in the square domain case, for this test, we take the corresponding reference value as λ ref = 15.543353940134807 which is obtained by manually removing the spurious pressure modes from the non-penalised problem. The variation of |λ 1 − λ ref |/λ ref with respect to ε is depicted in Figure 7, to examine the convergence behaviour of the approximated eigenvalue with respect to the penalty parameter. The convergence rate to the reference value with respect to the penalty parameter is clearly linear, as can be seen from the figure. Before passing to the results of the iterative procedure, we remark on the conditioning of the systems for both two-and three-dimensional cases. This is studied by means of the condition numbers of L and L ′ , denoted by κ and κ ′ , respectively. Figure 8 illustrates the variation of each condition number (calculated in the 2norm) with the penalty parameter. As can be inferred from this figure, in each case, the condition number grows linearly with the penalty parameter, and the order of magnitude does not depend on the spatial dimension.

The iterative penalty method
In analogy with the square domain case, a plot of the relative change in the approximated eigenvalue defined by |λ i − λ i−1 |/λ i against the number of the accumulated iterations is given in Figure 9 for ε = 10 −1 and ε = 10 −3 . In both cases corresponding to different penalty parameter, the quadratic dependence on the iteration number at initial steps shows that Newton's method governs the convergence behaviour of the procedure. In progressive iterations, the convergence is driven by the penalisation scheme where an almost linear profile is observed. The results are in accordance with the corresponding two-dimensional case (see, Figure 5). Finally, we present the results obtained from the application of the iterative penalisation procedure in the combined loop given in (21). As before, we repeat the experiments for ε = 10 −1 and ε = 10 −3 . The variation of the relative error with respect to the iteration number k is given in Figure 10 for each case. We see that the convergence properties are similar to those of the two-dimensional counterpart, especially in the case of smaller ǫ. The rate of convergence changes from quadratic to linear gradually, with a prominent alteration in the case of larger penalty parameter, as before. On the other hand, the convergence is noticeably slower compared to the solution on the square domain, for both penalty parameter cases. Not unexpectedly, a slight distortion in the convergence curve is observed, for the larger penalty parameter case. These arguments address the well known alteration in the convergence behaviour of Newton's method when approximating a multiple root, as it is the corresponding case for this problem.

Conclusions
We have presented two different methods, namely, the classical method and iterative method, based on the penalisation idea applied to the Stokes eigenproblem. We have shown that both procedures circumvent the difficulties related to the eigensystem solution, and further, provide an efficient means of approximating the Stokes eigenproblem directly with the use of CSCM. The numerical calculations suggest that the classical penalty method converges linearly to the reference values, for the two-and three-dimensional examples, namely a square and a cube domain, we considered. The iterative penalty method allows the use of large penalty parameters, leading to a system with better conditioning. On the other hand, it necessitates a novel algorithm to deal with the non-linearity inherited by the inhomogeneous nature of the eigenproblem. We have implemented this novel procedure based on Newton's method for approximating the eigenvalues of the resulting inhomogeneous generalised problem. Our results concerning two different applications, the nested loops and the single loop options, reflect the characteristic behaviours of the iterative penalisation. The convergence is governed by Newton's method initially as the quadratic dependence on the iteration number addresses for the nested loops option. For the single loop, numerical results have revealed that the quadratic convergence is achieved up to a certain iteration, after which, the profile tends to be linear. These behaviours are apparent for the larger penalty parameter value, whereas for the smaller one, a rapid decrease in the residual is observed for both cases considered. In the coupled iterative scheme, a decrease in convergence rate is observed for the larger penalty parameter in the cube domain where the approximated eigenvalue corresponds to a multiple root, as expected.