用Python找到矩阵特征多项式的零点

3 回答

• 1

  至少如果 I 具有正对角线条目，您可以简单地解决转换后的系统：

  # example problem
  >>> A = np.random.random((3, 3))
  >>> A = A.T @ A
  >>> I = np.identity(3) * np.random.random((3,))
  
  # transform 
  >>> J = np.sqrt(np.einsum('ii->i', I))
  >>> B = A / np.outer(J, J)
  
  # solve
  >>> eval_, evec = np.linalg.eigh(B)
  
  # back transform result
  >>> evec /= J[:, None]
  
  # check
  >>> A @ evec
  array([[ -1.43653725e-02,   4.14643550e-01,  -2.42340866e+00],
         [ -1.75615960e-03,  -4.17347693e-01,  -8.19546081e-01],
         [  1.90178603e-02,   1.34837899e-01,  -1.69999003e+00]])
  >>> eval_ * (I @ evec)
  array([[ -1.43653725e-02,   4.14643550e-01,  -2.42340866e+00],
         [ -1.75615960e-03,  -4.17347693e-01,  -8.19546081e-01],
         [  1.90178603e-02,   1.34837899e-01,  -1.69999003e+00]])
  

  OP的例子 . 重要提示：必须使用 np.array s用于 I 和 C ， np.matrix 将无效 .

  >>> I=np.array([[2,0,0],[0,6,0],[0,0,5]])
  >>> C=np.array([[4,7,0],[7,8,-4],[0,-4,1]])
  >>> 
  >>> J = np.sqrt(np.einsum('ii->i', I))
  >>> B = C / np.outer(J, J)
  >>> eval_, evec = np.linalg.eigh(B)
  >>> evec /= J[:, None]
  >>> 
  >>> evec
  array([[-0.35578356, -0.31094779, -0.52605088],
         [ 0.27778714,  0.1343625 , -0.267297  ],
         [ 0.23826117, -0.37371199,  0.05975754]])
  >>> eval_
  array([-0.73271478,  0.48762792,  3.7784202 ])
  

  如果 I 有正面和负面条目，则使用 eig 而不是 eigh ，然后将平方根转换为 complex dtype .

  回复于 2024-04-25T21:16:04+08:00
• 0

  与其他答案不同，我假设用符号I表示单位矩阵，Ix = x .

  您想要解决的是Cx =λIx，即所谓的标准特征值问题，并且大多数特征值求解器解决了该格式中描述的问题，因此Numpy函数具有签名 eig(C) .

  如果您的C矩阵是对称矩阵并且您的问题确实是标准的特征值问题，我建议使用 numpy.linalg.eigh ，这是针对此类问题进行优化的 .

  相反，如果你的问题确实是一个广义的特征值问题，例如，频率方程Kx =ω²Mx，你可以使用 scipy.linalg.eigh ，它支持对称矩阵的那种问题陈述 .

  eigvals, eigvecs = scipy.linalg.eigh(C, I)
  

  关于特征值的差异，Numpy实现不对它们的排序提供任何保证，因此它可能只是一个不同的排序，但如果你的问题确实是一个普遍的问题（我不是身份矩阵......）解决方案当然是不同的，你必须使用 eigh 的Scipy实现 .

  如果差异在特征向量内，请记住特征向量在任意比例因子中是已知的，并且，排序可能是未定义的（但是，当然，它们的顺序与您具有特征值的顺序相同） - scipy.linalg.eigh 的情况略有不同，因为在这种情况下，特征值被排序，并且特征向量相对于第二个矩阵参数（在您的示例中为I）进行归一化 .

  ---

  Ps： scipy.linalg.eigh 行为（即，排序的特征值和归一化的特征向量）对于我用来解决标准特征值问题的用例非常方便 .

  回复于 2024-04-25T21:16:04+08:00
• 1

  使用SymPy：

  >>> from sympy import *
  >>> t = Symbol('t')
  >>> D = diag(2,6,5)
  >>> S = Matrix([[ 4, 7, 0],
                  [ 7, 8,-4],
                  [ 0,-4, 1]])
  >>> (t*D - S).det()
  60*t**3 - 212*t**2 - 77*t + 81
  

  计算 exact 根：

  >>> roots = solve(60*t**3 - 212*t**2 - 77*t + 81,t)
  >>> roots
  [53/45 + (-1/2 - sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3) + 14701/(8100*(-1/2 - sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3)), 53/45 + 14701/(8100*(-1/2 + sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(312469/182250 + sqrt(797521629)*I/16200)**(1/3), 53/45 + 14701/(8100*(312469/182250 + sqrt(797521629)*I/16200)**(1/3)) + (312469/182250 + sqrt(797521629)*I/16200)**(1/3)]
  

  计算 floating-point 根的近似值：

  >>> for r in roots:
  ...     r.evalf()
  ... 
  0.487627918145732 + 0.e-22*I
  -0.73271478047926 - 0.e-22*I
  3.77842019566686 - 0.e-21*I
  

  请注意，根是真实的 .

  回复于 2024-04-25T21:16:04+08:00