行列の固有値計算のJacobi法をpythonで実装した

数値計算ライブラリのnumpyはほんとに便利。

import numpy as np

def jacobi(A,N,check):
    B = np.fabs(A - np.diag(list(np.diag(A))))
    nondiagmax = np.max(B)
    while nondiagmax > check:
        k = int(np.argmax(B) / 3)
        m = np.argmax(B) % 3
        cos2phi = np.fabs(A[k,k] - A[m,m]) / \
                np.sqrt(4.0 * A[k,m]**2 + (A[k,k] - A[m,m])**2)
        cosphi = np.sqrt((1.0 + cos2phi) / 2.0)
        sinphi = np.sign(-A[k,m] * (A[k,k] - A[m,m])) * \
                np.sqrt((1.0 - cos2phi) / 2.0)
        P = np.matrix(np.identity(N))
        P[k,k] = cosphi
        P[k,m] = sinphi
        P[m,k] = -sinphi
        P[m,m] = cosphi
        A = ((P.T).dot(A)).dot(P)
        B = np.fabs(A - np.diag(list(np.diag(A))))
        nondiagmax = np.max(B)
    print(A)
    print(np.diag(A))

if __name__ == "__main__":
    EPS = 1e-8
    N = 3
    A = np.matrix([[3.0,0.01,0.1],[0.01,2.0,0.1],[0.1,0.1,1.0]])
    print(A)
    jacobi(A,N,EPS)

結果より、固有値は、[ 3.00521156 2.00950971 0.98527872]