如何在Python中解决多项式特征值？

8 投票

2 回答

3235 浏览

提问于 2025-04-17 06:58

在我的Python代码中，我想解决一个多项式特征值问题：

A0 + lambda*A1 + lambda^2*A2 + lambda^3*A3 + .... = 0

这里的 An 是稠密矩阵，而 lambda 是一个常数。在Matlab中，可以使用 polyeig函数来解决这个问题。但在scipy中似乎没有类似的功能。目前我能想到的唯一方法是构造相应的伴随矩阵。这样就会形成一个等价的线性特征值问题，可以交给现有的scipy求解器来处理，不过这个问题会变得更大，而且我觉得它可能会比较不稳定。

有没有人能推荐一个现有的、开源的或免费的库来解决这个问题？我会很高兴使用一个可以通过f2py链接的Fortran库，或者一个可以通过C/C++链接的库来使用cython。

补充：如果有人对在Python中解决非线性特征值问题感兴趣，我自己写的代码可以在这里找到。请注意，我处理的是更一般的非线性特征值问题（因为它对lambda有非线性依赖）。要理解这个方法，请阅读代码注释中提到的论文。

2 个回答

这是我对在Python中实现polyeig的看法。作为测试案例，我使用了和Mathworks提供的例子一样的例子。

缩放：返回的特征向量和Matlab的结果不一样，因为它们是根据一个缩放因子来定义的。我不太确定Matlab使用了什么样的归一化方法，所以在这种情况下，我只是用每个向量的最大值进行了归一化。

排序：Matlab并不对特征值进行排序。在那个例子中，scipy似乎是按照和Matlab相同的顺序返回特征值的。你可以取消注释脚本中的几行代码来对特征值进行排序。

import numpy as np
from scipy import linalg

def polyeig(*A):
    """
    Solve the polynomial eigenvalue problem:
        (A0 + e A1 +...+  e**p Ap)x=0 

    Return the eigenvectors [x_i] and eigenvalues [e_i] that are solutions.

    Usage:
        X,e = polyeig(A0,A1,..,Ap)

    Most common usage, to solve a second order system: (K + C e + M e**2) x =0
        X,e = polyeig(K,C,M)

    """
    if len(A)<=0:
        raise Exception('Provide at least one matrix')
    for Ai in A:
        if Ai.shape[0] != Ai.shape[1]:
            raise Exception('Matrices must be square')
        if Ai.shape != A[0].shape:
            raise Exception('All matrices must have the same shapes');

    n = A[0].shape[0]
    l = len(A)-1 
    # Assemble matrices for generalized problem
    C = np.block([
        [np.zeros((n*(l-1),n)), np.eye(n*(l-1))],
        [-np.column_stack( A[0:-1])]
        ])
    D = np.block([
        [np.eye(n*(l-1)), np.zeros((n*(l-1), n))],
        [np.zeros((n, n*(l-1))), A[-1]          ]
        ]);
    # Solve generalized eigenvalue problem
    e, X = linalg.eig(C, D);
    if np.all(np.isreal(e)):
        e=np.real(e)
    X=X[:n,:]

    # Sort eigenvalues/vectors
    #I = np.argsort(e)
    #X = X[:,I]
    #e = e[I]

    # Scaling each mode by max
    X /= np.tile(np.max(np.abs(X),axis=0), (n,1))

    return X, e

if __name__=='__main__':
    M = np.diag([3,1,3,1])
    C = np.array([[0.4 , 0 , -0.3 , 0], [0 , 0  , 0 , 0], [-0.3 , 0 , 0.5 , -0.2 ], [ 0 , 0 , -0.2 , 0.2]])
    K = np.array([[-7  , 2 , 4    , 0], [2 , -4 , 2 , 0], [4    , 2 , -9  , 3    ], [ 0 , 0 , 3    , -3]])
    X,e = polyeig(K,C,M)
    print('X:\n',X)
    print('e:\n',e)
    # Test that first eigenvector and value satisfy eigenvalue problem:
    s = e[0];
    x = X[:,0];
    res = (M*s**2 + C*s + K).dot(x) # residuals
    assert(np.all(np.abs(res)<1e-12))

回答于 2025-04-17 由 Python大师

分享举报

这段讨论提到了一种通用的方法，可以把多项式特征值问题转化为广义特征值问题。之后，我们可以用 scipy的线性代数函数来解决这个问题。希望这对你有帮助！

回答于 2025-04-17 由 Python大师

分享举报

如何在Python中解决多项式特征值？

2 个回答

撰写回答