计算方法作业，在已给matlab 程序基础上，进行修改得到的python程序。
1.高斯消去法：
import numpy as np#注意本程序只对A矩阵为n*n的有效A = np.array([[10,-7,0,1],[-3,2.099999,6,-2],[5,-1,5,-1],[2,1,0,2]])b = np.array([[8],[5.900001],[5],[1]])n = A.shape[0]solution = np.zeros(n)#消元过程def gaosixiaoqu(A,b):for k in range(n-1):#遍历行,从首行开始到倒数第一行for i in range(k+1,n):#定义每一行的消元因子mik=A[i,k]/A[k,k]#每一行消元因子的计算方法for j in range(k+1,n):A[i,j]= A[i,j]-mik*A[k,j]#遍历每一列，注意，这一步相减的时候前面默认减为0的哪一个不参与运算b[i]=b[i]-mik*b[k]#从第二行开始算return A,b#回代过程def huidai(A,b):solution[n-1]=b[n-1]/A[n-1,n-1]#最后一行直接求for i in range(n-2,-1,-1):#从倒数第二行开始for j in range(i+1,n):solution[i] = solution[i]+A[i,j]*solution[j]solution[i]=(b[i]-solution[i])/A[i,i]return solutionA1,b1 = gaosixiaoqu(A,b)solution1 = huidai(A1,b1)print('高斯消去法结果',solution1)
2.列主元高斯消去法：
import numpy as np#注意本程序只对矩阵为n*n的有效A = np.array([[10,-7,0,1],[-3,2.099999,6,-2],[5,-1,5,-1],[2,1,0,2]],dtype=float)b = np.array([[8],[5.900001],[5],[1]],dtype=float)n = A.shape[0]solution = np.zeros(n)def find_max(B,k):#此时B传入的为一个行向量row = abs(B).argmax(axis=0)+k#为了找到其实际所在的行column = kreturn row,column#返回行和列for k in range(n-1):#从第一行开始到倒数第二行value=max(abs(A[k:n,k]))#主元的值row,column =find_max(abs(A[k:n,k]),k)if row!=k: #若A(k,k)不是绝对值最大的，换位置a_k_position=A[k,k:n].copy()#当前这一行的所有值,一定要用copy()!否则互换会不成功的b_k_position=b[k].copy()#同上A[k,k:n]=A[row,k:n]#把最大值所在的一行换上去A[row,k:n]=a_k_position#把原那一行换下来b[k]=b[row]#同上b[row]=b_k_position#换完一行就得消元for i in range(k + 1, n): # 定义每一行的消元因子mik = A[i, k] / A[k, k]# 每一行消元因子的计算方法for j in range(k + 1, n):t = A[i, j].copy() - mik * A[k, j]A[i, j] = t # 遍历每一列，注意，这一步相减的时候前面默认减为0的哪一个不参与运算b[i] = b[i] - mik * b[k] # 从第二行开始算def huidai(A,b):solution[n-1]=b[n-1]/A[n-1,n-1]#最后一行直接求for i in range(n-2,-1,-1):#从倒数第二行开始for j in range(i+1,n):solution[i] = solution[i]+A[i,j]*solution[j]solution[i]=(b[i]-solution[i])/A[i,i]return solutionsolution1 = huidai(A,b)print('列主元高斯消去法结果：',solution1)
3.直接三角形分解法（LU分解）
import numpy as np#注意本程序只对矩阵为n*n的有效A = np.array([[10,-7,0,1],[-3,2.099999,6,-2],[5,-1,5,-1],[2,1,0,2]],dtype=float)b = np.array([[8],[5.900001],[5],[1]],dtype=float)n = A.shape[0]y = np.zeros(n)x = np.zeros(n)U = np.zeros((n,n),dtype=np.float64)L = np.zeros((n,n),dtype=np.float64)U[0,:]=A[0,:]#第一行赋值L[1:n,0]=A[1:n,0]/U[0,0]for k in range(1,n):for j in range(k,n):t = np.matmul(L[k,0:k],U[0:k,j])#一定要用这个乘法U[k,j]=A[k,j]-np.matmul(L[k,0:k],U[0:k,j])for i in range(k+1,n):L[i,k]=(A[i,k]-np.matmul(L[i,0:k],U[0:k,k]))/U[k,k]#求yy[0]=b[0]for kk in range(1,n):t2=np.matmul(L[kk,0:kk],y[0:kk])y[kk]=b[kk]-np.matmul(L[kk,0:kk],y[0:kk])#求xx[n-1]=y[n-1]/U[n-1,n-1]for i in range(n-2,-1,-1):x[i]=(y[i]-np.matmul(U[i,i+1:n],x[i+1:n]))/U[i,i]print('LU分解结果',x)
