线性方程组求解:Jacobi 迭代法与 LU 分解法的应用
首先,我们给出线性方程组的系数矩阵 A 和常数向量 b:
A = np.array([[4, -1, 0],
[1, 6, -2],
[0, 2, 5]])
b = np.array([2, -4, 1])
接下来,我们来计算 Jacobi 迭代矩阵。Jacobi 迭代矩阵的定义如下:
D = np.diag(np.diag(A)) # A 的对角线元素构成的对角矩阵
L = -np.tril(A, k=-1) # A 的下三角矩阵
U = -np.triu(A, k=1) # A 的上三角矩阵
J = np.linalg.inv(D) @ (L + U)
其中,tril 和 triu 函数分别返回矩阵的下三角和上三角部分。
然后,我们可以计算 Jacobi 迭代矩阵的谱半径,判断迭代法是否收敛。谱半径的定义如下:
rho = np.max(np.abs(np.linalg.eigvals(J)))
如果 rho 小于 1,则迭代法收敛;如果 rho 大于等于 1,则迭代法不收敛。
最后,根据迭代法是否收敛,选择使用迭代法或 LU 分解法求解方程组。
下面是 Python 代码实现:
import numpy as np
# 线性方程组的系数矩阵 A 和常数向量 b
A = np.array([[4, -1, 0],
[1, 6, -2],
[0, 2, 5]])
b = np.array([2, -4, 1])
# 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)
# Jacobi 迭代矩阵的谱半径
rho = np.max(np.abs(np.linalg.eigvals(J)))
if rho < 1:
print('迭代法收敛')
# 使用迭代法求解方程组
x = np.zeros_like(b) # 初始化解向量
max_iter = 100 # 最大迭代次数
tol = 1e-6 # 迭代精度
for i in range(max_iter):
x_new = (b - (L + U) @ x) / np.diag(A)
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
print('方程组的解为:', x)
else:
print('迭代法不收敛')
# 使用 LU 分解法求解方程组
if np.linalg.det(A) != 0:
x = np.linalg.solve(A, b)
print('方程组的解为:', x)
else:
print('方程组无解或有无穷解')
注意,上述代码使用了 numpy 中的一些线性代数函数,需要提前导入 numpy 库。
原文地址: https://www.cveoy.top/t/topic/o9cW 著作权归作者所有。请勿转载和采集!