Generalized Nonsymmetric Eigensolver Python

Question

How do I solve a nonsymmetric eigenproblem. In terms of scipy.sparse.linalg.eigsh the matrix needs to be "real symmetric square matrix or complex Hermitian matrix A" (https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html). and the same for scipy.sparse.linalg.eigs. "M must represent a real symmetric matrix if A is real, and must represent a complex Hermitian matrix if A is complex. For best results, the data type of M should be the same as that of A" (https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs).

I need this for a vibroacoustic problem finite element problem. Which would have the matrix shape matrix description for vibroacoustic problem

I am using a windows computer.

ev, phi = scipy.sparse.linalg.eigs(A=stiff, k=nev, M=mass, which='SM') ev, phi = scipy.sparse.linalg.eigs(A=stiff, k=nev, M=mass, sigma=0) Before realizing the functions do not work for nonsymmetric problems.

The matrices are sparse.

Matlab Code:

k = load('k_global.mat');
m = load('m_global.mat');
[V,D] = eigs(k.array, m.array, 20,0);
D = diag(D)
natural_frequency = sqrt(D)/(2*pi)

Returns:

0.00000000000000 + 0.000356591992067188i
0.000668911454165071 + 0.00000000000000i
0.00000000000000 + 0.000973128785222090i
0.00222975851379527 + 0.00000000000000i
0.00246434216130016 + 0.00000000000000i
0.00000000000000 + 0.00372951940564144i
8.06883871646537 + 0.00000000000000i
64.7482150103242 + 0.00000000000000i
234.453670549319 + 0.00000000000000i
268.154072409059 + 0.00000000000000i
312.537263749716 + 0.00000000000000i
356.103849178590 + 0.00000000000000i
389.038117338274 + 0.00000000000000i
412.048267727649 + 0.00000000000000i
473.729345964820 + 0.00000000000000i
2996.35112385098 + 0.00000000000000i
3240.96766107255 + 0.00000000000000i
4186.42444133727 + 0.00000000000000i
4585.99172192305 + 0.00000000000000i
4794.52737053778 + 0.00000000000000i

Python code:

import pickle
import scipy
import numpy as np


if __name__ == '__main__':

    with open('../k_full.pickle', 'rb') as f:
        print('loading matrix K')
        k_global = pickle.load(f)
    with open('../m_full.pickle', 'rb') as f:
        print('loading matrix M')
        m_global = pickle.load(f)

    eigvalues, eigvectors = scipy.sparse.linalg.eigs(k_global, M=m_global, k = 20, which='SM')

    natural_frequency_Hz = np.sqrt(np.abs(eigvalues))/(2*np.pi)
    for i, nat_freq in enumerate(natural_frequency_Hz):
        print(f'[{i + 1: 3.0f}] : Freq = {nat_freq: 8.2f}')

Returns:

scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence: ARPACK error -1: No convergence (321 iterations, 4/20 eigenvectors converged)

Stiffness:

https://pastebin.com/Ri0rebyt

Mass

https://pastebin.com/mPKnTt8A

Answer 1

Despite your vibroacoustic link, you are (presumably) solving the undamped system

where M is the mass matrix and K is the stiffness matrix.

With solutions proportional to e^iωt this becomes

or

This is the system that you want to solve. However, M should be an invertible matrix, so you can rewrite this as

so looking for the eigensolutions of M^–1K.

If you look at the documentation ( https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs ), the routine scipy.sparse.linalg.eigs does NOT require the main matrix to be symmetric/hermitian and it doesn’t require the second matrix as an argument at all: that is optional.

For completeness, the final frequency f (in Hz) comes from

BTW, your matlab snippet is picking up some negative eigenvalues for lambda / complex values for omega. They are rather small in magnitude, so may be an artefact of floating-point operations. However, you should check for, and be careful how you deal with, such values.

EDIT - You have set k=20 (and which='SM'#allest). Thus, you aren't actually finding all the eigenvalues. If you increase k significantly you can find more, together with a warning about the scheme switching to the non-sparse version, scipy.linalg.eig. Given the relatively small size of your matrices, even if sparse, I recommend switching to that anyway. This then finds the particular eigenvalue that you want.

Code:

import numpy as np
import scipy

K = np.loadtxt( 'stiffness.txt' )
M = np.loadtxt( 'mass.txt' )
M1K = np.linalg.inv( M ) @ K
#eigenvalues, eigenvectors = scipy.sparse.linalg.eigs(M1K, k=20, which='SM')  # sparse version; 20 eigenvalues
eigenvalues, eigenvectors = scipy.linalg.eig(M1K)                             # general version -> 32 eigenvalues

natural_frequency_Hz = np.sqrt( np.abs( eigenvalues ) ) / ( 2 * np.pi )
for i, nat_freq in enumerate( natural_frequency_Hz ):
    print( f'[{i + 1: 3.0f}] : Freq = {nat_freq: 9.3f}' )

Output:

[  1] : Freq =  291049.779
[  2] : Freq =  291052.098
[  3] : Freq =  291056.220
[  4] : Freq =  291053.901
[  5] : Freq =  155654.905
[  6] : Freq =  155654.199
[  7] : Freq =  155632.590
[  8] : Freq =  155567.798
[  9] : Freq =  155596.123
[ 10] : Freq =  155588.053
[ 11] : Freq =  155630.579
[ 12] : Freq =  155578.601
[ 13] : Freq =  4794.527
[ 14] : Freq =  4585.992
[ 15] : Freq =  3240.968
[ 16] : Freq =  2996.351
[ 17] : Freq =  4186.424
[ 18] : Freq =   473.730
[ 19] : Freq =   411.719
[ 20] : Freq =   390.610
[ 21] : Freq =   356.104
[ 22] : Freq =   312.524
[ 23] : Freq =   268.154
[ 24] : Freq =   234.454
[ 25] : Freq =    64.748
[ 26] : Freq =     8.068
[ 27] : Freq =     0.022
[ 28] : Freq =     0.002
[ 29] : Freq =     0.002
[ 30] : Freq =     0.002
[ 31] : Freq =     0.002
[ 32] : Freq =     0.002