从零开始的数值分析--线性方程组直接解法-CSDN博客

本文链接：https://blog.csdn.net/qq_65711983/article/details/142619923

# 导入numpy库，用于进行科学计算
import numpy as np

# 导入scipy库，用于科学计算和工程计算
import scipy

# 导入matplotlib.pyplot，用于数据可视化
import matplotlib.pyplot as plt

# 导入sympy库，用于符号数学计算
import sympy as sp

一般来说,我们在解线性方程组是有两种办法:直接法和迭代法,在本章中,我们将介绍直接法.

在介绍Gauss消元法之前,我们先介绍一些基础的矩阵知识,我们将用代码来实现.

class matrixOfSteven:
    def __init__(self, a):
        self.array = a
        self.rows = len(a)
        self.columns = len(a[0])
        
    def __add__(self, other):
        if not isinstance(other, matrixOfSteven):
            raise TypeError("加法操作仅支持 matrixOfSteven 类型的对象")
        if self.rows != other.rows or self.columns != other.columns:
            raise ValueError("矩阵相加时，行列数必须相同")
        
        # 利用 NumPy 加法直接在原地修改，避免额外的拷贝操作
        self.array += other.array
        return self
    
    def __mul__(self, other):
        if not isinstance(other, matrixOfSteven):
            raise TypeError("乘法操作仅支持 matrixOfSteven 类型的对象")
        if self.columns != other.rows:
            raise ValueError(f"矩阵相乘时，{
     self.rows}x{
     self.columns}矩阵的列数必须等于{
     other.rows}x{
     other.columns}矩阵的行数")
        
        result_array = np.dot(self.array, other.array)
        return matrixOfSteven(result_array)
    
    def __eq__(self, other):
        if not isinstance(other, matrixOfSteven):
            return False
        if self.rows != other.rows or self.columns != other.columns:
            return False
        
        for i in range(self.rows):
            for j in range(self.columns):
                if self.array[i, j] != other.array[i, j]:
                    return False
        return True
    
    def __repr__(self) -> str:
        return str(self.array)
    
    def T(self):
        return matrixOfSteven(np.transpose(self.array))

上述例子供演示使用,在实际使用中我们仍然使用计算速度更快且功能更完善的numpy包中的矩阵来实现后续功能

Gauss消元法简介

Gauss消元法，又称高斯消元法，是一种在数学中广泛应用于解线性方程组的经典算法。该方法由著名数学家卡尔·弗里德里希·高斯命名，是线性代数中解决线性方程组问题的基础工具之一。Gauss消元法通过行变换将线性方程组转化为阶梯形矩阵（也称作行简化形式），进而求得方程组的解。下面是Gauss消元法的基本步骤和概念的Markdown格式介绍：

基本概念

线性方程组: 形如 $A x = b$ 的方程组，其中 $A$ 是系数矩阵， $x$ 是未知数向量， $b$ 是常数项向量。
增广矩阵: 将系数矩阵 $A$ 和常数项向量 $b$ 组合在一起形成的新矩阵 $[A ∣ b]$ 。
行简化: 通过交换行、倍乘行、加减行等基本行变换操作，使增广矩阵逐步转化为上三角形或阶梯形矩阵的过程。

Gauss消元法步骤

初始化:
构建增广矩阵: 将线性方程组的系数矩阵 $A$ 和常数项向量 $b$ 合并成增广矩阵 $[A ∣ b]$ 。
消元过程:
选择主元素: 对于每一列（从左至右），选择一个非零元素作为该列的“主元素”。理想情况下，选择绝对值最大的元素作为主元素，以减少舍入误差。
这一方法又称为“主元素消去法”
行变换:
1. 消零: 使用主元素所在行，通过适当的加减操作，将该列下方的所有元素变为0，这一过程称为“消元”。
2. 比例变换: 如有必要，可以通过乘以一个非零常数来调整行，使得主元素变为1，这有助于简化后续计算。
  重复步骤2和3，直至所有非主对角线位置的元素都变为0，形成上三角形矩阵。
回代求解:
回代求解: 利用上三角形矩阵，从最后一个方程开始向前，逐步解出未知数。具体来说，先解最后一个方程得到最后一个未知数的值，然后依次向前，利用已解出的未知数解出剩余未知数。

特殊情况处理

无解或无穷多解: 如果在消元过程中出现某一行全为0（除了最右边的常数项），则需根据常数项判断方程组是否有解。若所有常数项均为0，则方程组有无穷多解；若存在非零常数项，则方程组无解。
数值稳定性: 实际计算中，由于浮点运算的精度限制，应尽量避免将很小的数作为除数，以减少数值误差。

def gauss_solve(A_input, b_input, epsilon=1e-15):
    """
    高斯消元法求解方程组的解

    Args:
        A_input (np.ndarray): 输入的系数矩阵，应为方阵
        b_input (np.ndarray): 输入的常数向量,其长度应与A_input的行数或列数相等
        epsilon (float): 用于判断元素是否为0的阈值

    Returns:
        np.ndarray: 方程组的解向量

    Raises:
        ValueError: 当输入的矩阵不是方阵、输入向量的长度与矩阵的行数或列数不匹配、或者出现除零情况时抛出
    """
    
    # 校验输入类型和维度
    A = np.array(A_input).astype(float)
    b = np.array(b_input).astype(float)
    
    if not (A.ndim == 2 and b.ndim == 1):
        raise ValueError("输入的A和b的维度不正确")
    
    n = A.shape[0]
    if n != A.shape[1] or n != len(b):
        raise ValueError("输入的矩阵不是方阵或常数向量长度不匹配")
    
    # 复制矩阵以避免修改原始输入
    A = A.copy()
    b = b.copy()
    
    x = np.zeros_like(b)
    
    for i in range(n-1):
        
        # 检查主元是否接近0
        if np.abs(A[i, i]) < epsilon:
            raise ValueError(f"第{
     i+1}行第{
     i+1}列元素接近0,无法消元")
        
        # 高斯消元
        for j in range(i+1, n):
            factor = A[j, i] / A[i, i]
            A[j, :] -= factor * A[i, :]
            b[j] -= factor * b[i]
                
    # 后向 substitution
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]
        
    return x

A = np.array([[1, 2, 0], [4, 0, 6], [0, 8, 10]])
b = np.array([1, 2, 3])
x = np.linalg.solve(A, b)
x1 = gauss_solve(A, b)
print(x)
print(x1)

[0.40625  0.296875 0.0625  ]
[0.40625  0.296875 0.0625  ]

def column_main_gauss_elimination(A_in:np.array, b_in:np.array,output_key = False) -> np.array:
    
    """
    Summeration of the column-main Gauss elimination method for solving linear equations:
    采用列主元高斯消元法求解线性方程组
    
    Args:
        A_in (Matrix): 要求方阵
        b_in (Matrix): 要求长度与A_in行数相同的向量
    
    Returns:
        Matrix[]: 求解线性方程组的解向量

    """
    
    A = A_in.copy()
    
    b = b_in.copy()
    
    m = len(A)
    n = len(A[0])
    
    if m!=n:
        return "A is not a square matrix"
    
    x = np.zeros((n,1))
    
    i = 0
    while i < m-1:
        
        column_max = abs(A[i][i])
        column_max_index = i
        
        for j in range(i+1, m):
            
            if abs(A[j][i]) > abs(column_max):
                column_max = A[j][i]
                column_max_index = j
                
        if column_max == 0:
            return "A is singular"
                
        if i != j:
            for k in range(i,m):
            
                A[i][k], A[column_max_index][k] = A[column_max_index][k], A[i][k]
                
            b[i], b[column_max_index] = b[column_max_index], b[i]
            
        
        for j in range(i+1,m):
            
            param = A[j][i]/A[i][i]
            
            for k in range(i,m):
                A[j][k] = A[j][k] - param*A[i][k]
            b[j] = b[j] - param*b[i]
            
        if output_key:
            print(f"Item:{
     i+1}\nA:{
     A}\nb:{
     b}")
         
            
        i+=1
                
        
    x[n-1,0] = b[n-1] / A[n-1][n-1]    
        
    for i in range(m-2,-1,-1):
        t = 0
        for j in range(i+1,n):
            t = t + A[i][j]*x[j,0]
                
        x[i,0] = (b[i] - t) / A[i][i]
        
    return x


L2 = column_main_gauss_elimination(A, b)

print(L2)

[[ 0.5  ]
 [ 0.375]
 [-0.   ]]

LU分解介绍

LU分解是线性代数中一种重要的矩阵因式分解方法，广泛应用于数值分析领域，尤其是在解决线性方程组、计算逆矩阵和行列式值等问题上。该方法将一个满秩的方阵 $A$ 分解为一个下三角矩阵 $L$ 和一个上三角矩阵 $U$ 的乘积，即 $A = LU$ 。在某些情况下，为了保证分解过程的顺利进行，还会涉及到一个置换矩阵 $P$ ，这时分解形式变为 $P A = LU$ ，但这不影响基本思想。

基本概念

下三角矩阵(L)：

主对角线以上的元素全为0，主对角线上的元素通常为1（单位下三角矩阵），其余位置的元素可以是任意实数。
上三角矩阵(U)：

主对角线以下的元素全为0，其余位置的元素可以是任意实数。
置换矩阵(P )：

由单位矩阵经过行交换得到，作用是对矩阵的行进行重新排序，以确保在进行消元时每一步的主元不为0。

LU分解过程

初始化：设 $A$ 为 $\times n$ 的方阵，目标是找到 $L$ 和 $U$ ，使得 $A = LU$ 或 $P A = LU$ 。
消元过程：通过行变换将 $A$ 转换为上三角矩阵 $U$ ，同时记录每次行变换生成的下三角矩阵 $L$ 。具体操作包括选择主元素、消元和更新 $L$ 。
选择主元素：每一步选择当前子矩阵左上角的元素作为主元素，如果该元素为0或接近0，则可能需要进行行交换，这会引入置换矩阵 $P$ 。
消元：对于每个主元素，用它去除下面的列，确保该列下方的元素变为0，这一过程实质上是在构造 $U$ 。
更新 $L$ ：记录每次消元操作的缩放因子，形成 $L$ 的非对角线元素，而 $L$ 的对角线元素通常设为1。

应用

求解线性方程组：将方程组 $A x = b$ 转换为两个更简单的方程组 $L y = b$ 和 $Ux = y$ ，分别求解 $y$ 和 $x$ 。
求逆矩阵：一旦得到 $L$ 和 $U$ ，可以通过它们来快速计算原矩阵 $A$ 的逆矩阵。
计算行列式：矩阵 $A$ 的行列式等于 $L$ 和 $U$ 的行列式的乘积，而 $L$ 和 $U$ 的行列式又很容易计算，因为它们分别是下三角和上三角矩阵。

注意事项

LU分解只适用于方阵，并且要求矩阵为满秩，即行列式不为0。
在实际计算中，由于浮点运算的误差，直接应用高斯消元法可能不会严格得到下三角矩阵 $L$ ，因此通常采用修正的算法，如Crout分解或Doolittle算法，来确保 $L$ 保持单位下三角的形式。

def LU_decomposition(A, b=None, isCompute=False):
    """
    对矩阵A进行LU分解,如果提供了b,则同时解线性方程组Ax = b。
    
    参数:
    A: numpy.ndarray, 输入的系数矩阵，必须是方阵。
    b: numpy.ndarray, 可选,线性方程组的右侧常数向量,必须与A的列数匹配。
    isCompute: bool, 可选，指定是否计算线性方程组的解。
    
    返回:
    tuple: 包含L矩阵、U矩阵和解向量x的元组。如果isCompute为False,x为None。
    """


    # 参数检验
    # 检查A是否为方阵
    if not isinstance(A, np.ndarray) or A.ndim != 2 or A.shape[0] != A.shape[1]:
        raise ValueError("参数A必须是一个方阵。")
    # 如果提供了b，检查其维度是否与A匹配
    if b is not None and (not isinstance(b, np.ndarray) or A.shape[1] != b.shape[0]):
        raise ValueError("参数b的维度与A的列数不匹配。")
    
    n = len(A)
    
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            U[0][0] = A[0][0]
            for j in range(1, n):
                U[0][j] = A[0][j]
                L[j][0] = A[j][0] / U[0][0]
        else:
            for j in range(i, n):  # 生成U矩阵
                temp = 0
                for k in range(0, i):
                    temp += L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):  # 生成L矩阵
                temp = 0
                for k in range(0, i):
                    temp += L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]
    
    # 如果需要计算解
    # 数值解算
    if isCompute:
        
        if b is None:
            raise ValueError("当isCompute为True时,必须提供b。")
        
        Z = scipy.linalg.solve(L,b)
        
        x = scipy.linalg.solve(U,Z)
        
        return L, U, x
    
    else:
        
        return L, U, None

A = np.array([[3,1,0],[1,3,1],[0,1,3]]).astype(float)
b = np.array([6,12,9]).astype(float)
L,U,x = LU_decomposition(A,b,True)
print("\nL:\n",L,"\nU:\n",U,"\nx:\n",x)

L:
 [[1.         0.         0.        ]
 [0.33333333 1.         0.        ]
 [0.         0.375      1.        ]] 
U:
 [[3.         1.         0.        ]
 [0.         2.66666667 1.        ]
 [0.         0.         2.625     ]] 
x:
 [1. 3. 2.]