有精度要求的雅可比矩阵迭代代码

下面是一个简单的 Python 实现，它使用幂法来计算最大特征值和对应的特征向量，并使用反幂法来计算最小特征值和对应的特征向量。

import numpy as np

def jacobi_iteration(A, tol=1e-6, max_iter=100):
    """
    计算矩阵 A 的特征值和特征向量，使用雅可比迭代法。
    """
    n = A.shape[0]
    V = np.eye(n)
    iter_num = 0
    while True:
        # 找到绝对值最大的非对角线元素
        p, q = np.unravel_index(np.abs(A - np.diag(np.diag(A))).argmax(), (n, n))
        if np.abs(A[p, q]) < tol:
            break
        # 计算旋转角度
        theta = 0.5 * np.arctan2(2 * A[p, q], A[p, p] - A[q, q])
        # 构造旋转矩阵
        R = np.eye(n)
        R[p, p] = np.cos(theta)
        R[q, q] = np.cos(theta)
        R[p, q] = -np.sin(theta)
        R[q, p] = np.sin(theta)
        # 更新 A 和 V
        A = np.dot(np.transpose(R), np.dot(A, R))
        V = np.dot(V, R)
        iter_num += 1
        if iter_num >= max_iter:
            raise Exception("Reach maximum number of iterations.")
    # 提取特征值和特征向量
    eigenvalues = np.diag(A)
    eigenvectors = V
    return eigenvalues, eigenvectors

def power_iteration(A, tol=1e-6, max_iter=100):
    """
    计算矩阵 A 的最大特征值和对应的特征向量，使用幂法。
    """
    n = A.shape[0]
    x = np.random.rand(n)
    x /= np.linalg.norm(x)
    iter_num = 0
    while True:
        y = np.dot(A, x)
        lambda_ = np.dot(x, y)
        if np.linalg.norm(y - lambda_ * x) < tol:
            break
        x = y / np.linalg.norm(y)
        iter_num += 1
        if iter_num >= max_iter:
            raise Exception("Reach maximum number of iterations.")
    return lambda_, x

def inverse_power_iteration(A, sigma, tol=1e-6, max_iter=100):
    """
    计算矩阵 A 的最小特征值（与 sigma 的差的绝对值最小）和对应的特征向量，使用反幂法。
    """
    n = A.shape[0]
    x = np.random.rand(n)
    x /= np.linalg.norm(x)
    iter_num = 0
    while True:
        y = np.linalg.solve(A - sigma * np.eye(n), x)
        mu = 1 / np.linalg.norm(y)
        if np.linalg.norm(x - mu * y) < tol:
            break
        x = mu * y
        iter_num += 1
        if iter_num >= max_iter:
            raise Exception("Reach maximum number of iterations.")
    return sigma + 1 / mu, x

使用示例：

# 构造一个对称矩阵
A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])

# 计算矩阵的特征值和特征向量
eigenvalues, eigenvectors = jacobi_iteration(A)
print("Eigenvalues:", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

# 计算矩阵的最大特征值和对应的特征向量
lambda_, x = power_iteration(A)
print("Max Eigenvalue:", lambda_)
print("Max Eigenvector:", x)

# 计算矩阵的最小特征值和对应的特征向量
mu, y = inverse_power_iteration(A, 0.5)
print("Min Eigenvalue:", mu)
print("Min Eigenvector:", y)