【数值分析课程实验-03】线性代数方程组的数值解法

实验目标

使用编程语言实现以下算法：
1. 用 高斯 (Gauss) 消元法 求 n 阶线性方程组的解；
2. 用 列主元高斯 (Gauss) 消元法 求 n 阶线性方程组的解；
3. 用 列主元 LU 直接分解法 求 n 阶线性方程组的解；
4. 用 Jacobi 迭代法 求 n 阶线性方程组的解；
5. 用 Gauss-Seidel 迭代法 求 n 阶线性方程组的解。

编程语言与扩展库

编程语言：Python 语言
扩展模块：numpy

实现代码

高斯消元法 (Gauss)

# 导入模块
from numpy import *

# 高斯消元法  求解 Ax = b
def Gauss(A,b):

    # 获取 A 的行列
    row = A.shape[0]
    column = A.shape[1]
    # 非方阵情况
    if row != column:
        print('A不是方阵！')
        return

    n = row    # 未知数数量
    m = zeros((n,n))    # m 矩阵存储乘数

    # 消元过程
    for k in range(0,n-1):
        for i in range(k+1,n):
            # 第 i 行除第 k 行得乘数
            m[i,k] = A[i,k]/A[k,k]
            # 方程消元
            A[i,:] = A[i,:] - m[i,k]*A[k,:]
            b[i] = b[i] - m[i,k]*b[k]

    # 回代过程
    x = transpose(zeros((1,n)))
    x[n-1] = b[n-1]/A[n-1,n-1]    # Xn 的值

    # 计算 Xn-1, Xn-2, ... , X1 的解
    for k in range(0,n-1)[::-1]:    # 第 n-1 行到第 0 行逆序
        sum = 0
        for j in range(k+1,n):
            sum += A[k,j]*x[j]
        x[k] = (b[k]-sum)/A[k,k]

    return x


if __name__ == '__main__':
    A = array([[1,1,1],[0,4,-1],[2,-2,1]], dtype = double)
    b = transpose(array([[6,5,1]], dtype = double))
    x = Gauss(A,b)
    print(x)

列主元高斯消元法 (Gauss)

# 导入模块
from numpy import *

# 列主元高斯消元法  求解 Ax = b
def Gauss_column_pivoting(A,b):

    # 获取 A 的行列
    row = A.shape[0]
    column = A.shape[1]
    # 非方阵情况
    if row != column:
        print('A不是方阵！')
        return

    n = row    # 未知数数量
    m = zeros((n, n))    # m 矩阵存储乘数

    for k in range(0,n-1):
    # 寻找最大主元的行数
        h = argmax(abs(A[:,k]))
        # 交换行
        if k != h:
            temp1 = copy(A[k,:])    # copy 产生该行的复制
            A[k,:] = A[h,:]
            A[h,:] = temp1
            temp2 = copy(b[k])
            b[k] = b[h]
            b[h] = temp2

        # 消元过程
        for i in range(k + 1, n):
            # 第 i 行除第 k 行得乘数
            m[i, k] = A[i, k] / A[k, k]
            # 方程消元
            A[i, :] = A[i, :] - m[i, k] * A[k, :]
            b[i] = b[i] - m[i, k] * b[k]

    # 回代过程 
    x = transpose(zeros((1, n)))    # 构造 x 列向量
    x[n - 1] = b[n - 1] / A[n - 1, n - 1]    # Xn 的值

    # 计算 Xn-1, Xn-2, ... , X1 的解
    for k in range(0, n - 1)[::-1]:    # 第 n-1 行到第 0 行逆序
        sum = 0
        for j in range(k + 1, n):
            sum += A[k, j] * x[j]
        x[k] = (b[k] - sum) / A[k, k]

    return x

if __name__ == '__main__':
    A = array([[1, 1, 1], [0, 4, -1], [2, -2, 1]],  dtype = double)
    b = transpose(array([[6, 5, 1]], dtype = double))
    x = Gauss_column_pivoting(A,b)
    print(x)

列主元 LU 直接分解法

# 导入模块
from numpy import *

# 列主元 LU 分解法  求解 Ax = b
def LU_column_pivoting(A,b):

    # 获取 A 的行列
    row = A.shape[0]
    column = A.shape[1]
    # 非方阵情况
    if row != column:
        print('A不是方阵！')
        return

    n = row
    Ip = list(range(n))    # 标记主元交换次序

    for k in range(0,n):

        # 寻找最大主元的行数
        h = argmax(abs(A[:,k]))

        # 只能选择未交换的行进行交换
        if Ip[h] == h:
             temp = copy(A[k, :])
             A[k, :] = A[h, :]
             A[h, :] = temp
             Ip[k] = h

        # 对角线减去对应行列乘积和
        A[k,k] = A[k,k] - dot(A[k,0:k], A[0:k,k])
        if A[k,k] == 0:
            print('第{}个主元为0！'.format(k+1))
            return

        # 直接在 A 上保存 L 和 U
        for j in range(k+1,n):    # j > k
            A[k,j] = A[k,j] - dot(A[k,0:k], A[0:k,j])
            A[j,k] = (A[j,k] - dot(A[j,0:k], A[0:k,k])) / A[k,k]

    # 单位下三角 L 和上三角 U
    L = tril(A,-1) + identity((n))
    U = triu(A,0)
    L.dtype = double
    U.dtype = double

    # 交换 b 对应行
    for i in range(0,n):
        t = Ip[i]
        if t != i:
            temp = copy(b[i])
            b[i] = b[t]
            b[t] = temp

    # 求解 Ly = b
    y = transpose(array([arange(n)], dtype = double))
    y[0] = b[0]
    for i in range(1,n):
        sum = 0
        for t in range(0,i):
            sum += L[i,t]*y[t]
        y[i] = b[i] - sum

    # 求解 Ux = y
    x = transpose(array([arange(n)],dtype = double))
    x[n-1] = y[n-1]/U[n-1,n-1]
    for i in range(0,n-1)[::-1]:
        sum = 0
        for t in range(i+1,n):
            sum +=  U[i,t]*x[t]
        x[i] = (y[i]-sum)/U[i,i]

    return x

if __name__ == '__main__':
    A = array([[1,2,3],[2,5,2],[3,1,5]], dtype = double)
    b = transpose(array([[14,18,20]], dtype = double))
    x = LU_column_pivoting(A,b)
    print(x)

Jacobi 迭代法

# 导入模块
from numpy import *

# 初始解向量，迭代上限，误差范围
x_start = array([[0],[0],[0]])
N = 100
tol = 1e-7

# Jacobi 迭代
def Jacobi(A,b):

    # 获取 A 的行列
    row = A.shape[0]
    column = A.shape[1]
    # 非方阵情况
    if row != column:
        print('A不是方阵！')
        return

    n = row
    x_k = x_start

    # 删去 A 的对角线
    A_ = copy(A)
    for i in range(n):
        A_[i,i] = 0

    # 迭代求解
    for i in range(N):

        t = transpose([diag(A)])
        x_k_plus = (b - dot(A_, x_k))/t

        # 达到足够小的误差
        if (max(abs((x_k - x_k_plus)[:,0])) < tol):
            print('迭代次数: {}'.format(i+1))
            print('方程组解: ')
            print(x_k_plus)
            break
        x_k = x_k_plus

    if  i == N:
        print('错误！不收敛！')

if __name__ == '__main__':
    A = array([[8,-3,2],[4,11,-1],[6,3,12]], dtype = double)
    b = transpose(array([[20,33,36]], dtype = double))
    Jacobi(A,b)

Gauss-Seidel 迭代法

# 导入模块
from numpy import *

# 初始解向量，迭代上限，误差范围
x_start = array([[1],[2],[3]])
N = 1000
tol = 1e-10

# GS 迭代法  求解 Ax = b
def Gauss_Seidel(A,b):

    # 获取 A 的行列
    row = A.shape[0]
    column = A.shape[1]
    # 非方阵情况
    if row != column:
        print('A不是方阵！')
        return

    # 初值
    n = row
    x_k = x_start

    # 构造 DLU
    U = copy(A)    # U
    DL = copy(A)    # D-L
    for i in range(n):
        for j in range(n):
            if j<=i:
                U[i,j] = 0
            else:
                DL[i,j] = 0
    U = 0-U

    # DL 逆与 U,b 相乘
    DL_U = dot(linalg.inv(DL),U)
    DL_b = dot(linalg.inv(DL),b)

    # 迭代求解
    for i in range(N):
        x_k_plus = dot(DL_U,x_k) + DL_b

        # 达到误差要求
        if max(abs((x_k_plus - x_k)[:,0])) < tol:
            print('方程组解: ')
            print(x_k_plus)
            break
        x_k = x_k_plus

    if i == N:
        print('错误！不收敛！')

if __name__ == '__main__':
    A = array([[8, -3, 2], [4, 11, -1], [6, 3, 12]], dtype=double)
    b = transpose(array([[20, 33, 36]], dtype=double))
    Gauss_Seidel(A,b)

写在最后

声明：本专栏内容来源于武汉理工大学 2019-2020 学年数值分析方法课程实验，仅供学习参考。
如有错误或不足地方，还请指出，祝愿读者能够在编程之路上不断进步！