import numpy as np
'''
列主元高斯消元法
A:系数增广矩阵
n:未知数个数
'''defmain_element_gauss(A,n):for i in range(0,n-1):
if(np.max(A[i:,i])!=A[i,i]): #如果当前系数不是最大值,则列主元
temp_i=int(np.where(A==np.max(A[i:,i]))[0]) #temp_i为最大值所在的行索引
A[[i,temp_i],:]=A[[temp_i,i],:] #交换for j in range(1,(n-i)):
A[i+j,:]=A[i+j,:]-A[i,:]*(A[i+j,i]/A[i,i]) #消元
print("第%d次消元系数矩阵为:\n"%(i+1),A)
#回代得解
X=np.zeros((n,1))
for i in range(n-1,-1,-1):
temp=0for j in range(0,n):
temp+=A[i,j]*X[j,0]
X[i,0]=(A[i,n]-temp)/A[i,i]
print("方程组的解为:\n",X)
'''
LU三角分解法
A:系数增广矩阵
'''defLU_break_down(A):#生成维度为(3,3)的矩阵,初值为0
L=np.zeros((3,3))
U=np.zeros((3,3))
#根据公式得到L和U的具体值
U[0,:]=A[0,:3]
for i in range(3):
L[i,i]=1.0
L[1,0]=A[1,0]/U[0,0]
L[2,0]=A[2,0]/U[0,0]
U[1,1]=A[1,1]-L[1,0]*U[0,1]
U[1,2]=A[1,2]-L[1,0]*U[0,2]
L[2,1]=(A[2,1]-L[2,0]*U[0,1])/U[1,1]
U[2,2]=A[2,2]-(L[2,0]*U[0,2]+L[2,1]*U[1,2])
print("L为:\n",L)
print("U为:\n",U)
#将b添加到L,使之变为增广矩阵
L=np.hstack((L,A[:,3].reshape(3,1)))
#回代得解
Y=np.zeros((3,1))
for i in range(0,3):
temp=0for j in range(0,3):
temp+=L[i,j]*Y[j,0]
Y[i,0]=(L[i,3]-temp)/L[i,i]
#将Y添加到U,使之变为增广矩阵
U=np.hstack((U,Y))
#回代得解
X=np.zeros((3,1))
for i in range(2,-1,-1):
temp=0for j in range(0,3):
temp+=U[i,j]*X[j,0]
X[i,0]=(U[i,3]-temp)/U[i,i]
print("方程组的解为:\n",X)
defmain():#输入
n=int(input("请输入未知数个数:"))
temp=[]
print("请输入系数的增广矩阵:\n")
for i in range(n):
temp.append([float(i) for i in input().split()])
A=np.array(temp).reshape(n,n+1)
#只有n=3时,有两种方法。因为LU没有支持除3元以外的解法if(n==3):
#菜单选择
print("解方程组方法:\n\t1.列主元高斯消元法\n\t2.LU三角分解法")
choice=int(input("请选择一种:"))
if(choice==1):
main_element_gauss(A,n)
else:
LU_break_down(A)
else:
main_element_gauss(A,n)
if __name__=='__main__':
main()
一、实验内容二、代码(python)import numpy as np''' 列主元高斯消元法 A:系数增广矩阵 n:未知数个数'''def main_element_gauss(A,n): for i in range(0,n-1): if(np.max(A[i:,i])!=A[i,i]): #如果当前系数不是最大值,则列主元