solver_eigen_scipy_sym

[valvikby]

Hi, everyone.
Please explain, what is the general rule for choosing the parameters k and sigma in solver_eigen_scipy_sym()?

[kinnala]

It is basically just a simple wrapper to scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html

I think k is the number of eigenvalues. I cannot remember why sigma=10 by default, but you can check scipy docs to see for any hints on how sigma should be chosen (or set to None).

[gdmcbain]

It's the shift; sigma shifts the whole eigenvalue problem. Use your best guess at the eigenvalue of interest . The underlying solver is iterative, basically some refinement of the power iteration, so will work faster the greater the relative magnitude of the leading eigenvalue to the second. That basic iteration is no good if the eigenvalue of interest isn't the largest or if there are two almost equally large. See the Wikipedia article on inverse iteration, where the shift is denoted μ.

[gdmcbain]

I cannot remember why sigma=10 by default,

I believe that this was a good choice for one of the examples! :) Nothing really universal about it.

Well, sometimes if boundary conditions are applied sloppily we get spurious eigenvalues, usually at or near zero (or infinity). In structural mechanics, we're often interested in the lowest eigenvalues as representing the natural frequencies of vibration, but not the spurious modes corresponding to free translation if the piece isn't properly pinned down. A popular hack is to not bother with boundary conditions, i.e..leave the piece free in space, but shift the eigenvalue problem far enough away from zero so that the eigenvalue most favoured by the inverse iteration is the fundamental.

[gdmcbain]

An illustrative toy:

import numpy as np
from scipy.linalg import eigh
from skfem import *
from skfem.models.poisson import mass, laplace

mesh = MeshLine(np.linspace(0, 1, 2**3))
element = ElementLineP1()
basis = Basis(mesh, element)

L = asm(laplace, basis)
M = asm(mass, basis)

modes = 3
print("Dense spectrum:", eigh(L.A, M.A)[0][:modes])
print("Sparse spectrum:", solve(L, M)[0][:modes])
print("Shift 42:", solve(L, M, sigma=42)[0][:modes])

Output:

Dense spectrum: [9.48057857e-14 1.00363548e+01 4.21934166e+01]
Sparse spectrum: [1.00363548e+01+0.j 3.55271368e-15+0.j 4.21934166e+01+0.j]
Shift 42: [ 4.21934166e+01+0.j  1.00363548e+01+0.j -3.55271368e-14+0.j]

scipy.linalg.eigh computes the complete spectrum whereas the scipy.sparse.linalg eigensolves as wrapped by skfem.utils.solve only compute a few (k) leading modes, where "leading" means nearest sigma. In this toy example, without boundary conditions, the heat equation has an isothermal mode with zero eigenvalue as shown by eigh. This is pushed down the list by sigma.