Python LU 分解求解线性方程组和 Jacobi 迭代法

import numpy as np

def pyLU(A, b, flag=0): # LU分解 n = np.size(b, 0) u = np.zeros((n, n), dtype=float) L = np.identity(n, dtype=float) u[0, :] = A[0, :] L[1:n, 0] = A[1:n, 0] / u[0, 0] for k in range(1, n): u[k:k+1, k:n] = A[k:k+1, k:n] - np.dot(L[k:k+1, 0:k], u[0:k, k:n]) L[k+1:n, k:k+1] = (A[k+1:n, k:k+1] - np.dot(L[k+1:n, 0:k], u[0:k, k:k+1])) / u[k, k]

# 解下三角方程组Ly=b
y = np.zeros((n, 1), dtype=float)
y[0] = b[0]
for k in range(1, n):
    y[k] = b[k] - np.dot(L[k:k+1, 0:k], y[0:k])

# 解上三角方程组Ux=y
x = np.zeros((n, 1), dtype=float)
x[n-1] = y[n-1] / u[n-1, n-1]
for k in range(n-2, -1, -1):
    x[k] = (y[k] - np.dot(u[k, k+1:n], x[k+1:n])) / u[k, k]

return x, L, u

if name == 'main': # 例2.9 A = np.array([[2, -1, 4, -3, 1], [-1, 1, 2, 1, 3], [4, 2, 3, 3, -1], [-3, 1, 3, 2, 4], [1, 3, -1, 4, 4]], dtype=float) b = np.array([[11, 14, 4, 16, 18]], dtype=float).T x, L, u = pyLU(A, b) print('x= {}'.format(x)) print('L= {}'.format(L)) print('u= {}'.format(u))

# 定义线性方程组的系数矩阵A和常数向量b
A = np.array([[2, 3, 4, 5],
              [3, 5, 2, 1],
              [4, 3, 12, 5],
              [5, 6, 7, 8]])
b = np.array([24, -5, 34, 33])

# 计算Jacobi迭代矩阵
D = np.diag(np.diag(A))
L = -np.tril(A, k=-1)
U = -np.triu(A, k=1)
J = np.linalg.inv(D) @ (L + U)

# 计算迭代矩阵的谱半径
eigenvalues = np.linalg.eigvals(J)
spectral_radius = max(abs(eigenvalues))

# 判断迭代法是否收敛
if spectral_radius < 1:
    print('迭代法收敛')
    # 进行迭代求解
    x = np.zeros_like(b)
    for i in range(100): # 设置迭代次数上限
        x_new = J @ x + np.linalg.inv(D) @ b
        if np.linalg.norm(x_new - x) < 1e-6: # 设置迭代停止条件
            break
        x = x_new
    print('方程组的解为:', x)
else:
    print('迭代法不收敛')

# 采用LU分解法求解方程组
LU, piv = np.linalg.lu(A)
if np.linalg.det(LU) != 0: # 判断方程组是否有唯一解
    x = np.linalg.solve(A, b)
    print('方程组的解为:', x)
else:
    print('方程组无解或有无穷多解'