线性方程组的数值解法

*Major*-莙工科技有限公司

于 2025-02-04 10:02:55 发布

阅读量1.2k

点赞数 19

文章标签：机器学习算法人工智能

本文链接：https://blog.csdn.net/qq_41375318/article/details/145437533

版权

线性方程组的数值解法

一、直接方法

高斯消元法（Gaussian Elimination）：
高斯消元法是求解线性方程组的经典直接方法。通过对增广矩阵进行初等行变换，最终将矩阵转化为上三角矩阵，再通过回代求解。

步骤：
- 通过初等行变换将矩阵转换为上三角矩阵。
- 从最后一行开始回代，逐步求解未知数。
例子：
求解以下线性方程组：
$2 x + y + 3 z = 9$
$4 x + 2 y + 5 z = 12$
$3 x + 3 y + 2 z = 8$

步骤：
1. 构造增广矩阵：
  $\left[ \begin{array}{ccc|c} 2 & 1 & 3 & 9 \\ 4 & 2 & 5 & 12 \\ 3 & 3 & 2 & 8 \\ \end{array} \right]$
2. 通过行变换消去下三角部分，得到上三角矩阵：
  $\left[ \begin{array}{ccc|c} 2 & 1 & 3 & 9 \\ 0 & 0 & -1 & -6 \\ 0 & 0 & -3 & -8 \\ \end{array} \right]$
3. 进行回代，得到解。
LU分解（LU Decomposition）：
LU分解是将一个方阵分解为下三角矩阵 $L$ 和上三角矩阵 $U$ ，然后通过解两个三角方程来求解线性方程组。

公式：
$A = LU$
其中 $L$ 是下三角矩阵， $U$ 是上三角矩阵。

通过解 $L y = b$ 和 $Ux = y$ ，得到线性方程组的解。
矩阵求逆（Matrix Inversion）：
通过求解矩阵的逆矩阵 $A^{-1}$ ，可以得到线性方程组的解：
$x = A^{-1}b$

但该方法在数值上通常不太稳定，尤其是当矩阵的条件数较大时。

二、迭代方法

Jacobi迭代法（Jacobi Iteration）：
Jacobi迭代法是一种基于分解矩阵的迭代方法。在每次迭代中，通过使用上一轮的解来计算下一轮的解。

公式：
$x^{(k+1)}_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j\neq i} a_{ij} x^{(k)}_j \right)$
其中， $x^{(k)}_i$ 是第 $k$ 次迭代的第 $i$ 个分量， $a_{ii}$ 是矩阵 $A$ 的对角线元素。

例子：
求解以下线性方程组：
$4 x - y = 3$
$- x + 3 y = 5$

初始猜测： $x^{(0)} = 0, y^{(0)} = 0$

通过迭代更新公式，求解最终解。
Gauss-Seidel迭代法（Gauss-Seidel Iteration）：
Gauss-Seidel方法类似于Jacobi方法，但在每次迭代中，计算使用的是当前迭代的最新值，而不是上一轮的值。

公式：
$x^{(k+1)}_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j<i} a_{ij} x^{(k+1)}_j - \sum_{j>i} a_{ij} x^{(k)}_j \right)$

例子：
求解以下线性方程组：
$5 x + y = 4$
$2 x + 3 y = 5$
共轭梯度法（Conjugate Gradient Method）：
共轭梯度法是专门用于求解对称正定矩阵的迭代方法，常用于大型稀疏系统的求解。

三、条件数与解的稳定性

条件数：用于衡量线性方程组解的稳定性。条件数越大，解的误差会越大。

公式：
$\kappa(A) = \|A\| \cdot \|A^{-1}\|$
其中， $A$ 是矩阵， $A^{-1}$ 是矩阵的逆， $\| \cdot \|$ 是矩阵的范数。
解的稳定性：数值解法的稳定性与问题的条件数密切相关。如果条件数很大，使用直接方法（如高斯消元法）可能会导致误差放大，而迭代方法则可能需要更多迭代次数。

课堂活动与案例演示

课堂活动一：编程实现高斯消元法

例题：

给定以下线性方程组：
$2 x + 3 y = 5$
$4 x + y = 6$

Python代码实现：

import numpy as np

def gaussian_elimination(A, b):
    n = len(b)
    # 将增广矩阵构造为 A|b
    Augmented_matrix = np.hstack([A, b.reshape(-1, 1)])

    for i in range(n):
        # 找到最大主元进行交换
        max_row = np.argmax(np.abs(Augmented_matrix[i:n, i])) + i
        Augmented_matrix[[i, max_row]] = Augmented_matrix[[max_row, i]]

        # 对于当前行的元素进行消元
        for j in range(i + 1, n):
            factor = Augmented_matrix[j, i] / Augmented_matrix[i, i]
            Augmented_matrix[j, i:] -= factor * Augmented_matrix[i, i:]

    # 回代解方程
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Augmented_matrix[i, -1] - np.dot(Augmented_matrix[i, i + 1:n], x[i + 1:n])) / Augmented_matrix[i, i]
    
    return x

# 输入矩阵 A 和常数向量 b
A = np.array([[2, 3], [4, 1]], dtype=float)
b = np.array([5, 6], dtype=float)

# 计算解
solution = gaussian_elimination(A, b)
print("解为:", solution)

输出：

解为: [1. 1.]

课堂活动二：编程实现Jacobi迭代法

例题：

求解以下线性方程组：
$4 x - y = 3$
$- x + 3 y = 5$
初始猜测： $x^{(0)} = 0, y^{(0)} = 0$

Python代码实现：

import numpy as np

def jacobi_iteration(A, b, x0, tol=1e-6, max_iter=100):
    n = len(b)
    x = np.copy(x0)
    for k in range(max_iter):
        x_new = np.copy(x)
        for i in range(n):
            sum_ = b[i] - np.dot(A[i, :], x) + A[i, i] * x[i]
            x_new[i] = sum_ / A[i, i]
        if np.linalg.norm(x_new - x, ord=np.inf) < tol:
            break
        x = np.copy(x_new)
    return x

A = np.array([[4, -1], [-1, 3]], dtype=float)
b = np.array([3, 5], dtype=float)
x0 = np.array([0, 0], dtype=float)

solution = jacobi_iteration(A, b, x0)
print("解为:", solution)