数值分析-线性方程组的直接解法

数值积分上机课写的 Python 代码，实现线性方程组的直接解法。

import numpy as np
import time


def gaussian_elimination(A, b):
    if A.shape[0] != A.shape[1]:
        raise ValueError('A must be a square matirx', A)
    if b.shape[1] != 1:
        raise ValueError('B must be a column vector', b)
    if A.shape[0] != b.shape[0]:
        raise ValueError('A and b must have the same number of rows', A, b)
    
    n = A.shape[0]
    if np.linalg.matrix_rank(A) != n:
        raise ValueError('A is a singular matrix', A)
    
    X = np.matrix(np.empty(n)).T
    
    if A[0, 0] == 0:
        for i in range(1, len(A)):
            if A[i, 0] != 0:
                A[[0, i], :] = A[[i, 0], :]
                break
    
    for j in range(n):
        for i in range(j + 1, n):
            if A[j, j] == 0:
                for i in range(j + 1, len(A)):
                    if A[i, j] != 0:
                        A[[j, i], :] = A[[i, j], :]
                        break
            m = A[i, j] / A[j, j]
            A[i] = A[i] - A[j] * m
            b[i] = b[i] - b[j] * m
            
    X[-1] = b[-1] / A[-1, -1]
    for i in reversed(range(n - 1)):
        sum = 0
        for j in range(i + 1, n):
            sum += A[i, j] * X[j]
        X[i] = (b[i] - sum) / A[i, i]
    return X


def lu_decomposition(A, b):
    def LU(A):
        L = np.eye(len(A))
        U = np.zeros(np.shape(A))
        for r in range(1, len(A)):
            U[0, r - 1] = A[0, r - 1]
            L[r, 0] = A[r, 0] / A[0, 0]
        U[0, -1] = A[0, -1]
        for r in range(1, len(A)):
            for i in range(r, len(A)):
                delta = 0
                for k in range(0, r):
                    delta += L[r, k] * U[k, i]
                U[r, i] = A[r, i] - delta

                for i in range(r + 1, len(A)):
                    theta = 0
                    for k in range(0, r):
                        theta += L[i, k] * U[k, r]
                    L[i, r] = (A[i, r] - theta) / U[r, r]
        return L,U
    
    
    if A.shape[0] != A.shape[1]:
        raise ValueError('A must be a square matirx', A)
    if b.shape[1] != 1:
        raise ValueError('B must be a column vector', b)
    if A.shape[0] != b.shape[0]:
        raise ValueError('A and b must have the same number of rows', A, b)
    
    n = A.shape[0]
    if np.linalg.matrix_rank(A) != n:
        raise ValueError('A is a singular matrix', A)
    
    L, U = LU(A)
    n = len(A)
    y = np.zeros((n, 1))
    b = np.array(b).reshape(n, 1)
    for i in range(len(A)):
        t = 0
        for j in range(i):
            t += L[i][j] * y[j][0]
        y[i][0] = b[i][0] - t
    X = np.zeros((n, 1))
    for i in range(len(A) - 1, -1, -1):
        t = 0
        for j in range(i + 1, len(A)):
            t += U[i][j] * X[j][0]
        t = y[i][0] - t
        if t != 0 and U[i][i] == 0:
            return 0
        X[i] = t / U[i][i]
    return X