【数值分析】高斯-赛德尔方法、规范化幂法、原点移位法
题目
要求
代码实现过程不能调用任何库函数自带的“线性 方程组求解、特征值求解库函数”
利用高斯-赛德尔方法求解上述线性方程组
使用Python编程求解矩阵A与列向量b
import numpy as np
import sympy as spdef create_matrix(n = 12):"""计算n阶矩阵(n为偶数) n默认为12"""matrix = np.zeros((n, n)) # 初始化一个n*n矩阵# 填充矩阵for i in range(1, n+1):# 数组下标从0开始 matrix[i-1, i-1] = 3 # 主对角元元素if i < n:matrix[i-1, i] = -1 # 超对角元元素if i > 1:matrix[i-1,i-2] = -1 # 次对角元元素if i != n/2 and i != n/2 + 1:matrix[i-1, n-i] = 1/2 # 副对角元元素return matrixprint("--------- 矩阵A ---------")
matrix_A = create_matrix()
print(matrix_A)# 将 numpy 数组转换为 sympy 矩阵
# sympy_matrix_A = sp.Matrix(matrix_A)
# latex_A = sp.latex(sympy_matrix_A)
# print(latex_A)def create_matrix_b(n = 12):"""计算n维列向量b(n为偶数) 默认n=12 """matrix_b = np.zeros((n, 1)) # 初始化一个n*1矩阵matrix_b[0, 0] = 2.5matrix_b[-1, 0] = 2.5for i in range(1, n-1):matrix_b[i, 0] = 1.5matrix_b[n//2-1, 0] = 1 # 整除 有效索引matrix_b[n//2, 0] = 1 # 整除 有效索引return matrix_bprint("-----------向量 b-----------------")
matrix_b = create_matrix_b()
print(matrix_b)# 将 numpy 数组转换为 sympy 矩阵
# sympy_matrix_b = sp.Matrix(matrix_b)
# latex_b = sp.latex(sympy_matrix_b)
# print(latex_b)
#
运行上述程序,可得
高斯—赛德尔方法求解上述线性方程组
def gauss_seidel(A, b, x0, tol=1e-5, max_iter=100):"""高斯-赛德尔迭代法求解线性方程组 Ax = b参数:\nA -- 系数矩阵\nb -- 常数列向量\nx0 -- 初始解向量\ntol -- 容忍误差(误差限)\nmax_iter -- 最大迭代次数(防止不收敛)\n返回:x -- 迭代得到的解向量n -- 迭代次数"""D = np.diag(np.diag(A)) # 主对角元矩阵L = -np.tril(A, -1) # 下三角矩阵(不包含主对角元) lower triangular matrixU = -np.triu(A, 1) # 上三角矩阵(不包含主对角元) upper triangular matrixB = np.linalg.inv(D-L) @ U # 迭代矩阵# print(np.round(B, 4))# print(sp.latex(sp.Matrix(np.round(B, 4))))f = np.linalg.inv(D-L) @ bx = x0.copy() # 复制初始解向量for iter_num in range(max_iter):x_new = B @ x + f# 检查是否满足容忍误差if np.linalg.norm(x_new - x, np.inf) < tol:n = iter_num # 迭代次数return x_new, nx = x_new # 更新解向量print("不收敛")# 调用函数A = create_matrix()
b = create_matrix_b()x0 = np.zeros((12, 1))
x, n = gauss_seidel(A, b, x0)
print("迭代次数为: {},近似解为{}".format(n, x))
计算结果
规范化幂法求解迭代矩阵的模最大特征值
def normalized_power_method(A, x0, tol=1e-5, max_iter=100):"""规范化幂法求解矩阵A的按模最大特征值参数:\nA -- 迭代矩阵\nx0 -- 初始非零向量\ntol -- 容忍误差(误差限)\nmax_iter -- 最大迭代次数(防止不收敛)\n返回:lambda_max -- 模最大特征值\nx -- 迭代得到的特征向量\nn -- 迭代次数\n"""x = x0.copy()for iter_num in range(max_iter):m = np.linalg.norm(x, np.inf)x /= m # 归一化x_new = A @ x# 检查是否满足误差限if np.linalg.norm((x_new / np.linalg.norm(x_new, np.inf) - x), np.inf) < tol:lambda_max = np.linalg.norm(x_new, np.inf)return lambda_max, x, iter_numx = x_newprint("不收敛")# 调用函数
A = create_matrix()
D = np.diag(np.diag(A)) # 主对角元矩阵
L = -np.tril(A, -1) # 下三角矩阵(不包含主对角元) lower triangular matrix
U = -np.triu(A, 1) # 上三角矩阵(不包含主对角元) upper triangular matrixB = np.linalg.inv(D-L) @ U # 迭代矩阵x0 = np.ones((12, 1))
lambda_max, x, n = normalized_power_method(B, x0)
print("迭代矩阵模最大特征值: ",lambda_max)
print("模最大特征值的近似解: ", x)
print("迭代步数: ", n)
计算结果
原点移位法加速
# origin accelerate
def origin_acc_normalized_power_method(A, x0, lambda_0, tol=1e-5, max_iter=1000):"""使用原点平移法的加速算法\n规范化幂法求解矩阵A的按模最大特征值参数:\nA -- 迭代矩阵\nx0 -- 初始非零向量\nlambda_0 -- 试算参数\ntol -- 容忍误差(误差限)\nmax_iter -- 最大迭代次数(防止不收敛)\n返回:lambda_max -- 模最大特征值\nx -- 迭代得到的特征向量\nn -- 迭代次数\n"""A = A - lambda_0 * np.eye(12)x = x0.copy()list = []for iter_num in range(max_iter):m = np.linalg.norm(x, np.inf)x /= m # 归一化x_new = A @ xlist.append(np.linalg.norm(x_new, np.inf) + lambda_0)# 检查是否满足误差限if np.linalg.norm((x_new / np.linalg.norm(x_new, np.inf) - x), np.inf) < tol:lambda_max = np.linalg.norm(x_new, np.inf)for item in list:print(item)return lambda_max + lambda_0, x, iter_numx = x_newprint("不收敛")A = create_matrix()
D = np.diag(np.diag(A)) # 主对角元矩阵
L = -np.tril(A, -1) # 下三角矩阵(不包含主对角元) lower triangular matrix
U = -np.triu(A, 1) # 上三角矩阵(不包含主对角元) upper triangular matrixB = np.linalg.inv(D-L) @ U # 迭代矩阵x0 = np.ones((12, 1))
lambda_max, x, n = origin_acc_normalized_power_method(B, x0, 0.1)
print("迭代矩阵模最大特征值: ",lambda_max)
print("模最大特征值的近似解: ", x)
print("原点移位法迭代步数: ", n)
运行结果