反幂法

import numpy as np


def lufenjie(a):
    n = a.shape[0]
    L = np.mat(np.zeros((n, n), dtype=float))
    for i in range(n):
        L[i, i] = 1
    u = np.mat(np.zeros((n, n), dtype=float))
    u[0, :] = a[0, :]
    L[1:n, 0] = a[1:n, 0]/u[0, 0]
    for item in range(1, n-1):
        u[item, item] = a[item, item] - np.sum(np.multiply(L[item, 0:item], u[0:item, item].T))
        for j in range(item+1, n):
            u[item, j] = a[item, j] - np.sum(np.multiply(L[item, 0:item], u[0:item, j].T))
            L[j, item] = (a[j, item] - np.sum(np.multiply(L[j, 0:item], u[0:item, item].T)))/u[item, item]
    u[n-1, n-1] = a[n-1, n-1] - np.sum(np.multiply(L[n-1, 0:n-1], u[0:n-1, n-1].T))
    return L, u

def fgauss(A, b):
    n = A.shape[0]
    zengguan = np.hstack((A, b.T))
    ra = np.linalg.matrix_rank(A)
    rz = np.linalg.matrix_rank(zengguan)
    temp1 = rz - ra
    if temp1 > 0:
        print("无一般意义下的解，系数矩阵与增广矩阵的秩不同")
        return
    if ra == rz:
        if ra == n:
            x = np.mat(np.zeros(zengguan.shape[0], dtype=float))
            for p in range(0, n):
                for k in range(p+1, n):
                    m = zengguan[k, p]/zengguan[p, p]
                    zengguan[k, p:n+1] = zengguan[k, p:n+1] - m * zengguan[p, p:n+1]
            b1 = zengguan[0:n, n]
            a1 = zengguan[0:n, 0:n]
            x[0, n-1] = b1[n-1]/a1[n-1, n-1]
            for i in range(n - 2, -1, -1):
                try:
                    x[0, i] = (b1[i] - np.sum(np.multiply(a1[i, i+1:n], x[0, i+1:n]))) / (a1[i, i])
                except:
                    print("错误")
    return x

A = np.mat([[-1, 2, 1],
            [2, -4, 1],
            [1, 1, -6]], dtype=float)
ladam = -6.42
I = np.mat(np.eye(3))
b = A - ladam*I
L, U = lufenjie(b)
err = 1
m_r = []
k = 1
while err > 0.0001:
    print(k,end=' ')
    if k==1:
        u = np.mat([1, 1, 1])
    else:
        u = fgauss(L, z)
    print('u=', u, end=' ')
    y1 =fgauss(U, u)
    y2 = abs(y1.copy())
    a = y2.argmax()
    m1 = y1[0, a]
    z = y1/m1
    print('m=', m1, end=' ')
    print('z=', z)
    m_r.append(m1)
    if k > 2:
        err=abs(m_r[-1]-m_r[-2])
    k = k+1
print('特征值λ=', ladam+1/m1)
print('特征向量为：', z)