How to get eigenvectors from eigenvalues?


87
views
0
7 weeks ago by
CD  
Hi,

I was using eigensolver to solve a Neumann eigenfunction problem, and it didn't converge for most 1<=N<=50, but for N=37, it printed out a list of 37 eigenvalues. I am trying to get the .pvd files which come associated with each eigenfunction on my domain, but I don't have all 37 now. Does anybody know how to get a particular eigenvector for FENICS reading in a particular eigenvalue? I feel like this should be possible, but all the functions I see in the documentation seem like they accept the solving conditions and return an eigenvalue and update the vector within the working directory to the eigenvalue.

Thanks in advance,
-C
Community: FEniCS Project

1 Answer


4
7 weeks ago by
Here's a piece of a code that creates a list of eigenvalues and a list of eigenfunctions.  You can write the eigenfunctions to PVD in the usual way (File('eigenfun3.pvd') << eigenfunctions[3]).  I do a normalization of my eigenfunction which you may, or may not, want to do.  Also, I am doing a symmetric problem, so ignore the imaginary part of the eigenvalues and eigenvectors.  Hope this helps.

A = assemble(a)
M = assemble(m)
# downcast to PETSc matrices
A = as_backend_type(A)
M = as_backend_type(M)
# Create eigensolver
eigensolver = SLEPcEigenSolver(A, M)
eigensolver.parameters["solver"] = "krylov-schur"
eigensolver.parameters["spectrum"] = "smallest magnitude"
eigensolver.parameters["problem_type"] = "gen_hermitian"
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = 0.3
eigensolver.parameters["tolerance"] = 1.0e-5
eigensolver.parameters["maximum_iterations"] = 1000
eigensolver.solve(20)
nconv = eigensolver.get_number_converged()
eigenvalues = []
eigenfunctions = []
for i in range(nconv):
    r, _, rx, _ = eigensolver.get_eigenpair(i)
    eigenvalues.append(r)
    # normalize eigenfunction so max = max magnitude = 1
    rxmax = rx.max()
    rxmin = rx.min()
    if rxmax < -rxmin:
        rx = rx/rxmin
    else:
        rx = rx/rxmax
    eigenfun = Function(V)
    eigenfun.vector()[:] = rx
    eigenfunctions.append(eigenfun)
​
Thanks a lot for the detailed reply--I had found success printing the .pvds of the eigenvectors by putting the declaration within the for loop. I see you're normalizing eigenvectors, if you don't mind would you mind telling me  why you were doing this? Is there a notion of normalizing eigenvalues as well? I ask because the eigenvalues which were returned from my solver seemed too regular to be accurate without some funky normalization going on.
written 7 weeks ago by CD  
1
In an eigenvalue problem, the eigenvalues are well-defined numbers and there is no normalizing them.  But the eigenfunctions are only determined up to a constant multiples.  If  $Lu=\lambda u$Lu=λu, then the equation remains true if you replace  $u$u by $cu$cu  for any constant $c$c. So you can normalize by choosing a convenient constant.  A common choice is to pick a norm of interest and set the constant c so that the norm of the eigenfunction is 1.  E.g., you might want the  $L_2$L2  norm or the  $H^1$H1  to be equal to 1.  I normalized so that the  $L_{\infty}$L norm (i.e., the maximum absolute value) was equal to 1.

What do you mean by saying that the eigenvalues returned by your solver were too regular?

written 7 weeks ago by Douglas N Arnold  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »