Python实现列主元高斯消去法与LU分解法

数值分析：Python实现列主元高斯消去法与LU分解法求解线性方程组

一、矩阵形式的线性代数方程组

二、高斯消去法

三、高斯列主元消去法

四、矩阵三角分解法（LU分解）

这里只简单介绍Doolittle分解法。

题目：编写列主元高斯消去法与LU分解法解线性方程组Ax=b。

列主元消去法代码实现：

import mathimport numpy as np#目的：熟悉列主元消去法，以及三角分解法等直接求解线性方程组的算法#列主元消元法def CME(a,b,x):isdet0 = 0m, n = a.shape #矩阵a的行数和列数# j表示列for k in range(n - 1):# k表示第一层循环，(0，n-1)行#在每次计算前，找到最大主元，进行换行ans = np.fabs(a[k][k])ik = kfor i in range(k+1, n):if ans < np.fabs(a[i][k]): # fabs是绝对值，将a中绝对值最大的找出来ik = ians = np.fabs(a[i][k])if np.fabs(ans) < 1e-10:isdet0 = 1breakif ik != k :for i in range(k,m):temp = a[k][i]a[k][i] = a[ik][i]a[ik][i] = temptemp = b[k]b[k] = b[ik]b[ik] = tempfor i in range(k + 1, n): # i表示第二层循环,(k+1,n)行,计算该行消元的系数temp = a[i][k] / a[k][k]#计算for j in range(k,m):# j表示列，对每一列进行运算a[i][j] = a[i][j] - temp * a[k][j]b[i] = b[i] - temp * b[k]# 回代求出方程解if np.fabs(a[n-1][n-1]) < 1e-10 :isdet0 = 1if isdet0 == 0:# x = np.zeros(n)x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解for i in range(n - 2, -1, -1): #依次回代倒着算每一个解temp = 0for j in range(n - 1, i,-1):temp = temp + a[i][j]*x[j]x[i] = (b[i]-temp) / a[i][i]for i in range(n):print("x" + str(i + 1) + " = ", x[i])print("x" " = ", x)if __name__ == '__main__':#当模块被直接运行时，以下代码块将被运行，当模块是被导入时，代码块不被运行。a = np.array([[3.01, 6.03, 1.99], [1.27, 4.16, -1.23], [0.987, -4.81, 9.34]])b = np.array([1.0, 1.0, 1.0])m,n = a.shapex = np.zeros(n)B = np.zeros((n, n))for i in range(n):for j in range(n):B[i][j] = a[i][j]CME(a,b,x)#验证for i in range(0, n):temp = 0for j in range(0, n):temp = temp + B[i][j] * x[j]print("%f ", temp)if __name__ == '__main__':main()

LU分解法代码实现：

import mathimport numpy as np#目的：熟悉列主元消去法，以及三角分解法等直接求解线性方程组的算法#列主元消元法def CME(a,b,x):isdet0 = 0m, n = a.shape #矩阵a的行数和列数# j表示列for k in range(n - 1):# k表示第一层循环，(0，n-1)行#在每次计算前，找到最大主元，进行换行ans = np.fabs(a[k][k])ik = kfor i in range(k+1, n):if ans < np.fabs(a[i][k]): # fabs是绝对值，将a中绝对值最大的找出来ik = ians = np.fabs(a[i][k])if np.fabs(ans) < 1e-10:isdet0 = 1breakif ik != k :for i in range(k,m):temp = a[k][i]a[k][i] = a[ik][i]a[ik][i] = temptemp = b[k]b[k] = b[ik]b[ik] = tempfor i in range(k + 1, n): # i表示第二层循环,(k+1,n)行,计算该行消元的系数temp = a[i][k] / a[k][k]#计算for j in range(k,m):# j表示列，对每一列进行运算a[i][j] = a[i][j] - temp * a[k][j]b[i] = b[i] - temp * b[k]# 回代求出方程解if np.fabs(a[n-1][n-1]) < 1e-10 :isdet0 = 1if isdet0 == 0:# x = np.zeros(n)x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解for i in range(n - 2, -1, -1):#依次回代倒着算每一个解temp = 0for j in range(n - 1, i,-1):temp = temp + a[i][j]*x[j]x[i] = (b[i]-temp) / a[i][i]for i in range(n):print("x" + str(i + 1) + " = ", x[i])print("x" " = ", x)#三角消元法def LU(a,b,x):m, n = a.shape # 矩阵a的行数和列数y = np.array([0.0, 0.0, 0.0])for j in range(1,n):# L的第0列a[j][0] = a[j][0] / a[0][0]for i in range(1,n-1):# 求U的第i行 L的第i行for j in range(i,n):#求U的第i行的第j个元素sum = 0.0 #求和for s in range(0,i):sum = sum +a[i][s] * a[s][j]a[i][j] = a[i][j] - sum#求L的第i列的第j个元素 在j行i列for j in range(i+1,n):sum = 0.0for s in range(0,i):sum = sum + a[j][s] * a[s][i]a[j][i] = ( a[j][i] - sum ) / a[i][i]#求U[n-1][n-1]sum = 0.0 #求和for s in range(0,n-1):sum = sum + a[n-1][s] * a[s][n-1]a[n-1][n-1] = a[n-1][n-1] - sumy[0] = b[0]for i in range(1,n):sum = 0.0for j in range(0,i):sum = sum + a[i][j] * y[j]y[i] = b[i] - sumx[n-1] = y[n-1] / a[n-1][n-1]for i in range(n-2,-1,-1):#求x[i]sum = 0.0for j in range(n-1,i,-1):sum = sum + a[i][j] * x[j]x[i] = ( y[i] - sum ) / a[i][i]for i in range(n):print("x" + str(i + 1) + " = ", x[i])print("x" " = ", x)if __name__ == '__main__':#当模块被直接运行时，以下代码块将被运行，当模块是被导入时，代码块不被运行。a = np.array([[3.01, 6.03, 1.99], [1.27, 4.16, -1.23], [0.987, -4.81, 9.34]])b = np.array([1.0, 1.0, 1.0])m,n = a.shapex = np.zeros(n)B = np.zeros((n, n))for i in range(n):for j in range(n):B[i][j] = a[i][j]# CME(a,b,x)LU(a,b,x)#验证for i in range(0, n):temp = 0for j in range(0, n):temp = temp + B[i][j] * x[j]print("%f ", temp)