高斯消去法，列主元法，LU分解法python程序

计算方法作业，在已给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 solution
A1,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 = k
    return 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 solution
solution1 = 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]

#求y
y[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])
#求x
x[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)