Hello all, One more time, thanks for your support. I've *managed to identify and resolve the issue I was facing*, and I'd like to share the solution here for the benefit of others who might encounter similar challenges.
The problem stemmed from my approach of solving for the reciprocals of the eigenvalues, which I did just because the inverse of the mass matrix M was easier. Specifically, I was solving M x = 1/\lambda AA x. The question of why this was the problem remains unsolved. I decided to solve the exact problem directly, namely AA x = \lambda M x. To achieve this, I computed the inverse of the matrix AA= Bt * A^inv * B using a Conjugate Gradient (CG) solver. Subsequently, I solved for the exact eigenvalues, and to my satisfaction, I obtained results that aligned with those obtained from MATLAB. If anyone has insights into why solving for the reciprocal yields different results, I would greatly appreciate your input on this matter. Best, Najwa On Friday, May 31, 2024 at 11:41:27 AM UTC+3 Najwa Alshehri wrote: > Dear Wolfgang: > > Thank you for your answer. I have checked the matrices and I printed the > case for a 2 by 2 mesh from both Matlab and deal.II. With a simple change > of rows and columns, they are exactly the same ( see the files attached). > Not sure yet why they result in different results of the eigenvalue > problem. There is a very importent detail that I want to point out here. > When I solved on the very course mesh (2 by 2) I got the following error, > which is not the case when I solved on a 4 by 4 mesh for example. > > An error occurred in line <791> of file > </usr/local/include/deal.II/lac/arpack_solver.h> in function > void dealii::ArpackSolver::solve(const MatrixType1&, const > MatrixType2&, const INVERSE&, std::vector<std::complex<double> >&, > std::vector<ValueType>&, unsigned int) [with VectorType = > dealii::Vector<double>; MatrixType1 = > dealii::LinearOperator<dealii::Vector<double>, dealii::Vector<double>, > dealii::internal::LinearOperatorImplementation::EmptyPayload>; MatrixType2 > = dealii::LinearOperator<dealii::Vector<double>, dealii::Vector<double>, > dealii::internal::LinearOperatorImplementation::EmptyPayload>; INVERSE = > dealii::SparseDirectUMFPACK] > The violated condition was: > false > Additional information: > Error with dnaupd, info -2. Check documentation of ARPACK > > > Dear Mathias: > > Thank you for your answer. The matrices are exactly the same in both > codes. The different is that I solve for the recipricals in deal.II. To be > precise, > > let AA= Bt A^inv B > and let M be a mass matrix > > where B, A, M are the same matrices in both Matlab and deal.II. > > In Matlab: I solve eig(AA, M) (AA x = \lambda M x) then I seeks the > square root of the smalles \lambda > In dealii: I solve using ArpackSolver solving M x = 1/\Lambda AA x) > so I do first > M_inv_umfpack.initialize(M); > auto M = linear_operator(M); > auto M_inv = linear_operator(M, M_inv_umfpack); > then, I create a linear operator of the other matrices and the inverse of > A and perform AA= Bt * A^inv * B > lastly I solve : arpack_solver.solve( M, AA, M_inv_umfpack, > eigenvalues,eigenvectors); > > So I am seeking the largest eigenvalue and \lambda = 1/largest eigenvalue > and the infsup constant is the square root of that. > > Could solving for the recipricals is the issue? if yes what should I fill > in the third intery of the solving if Iam doing > arpack_solver.solve( AA, M, ????????, eigenvalues,eigenvectors); > > One more time, your questions and help are highly appreciated. > Best, > Najwa > On Friday, May 31, 2024 at 12:40:04 AM UTC+3 Matthias Maier wrote: > >> Do you use the same "mass matrix" in both cases, meaning do you call >> eig(A, M) where A is the stiffness matrix and M is the mass matrix in >> matlab? >> >> >> You will need to compute the spectrum of the eigenvalue problem >> >> A x = \lambda M x >> >> >> Best, >> Matthias >> >> >> On Thu, May 30, 2024, at 10:55 CDT, Najwa Alshehri <[email protected]> >> wrote: >> >> > Dear all, >> > >> > >> > Thank you for your always support and help. >> > >> > >> > I am encountering an issue while solving an eigenvalue problem related >> to >> > the computation of the discrete inf-sup constant for a saddle point >> > problem. Specifically, I am solving the following system: >> > >> > >> > Bt A^-1 B eigenvector = eigenvalue M eigenvector, >> > >> > >> > where Bt A^-1 B is clearly symmetric and at least semi-positive >> definite, >> > and M is a symmetric positive definite matrix. >> > >> > >> > The problem arises when comparing the inf-sup constants obtained using >> > Matlab and dealii (solving on the same mesh ): >> > >> > >> > 1. Matlab: Using the "eig" solver, the inf-sup constant decays to zero >> > as the mesh is refined. >> > 2. deal.II: Using the "ArpackSolver", the inf-sup constant is bounded >> > away from zero. >> > >> > The discrepancy is puzzling, as both methods should give consistent >> > results. To illustrate this issue, I have attached a figure comparing >> the >> > two results and included in this post the two data files obtained from >> > Matlab and deal.II. >> > >> > >> > I would greatly appreciate any insights or suggestions on what might be >> > causing the discrepancy. >> > >> > >> > Your help is highly appreciated, >> > >> > Najwa >> > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/6b6301a3-ef28-4410-8a0c-b8ea764486a7n%40googlegroups.com.
