一个方程组，等式左边为一个系数矩阵（4×4）和一个列向量相乘，系数矩阵第一行为0.001，2，0，1，第二行为-3，2.099999，6，2，第三行为5，-1，5，-1，第四行为2，1，0，2，列向量为x1,x2,x3,x4,x5，等式右边为一个列向量，从上到下为8，5.900001，5，1，使用LU分解和列主元高斯消去法求解方程组，用python代码实现

首先，我们需要导入numpy库来进行矩阵运算。

import numpy as np

定义系数矩阵A和列向量B，以及初始化LU分解的L和U矩阵。

A = np.array([[0.001, 2, 0, 1], 
              [-3, 2.099999, 6, 2], 
              [5, -1, 5, -1], 
              [2, 1, 0, 2]])

B = np.array([[8], [5.900001], [5], [1]])

n = len(A)

L = np.zeros((n, n))
U = np.zeros((n, n))

接下来，我们使用列主元高斯消去法求解方程组。由于需要进行列主元高斯消去法，我们需要定义一个函数来交换矩阵的行。

def swap_rows(matrix, i, j):
    temp = np.copy(matrix[i,:])
    matrix[i,:] = matrix[j,:]
    matrix[j,:] = temp

然后，我们需要对A矩阵进行列主元高斯消去法操作。具体步骤如下：

1. 遍历每一列，找到该列中绝对值最大的元素所在的行。
2. 如果该元素所在的行不是当前列的第一行，就交换该行和当前列的第一行。
3. 对当前列进行高斯消元操作，使下面的元素都变为0。

for j in range(n):
    swap_index = j
    max_value = abs(A[j,j])

    for i in range(j+1, n):
        if abs(A[i,j]) > max_value:
            max_value = abs(A[i,j])
            swap_index = i

    if swap_index != j:
        swap_rows(A, j, swap_index)
        swap_rows(B, j, swap_index)

    for i in range(j+1, n):
        factor = A[i,j] / A[j,j]
        A[i,j:] = A[i,j:] - factor * A[j,j:]
        B[i] = B[i] - factor * B[j]

此时，A矩阵的下三角部分就是L矩阵，而A矩阵的上三角部分就是U矩阵。我们可以将其分别赋值给L和U。

for i in range(n):
    L[i,i] = 1
    for j in range(i):
        L[i,j] = A[i,j]
        U[j,i] = 0
    for j in range(i, n):
        U[i,j] = A[i,j]

接下来，我们需要解出Ly=B和Ux=y两个方程组。首先，我们解出Ly=B。

Y = np.zeros((n,1))

for i in range(n):
    Y[i] = B[i] / L[i,i]
    for j in range(i):
        Y[i] = Y[i] - L[i,j] / L[i,i] * Y[j]

然后，我们解出Ux=Y。

X = np.zeros((n,1))

for i in range(n-1, -1, -1):
    X[i] = Y[i] / U[i,i]
    for j in range(i+1, n):
        X[i] = X[i] - U[i,j] / U[i,i] * X[j]

最后，我们输出解向量X。

print(X)

完整代码如下：

import numpy as np

def swap_rows(matrix, i, j):
    temp = np.copy(matrix[i,:])
    matrix[i,:] = matrix[j,:]
    matrix[j,:] = temp

A = np.array([[0.001, 2, 0, 1], 
              [-3, 2.099999, 6, 2], 
              [5, -1, 5, -1], 
              [2, 1, 0, 2]])

B = np.array([[8], [5.900001], [5], [1]])

n = len(A)

L = np.zeros((n, n))
U = np.zeros((n, n))

for j in range(n):
    swap_index = j
    max_value = abs(A[j,j])

    for i in range(j+1, n):
        if abs(A[i,j]) > max_value:
            max_value = abs(A[i,j])
            swap_index = i

    if swap_index != j:
        swap_rows(A, j, swap_index)
        swap_rows(B, j, swap_index)

    for i in range(j+1, n):
        factor = A[i,j] / A[j,j]
        A[i,j:] = A[i,j:] - factor * A[j,j:]
        B[i] = B[i] - factor * B[j]

for i in range(n):
    L[i,i] = 1
    for j in range(i):
        L[i,j] = A[i,j]
        U[j,i] = 0
    for j in range(i, n):
        U[i,j] = A[i,j]

Y = np.zeros((n,1))

for i in range(n):
    Y[i] = B[i] / L[i,i]
    for j in range(i):
        Y[i] = Y[i] - L[i,j] / L[i,i] * Y[j]

X = np.zeros((n,1))

for i in range(n-1, -1, -1):
    X[i] = Y[i] / U[i,i]
    for j in range(i+1, n):
        X[i] = X[i] - U[i,j] / U[i,i] * X[j]

print(X)