Find zeros of a characteristic polynomial matrix with Python

Question

Given the symmetric matrix N x N Cand the diagonal matrix N x N I, we find solutions to the equation det(λI-C)=0. In other words, (generalized) eigenvalues ​​must be found C.

I know several ways to solve this problem in MATLAB using built-in functions:

First way:

function lambdas=eigenValues(C,I)
    syms x;
    lambdas=sort(roots(double(fliplr(coeffs(det(C-I*x))))));

The second way:

[V,D]=eig(C,I);

However, I need to use Python. Similar functions exist in NumPy and SymPy, but according to the docs ( numpy , sympy ), they accept only one C matrix as input. Although, the result is different from the result obtained by Matlab. In addition, symbolic solutions created by SymPy do not help. Maybe I'm doing something wrong? How to find a solution?

---

Example

• MATLAB:

% INPUT

I =

     2     0     0
     0     6     0
     0     0     5
C =

     4     7     0
     7     8    -4
     0    -4     1

[v,d]=eig(C,I)

% RESULT

v =

    -0.3558   -0.3109   -0.5261
     0.2778    0.1344   -0.2673
     0.2383   -0.3737    0.0598


d =

       -0.7327    0         0
        0    0.4876         0
        0         0    3.7784

• Python 3.5:

% INPUT

I=np.matrix([[2,0,0],
                 [0,6,0],
                 [0,0,5]])

C=np.matrix([[4,7,0],[7,8,-4],[0,-4,1]])


np.linalg.eigh(C)

% RESULT

(array([-3., 1.91723747, 14.08276253]),


matrix(

       [[-0.57735027,  0.60061066, -0.55311256],

        [ 0.57735027, -0.1787042 , -0.79670037],

        [ 0.57735027,  0.77931486,  0.24358781]]))

Answer 1

At least if it Ihas positive diagonal entries, you can simply solve the converted system:

# 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 example. IMPORTANT: The use np.arrayof Iand C, np.matrixit will not work.

>>> 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.

Answer 2

, , , Ix = x.

, , Cx = λIx, , , , , Numpy eig(C).

C , , numpy.linalg.eigh, .

, , , , Kx = ω²Mx, scipy.linalg.eigh, .

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

Numpy w/r , , ( ...) , , , Scipy eigh.

, , , , , undefined (, , - , ) - scipy.linalg.eigh, (I ).

---

Ps: scipy.linalg.eigh (.. ) my, .

Answer 3

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

:

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

:

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

, .