It is well known that the approximation of eigenvalues and associated eigenfunctions of a linear operator under constraint is a difficult problem. One of the difficulties is to propose methods of approximation which satisfy in a stable and accurate way the eigenvalues equations, the constraint one and the boundary conditions. Using any non-stable method leads to the presence of non-physical eigenvalues: a multiple zero one called spurious modes and non-zero one called pollution modes. One way to eliminate these two families is to favor the constraint equations by satisfying it exactly and to verify the equations of the eigenvalues equations in weak ways. To illustrate our contribution in this field we consider in this paper the case of Stokes operator. We describe several methods that produce the correct number of eigenvalues. We numerically prove how these methods are adequate to correctly solve the 2D Stokes eigenvalue problem.