左同学在这
Posts Tags Categories About
左同学在这

Python数值计算：方程组求解算法实现

ZuoPeng  included in 编程实践
   1064 words   5 minutes 
Contents

理论推导请参考《数值分析》李庆扬

完整代码放置在GitHub

线性方程组的求解是数值计算中的基础问题。本文将介绍并实现几种常用的数值方法，包括直接方法和迭代方法。每种方法都有其特点和适用场景。

主要方法分类：

  • 直接方法：Gauss消去法、列主元素消去法、追赶法
  • 迭代方法：Jacobi迭代法、Gauss-Seidel迭代法

Gauss消去法是求解线性方程组最基本的直接方法，通过消元过程将系数矩阵化为上三角矩阵，然后回代求解。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

import numpy as np

# 测试矩阵
A = np.array([[0.4096,0.1234,0.3678,0.2943,0.4043],
              [0.2246,0.3872,0.4015,0.1129,0.1550],
              [0.3645,0.1920,0.3781,0.0643,0.4240],
              [0.1784,0.4002,0.2786,0.3927,-0.2557]])

print("原始增广矩阵:")
print(A)

rows, columns = A.shape
ans = np.zeros(shape=(rows,1))

# 前向消元过程
for c1 in range(rows-1):
    for c2 in range(rows-c1-1):
        if A[c1+c2+1, c1] != 0:
            # 计算消元因子
            a = A[c1+c2+1, c1] / A[c1, c1]
            # 消元操作
            A[c1+c2+1:c1+c2+2, :] = A[c1+c2+1:c1+c2+2, :] - A[c1:c1+1, :] * a

print("\\n消元后的上三角矩阵:")
print(A)

# 回代求解过程
for i in range(rows):
    ans[rows-i-1] = A[rows-i-1, columns-1] / A[rows-i-1, columns-2-i]
    A[:, columns-1:columns] = A[:, columns-1:columns] - ans[rows-i-1] * A[:, rows-i-1:rows-i]

print("\\n解向量:")
print(ans)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

原始增广矩阵:
[[ 0.4096  0.1234  0.3678  0.2943  0.4043]
 [ 0.2246  0.3872  0.4015  0.1129  0.155 ]
 [ 0.3645  0.192   0.3781  0.0643  0.424 ]
 [ 0.1784  0.4002  0.2786  0.3927 -0.2557]]

消元后的上三角矩阵:
[[ 4.09600000e-01  1.23400000e-01  3.67800000e-01  2.94300000e-01  4.04300000e-01]
 [-2.77555756e-17  3.19534863e-01  1.99820605e-01 -4.84764160e-02 -6.66937988e-02]
 [ 7.13898500e-18  0.00000000e+00 -5.98156607e-04 -1.85126813e-01  8.13706787e-02]
 [-1.14249713e-15  0.00000000e+00  1.38777878e-17  3.07244855e+01 -1.37247573e+01]]

解向量:
[[-0.18191778]
 [-1.66303081]
 [ 2.21722832]
 [-0.44670422]]
  • 优点：简单直接，适用于大多数情况
  • 缺点：对主元素敏感，可能导致数值不稳定
  • 时间复杂度：O(n³)

为了改善Gauss消去法的数值稳定性，引入列主元素选取策略，在每一步消元前选择绝对值最大的元素作为主元。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import numpy as np

A = np.array([[0.4096,0.1234,0.3678,0.2943,0.4043],
              [0.2246,0.3872,0.4015,0.1129,0.1550],
              [0.3645,0.1920,0.3781,0.0643,0.4240],
              [0.1784,0.4002,0.2786,0.3927,-0.2557]])

np.set_printoptions(precision=9)
print("原始增广矩阵:")
print(A)

rows, columns = A.shape
ans = np.zeros(shape=(rows,1))

# 前向消元过程（含列主元选取）
for c1 in range(rows-1):
    # 寻找列主元
    k = c1
    for j in range(c1+1, rows):
        if abs(A[k,c1]) < abs(A[j,c1]):
            k = j
    
    # 交换行
    if k != c1:
        A[[c1, k]] = A[[k, c1]]
    
    # 消元操作
    for c2 in range(rows-c1-1):
        if A[c1+c2+1, c1] != 0:
            a = A[c1+c2+1, c1] / A[c1, c1]
            A[c1+c2+1] = A[c1+c2+1] - A[c1] * a

# 回代求解
for i in range(rows):
    ans[rows-i-1] = A[rows-i-1, columns-1] / A[rows-i-1, columns-2-i]
    A[:, columns-1:columns] = A[:, columns-1:columns] - ans[rows-i-1] * A[:, rows-i-1:rows-i]

print("\\n解向量:")
print(ans)

1
2
3
4
5

解向量:
[[-0.181917777]
 [-1.663030813]
 [ 2.217228318]
 [-0.446704219]]
  • 数值稳定性：通过选择较大的主元，减少舍入误差
  • 适用性：适合处理系数矩阵条件数较大的情况
  • 可靠性：在工程计算中广泛应用

追赶法专门用于求解三对角线性方程组，是一种高效的直接方法。通过LU分解将三对角矩阵分解为下三角和上三角矩阵。

对于三对角矩阵： [b1c1a2b2c2a3b3c3⋱⋱]\begin{bmatrix} b_1 & c_1 & & \\ a_2 & b_2 & c_2 & \\ & a_3 & b_3 & c_3 \\ & & \ddots & \ddots \end{bmatrix}

进行LU分解，然后通过前代和回代求解。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

import numpy as np

# 三对角矩阵的LU分解
n = 10  # 矩阵维数
off0_i = np.arange(n)     # 主对角线索引
off0_j = np.arange(n)
off1_i = np.arange(n-1)   # 上对角线索引
off1_j = np.arange(1, n)
offd1_i = np.arange(1, n) # 下对角线索引
offd1_j = np.arange(n-1)

# 矩阵分解
L = np.zeros(shape=(n,n))
U = np.eye(n)
b = A[off0_i, off0_j].reshape((n,1))    # 主对角线元素
c = A[off1_i, off1_j].reshape((n-1,1))  # 上对角线元素
gamma = A[offd1_i, offd1_j].reshape((n-1,1))  # 下对角线元素

alpha = np.zeros(shape=(n, 1))
beta = np.zeros(shape=(n-1, 1))

# 设置L矩阵的下对角线
L[offd1_i, offd1_j] = A[offd1_i, offd1_j]

# 计算分解系数
alpha[0, 0] = b[0, 0]
beta[0, 0] = c[0, 0] / alpha[0, 0]

for i in range(1, n):
    alpha[i, 0] = b[i, 0] - gamma[i-1, 0] * c[i-1, 0] / alpha[i-1, 0]
    if i != n-1:
        beta[i, 0] = c[i, 0] / alpha[i, 0]

# 构造L和U矩阵
L[off0_i, off0_j] = alpha.reshape((n,))
U[off1_i, off1_j] = beta.reshape((n-1,))

# 验证分解的正确性
print("L*U矩阵:")
print(np.dot(L,U))

# 前代求解Ly=b
y = np.zeros(shape=(n,1))
x = np.zeros(shape=(n,1))
y[0,0] = B[0,0] / alpha[0,0]
for i in range(1,n):
    y[i,0] = (B[i,0] - gamma[i-1,0] * y[i-1,0]) / alpha[i,0]

# 回代求解Ux=y
x[n-1, 0] = y[n-1, 0]
for i in range(n-2, -1, -1):
    x[i,0] = y[i,0] - beta[i,0] * x[i+1,0]

print("解向量:")
print(x.reshape((n,-1)))
  • 高效性：时间复杂度为O(n)，远优于一般方法的O(n³)
  • 专用性：专门针对三对角矩阵设计
  • 应用场景：差分方程、样条插值等问题

Jacobi迭代法是一种经典的迭代方法，通过不断迭代逼近真解。适用于大型稀疏矩阵方程组的求解。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

import numpy as np

# 测试方程组
AB = np.array([[8,-3,2,20],[4,11,-1,33],[6,3,12,36]])
np.set_printoptions(precision=12)

rows, columns = AB.shape
A = AB[:, 0:rows]  # 系数矩阵
B = AB[:, rows:columns]  # 常数向量

ans1 = np.zeros(shape=(rows,1))
ans2 = np.zeros(shape=(rows,1))

# 提取对角线元素
a = np.diagonal(A).reshape((rows, 1))

# 构造迭代矩阵
A_new = -A
i, j = np.diag_indices_from(A_new)
A_new[i,j] = 0  # 对角线置零
A_new = A_new / a  # 归一化
B = B / a

# 迭代求解
err = 1
iteration_count = 0
while err >= 1.e-10:
    ans2 = np.dot(A_new, ans1) + B
    err = np.max(np.abs(ans1 - ans2))
    ans1 = ans2.copy()
    iteration_count += 1

print(f"经过{iteration_count}次迭代后的解:")
print(ans2)

1
2
3
4

经过XX次迭代后的解:
[[3.000000000013]
 [2.000000000013]
 [0.999999999992]]

Jacobi迭代法收敛的充分条件：

  • 对角占优：∣aii∣>∑j≠i∣aij∣|a_{ii}| > \sum_{j≠i}|a_{ij}|（严格对角占优）
  • 正定矩阵：系数矩阵为对称正定矩阵

Gauss-Seidel迭代法是Jacobi迭代法的改进版本，使用最新计算出的分量值进行迭代，通常具有更好的收敛性。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

import numpy as np

def gauss_seidel_iteration_matrix(A):
    """计算Gauss-Seidel迭代矩阵"""
    n = len(A)
    D = np.diag(np.diag(A))      # 对角矩阵
    L = np.tril(A, k=-1)         # 下三角矩阵
    U = np.triu(A, k=1)          # 上三角矩阵
    D_inv = np.linalg.inv(D)
    M = np.dot(D_inv, L + U)
    return M

def spectral_radius(matrix):
    """计算谱半径"""
    eigenvalues, _ = np.linalg.eig(matrix)
    return max(abs(eigenvalues))

def analyze_convergence(A):
    """分析收敛性"""
    iteration_matrix = gauss_seidel_iteration_matrix(A)
    radius = spectral_radius(iteration_matrix)
    
    print("Gauss-Seidel 迭代矩阵:")
    print(iteration_matrix)
    print(f"谱半径: {radius:.6f}")
    
    if radius < 1:
        print("收敛条件满足：谱半径 < 1")
    else:
        print("警告：可能不收敛，谱半径 >= 1")

# 收敛性分析
analyze_convergence(A)

# 迭代求解实现
ans1 = np.zeros(shape=(n,1))
ans2 = np.zeros(shape=(n,1))

# 处理系数矩阵
a = np.diagonal(A).reshape((n, 1))
A_iter = -A
i, j = np.diag_indices_from(A_iter)
A_iter[i,j] = 0
A_iter = A_iter / a
B = B / a

err = 1
k = 0
while err >= 1.e-5:
    ans_temp = ans1.copy()
    # 使用最新的分量值进行迭代
    for i in range(n):
        ans2[i,0] = np.dot(A_iter[i:i+1], ans1) + B[i,0]
        ans1[i,0] = ans2[i,0]  # 立即使用新值
    
    err = np.max(np.abs(ans_temp - ans2))
    k += 1

print(f"经过{k}次迭代后的解:")
x = ans2.reshape((n,-1))
print(x)
  • 收敛速度：通常比Jacobi迭代法收敛更快
  • 内存效率：只需要一个解向量存储空间
  • 适用性：特别适合大型稀疏矩阵

通过计算迭代矩阵的谱半径来判断收敛性：

  • 谱半径 < 1：保证收敛
  • 谱半径 ≥ 1：可能不收敛
方法类型时间复杂度适用场景主要优点主要缺点
Gauss消去法直接法O(n³)一般线性方程组简单直接数值不稳定
列主元消去法直接法O(n³)一般线性方程组数值稳定计算量大
追赶法直接法O(n)三对角矩阵极高效率应用受限
Jacobi迭代法迭代法O(kn²)大型稀疏矩阵并行性好收敛慢
Gauss-Seidel迭代法迭代法O(kn²)大型稀疏矩阵收敛较快难并行化
  1. 小型密集矩阵：使用列主元Gauss消去法
  2. 大型稀疏矩阵：优先考虑迭代方法
  3. 三对角矩阵：首选追赶法
  4. 对角占优矩阵：迭代方法效果好
  5. 病态矩阵：需要特殊处理或正则化

通过合理选择算法，可以在保证精度的同时提高计算效率。在实际应用中，还需要考虑矩阵的条件数、稀疏性、对称性等特征来选择最适合的求解方法。

Updated on 2025-09-11
Read Markdown
 Python, 数值计算, 线性代数, 算法实现, 数值分析
Back | Home