广义非对称特征求解器 Python

Question

如何解决非对称特征值问题。就 scipy.sparse.linalg.eigsh 而言，矩阵需要是“实对称方阵或复埃尔米特矩阵 A”（https://docs.scipy.org/doc/scipy/reference/ generated/scipy.sparse。 linalg.eigsh.html）。 scipy.sparse.linalg.eigs 也是如此。 “如果 A 是实数，则 M 必须表示实数对称矩阵；如果 A 是复数，则 M 必须表示复厄米特矩阵。为了获得最佳结果，M 的数据类型应与 A 的数据类型相同”（https://docs .scipy.org/doc/scipy/reference/ generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs）。

我需要这个来解决振动声学问题和有限元问题。其矩阵形状为振动声学问题的矩阵描述

我使用的是 Windows 电脑。

ev, phi = scipy.sparse.linalg.eigs(A=stiff, k=nev, M=mass,which='SM') ev, phi = scipy.sparse.linalg.eigs(A=刚性, k=nev, M=质量, sigma=0) 在意识到这些函数不适用于非对称问题之前。

矩阵是稀疏的。

Matlab代码：

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)

退货：

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代码：

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}')

退货：

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

刚度：

https://pastebin.com/Ri0rebyt

弥撒

https://pastebin.com/mPKnTt8A

Answer 1

你（大概）正在解决系统问题

其中 M 是质量矩阵，K 是刚度矩阵。

解与 e^iωt 成正比，这就变成了

或

这就是您要解决的系统。然而，M应该是一个可逆矩阵，所以你可以将其重写为

因此寻找 M^–1K 的特征解。

如果您查看文档（https://docs.scipy.org/doc/scipy/reference/ generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs）例程 scipy.稀疏.linalg.eigs不要求主矩阵是对称/埃尔米特的，并且不需要第二个矩阵作为参数完全：这是可选的。

为了完整起见，最终频率 f（以 Hz 为单位）来自

代码：

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')

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: 8.2f}' )

输出：

[  1] : Freq =  4794.55
[  2] : Freq =  4585.98
[  3] : Freq =  4186.41
[  4] : Freq =  3240.97
[  5] : Freq =  2996.36
[  6] : Freq =   459.12
[  7] : Freq =   422.46
[  8] : Freq =   386.07
[  9] : Freq =   386.07
[ 10] : Freq =   305.43
[ 11] : Freq =   208.57
[ 12] : Freq =   208.57
[ 13] : Freq =   198.29
[ 14] : Freq =    65.21
[ 15] : Freq =     0.51
[ 16] : Freq =     0.29
[ 17] : Freq =     0.24
[ 18] : Freq =     0.24
[ 19] : Freq =     0.07
[ 20] : Freq =     0.11