pyyeti.ytools.eig_si

pyyeti.ytools.eig_si(K, M, *, Xk=None, f=None, p=10, mu=0, tol=1e-06, pmax=None, maxiter=50, verbose=True, rng=None)[source]

Perform subspace iteration to calculate eigenvalues and eigenvectors.

Parameters:

  • K (ndarray) – The stiffness (assumed symmetric).

  • M (ndarray) – The mass (assumed positive-definite).

  • Xk (ndarray or None) – Initial guess @ eigenvectors; # columns > p. If None, random vectors are generated internally; see rng below.

  • f (scalar or None) – Desired cutoff frequency in Hz. pmax will override this if set. Takes precedence over p if both are input.

  • p (scalar or None) – Number of desired eigenpairs (eigenvalues and eigenvectors). pmax will limit this if set. If f is input, p is calculated internally (from sturm()).

  • mu (scalar) – Shift value in (rad/sec)^2 units. See notes.

  • tol (scalar) – Eigenvalue convergence tolerance.

  • pmax (scalar or None) – Maximum number of eigenpairs; no limit if None.

  • maxiter (scalar) – Maximum number of iterations.

  • verbose (bool) – If True, print status message for each iteration.

  • rng (numpy.random.Generator object or None; optional) – Random number generator. If None, a new generator is created via numpy.random.default_rng(). Uniform deviates are generated via rng.random(). Supplying your own rng can be handy for parallel applications, for example, when you need repeatability. For illustration, the following creates a PCG-64 DXSM generator and initializes it with a seed of 1:

    from numpy.random import Generator, PCG64DXSM
rng = Generator(PCG64DXSM(seed=1))
Returns:

  • lam (ndarray) – Ideally, p converged eigenvalues.

  • phi (ndarray) – Ideally, p converged eigenvectors.

  • phiv (ndarray) – First p columns are phi, others are leftover iteration vectors which may be a good starting point for a second call.

Notes

The routine solves the eigenvalue problem:

\[K \Phi = M \Phi \Lambda\]

Where \(\Phi\) is a matrix of right eigenvectors and \(\Lambda\) is a diagonal matrix of eigenvalues.

This routine works well for relatively small p. Trying to recover a large portion of modes may fail. Craig-Bampton models with residual flexibility modes also cause trouble.

mu must not equal any eigenvalue. For systems with rigid-body modes, mu must be non-zero. Recommendations:

  • If you have eigenvalue estimates, set mu to be average of two widely spaced, low frequency eigenvalues. For example, mu = 5000 worked well when the actual eigenvalues were: [0, 0, 0, 0, .05, 15.8, 27.8, 10745.4, …]

  • mu = -10 has worked well.

  • mu = 1/10 of the first flexible eigenvalue has worked well.

It may be temping to set mu to a higher value so a few higher frequency modes can be calculated. This might work, especially if you have good estimates for Xk. Otherwise, it is probably better to set mu to a lower value (as recommended above) and recover more modes to span the range of interest.

In practice, unless you have truly good estimates for the eigenvectors (such as the output phiv may be), letting Xk start as random seems to work well.

Routine follows the basic algorithm as outlined in [1].

References

Examples

>>> from pyyeti import ytools
>>> import numpy as np
>>> k = np.array([[5, -5, 0], [-5, 10, -5], [0, -5, 5]])
>>> m = np.eye(3)
>>> np.set_printoptions(precision=4, suppress=True)
>>> w, phi, phiv = ytools.eig_si(k, m, mu=-1)
Iteration 1 completed
Convergence: 3 of 3, tolerance range after 2 iterations is [...
>>> print(abs(w))
[  0.   5.  15.]
>>> import scipy.linalg as linalg
>>> rng = np.random.default_rng()
>>> k = rng.normal(size=(40, 40))
>>> m = rng.normal(size=(40, 40))
>>> k = np.dot(k.T, k) * 1000
>>> m = np.dot(m.T, m) * 10
>>> w1, phi1 = linalg.eigh(k, m, subset_by_index=(0, 14))
>>> w2, phi2, phiv2 = ytools.eig_si(
...     k, m, p=15, mu=-1, tol=1e-12, verbose=False
... )
>>> fcut = np.sqrt(w2.max())/2/np.pi * 1.001
>>> w3, phi3, phiv3 = ytools.eig_si(
...     k, m, f=fcut, tol=1e-12, verbose=False
... )
>>> print(np.allclose(w1, w2))
True
>>> print(np.allclose(np.abs(phi1), np.abs(phi2)))
True
>>> print(np.allclose(w1, w3))
True
>>> print(np.allclose(np.abs(phi1), np.abs(phi3)))
True