Python实现部分主元法下LU分解

Python实现部分主元法下LU分解

''' 《矩阵分析与应用》小作业1 实现部分主元法下的LU分解 by苗栋 程序大体介绍： 引入了numpy便于对数组的操作 ①寻找出一列中绝对值最大的元素作为主元并进行数组的行交换，并将L主对角线元素置为1 ②构造两个for循环，取出L和A对应的值并保存 L[j, i] = A[j, i] / A[i, i] A[j, :] = A[j, :] - L[j, i] * A[i, :] ③最后打印结果 ''' import numpy as np def LU_with_partial_pivoting(A): #转换为float类型的数组并初始化L A = np.array(A, dtype=np.float64) L = np.zeros_like(A) m = A.shape[0] #新建一个单位矩阵P P = np.identity(m) #构造一个双重循环分别取出L和U的元素 for i in range(m-1): #得到这一列中绝对值元素最大所在的行数 idx = i + np.argmax(abs(A[i:, i]))#每列比较的元素递减，因此加上了A[i:]后面的元素 #交换A，L，P的行 if idx != i & idx != m-1: A[[i, idx], :] = A[[idx, i], :] L[[i, idx], :] = L[[idx, i], :] P[[i, idx], :] = P[[idx, i], :] #主对角线元素至为1 L[i, i] = 1 for j in range(i+1, m): #计算L和U元素的值 L[j, i] = A[j, i] / A[i, i] A[j, :] = A[j, :] - L[j, i] * A[i, :] #最后处理一下主对角线元素 L[i+1, i+1] = 1 #得到U U = np.array(A) print('L矩阵为：\n', L) print('U矩阵为：\n', U) print('P矩阵为：\n', P) #测试程序 A = np.array([[1, 2, -3, 4], [4, 8, 12, -8], [2, 3, 2, 1], [-3, -1, 1, -4]]) LU_with_partial_pivoting(A)