Jacobi迭代法思想及C编程如何实现？-C语言-盈思创科技

雅可比迭代法的思想

雅可比迭代法是一种用于求解线性方程组 Ax = b 的迭代算法，它特别适用于系数矩阵 A 是大型稀疏矩阵（即矩阵中大部分元素为零）的情况。

（图片来源网络，侵删）

适用条件

该方法要求系数矩阵 A 满足以下条件之一（通常是对角占优）：

严格对角占优：矩阵的每一行，对角线元素的绝对值都大于该行其他元素绝对值之和。 |a_ii| > Σ(|a_ij|) for j ≠ i (对所有行i成立)
不可约对角占优：这是一个稍弱的条件,但对于收敛性同样重要。

如果矩阵不满足这些条件,雅可比迭代法可能不收敛。

核心思想

雅可比迭代法的基本思想是“逐步逼近”，它从一个初始猜测值出发，通过迭代公式不断更新解向量的值,直到解的精度达到要求为止。

推导过程：

（图片来源网络，侵删）

对于一个 n×n 的线性方程组 Ax = b,我们可以将其展开：

a₁₁x₁ + a₁₂x₂ + ... + a₁ₙxₙ = b₁
a₂₁x₁ + a₂₂x₂ + ... + a₂ₙxₙ = b₂
...
aₙ₁x₁ + aₙ₂x₂ + ... + aₙₙxₙ = bₙ

雅可比迭代法的核心是：在第 k+1 次迭代时，用第 k 次迭代得到的所有其他变量的值来计算当前变量的新值。

我们从第 i 个方程中解出 x_i：

a_ii * x_i = b_i - (a_i1*x₁ + ... + a_i,i-1*x_{i-1} + a_i,i+1*x_{i+1} + ... + a_in*x_n)

（图片来源网络，侵删）

然后得到迭代公式：

x_i^(k+1) = (1 / a_ii) * [ b_i - (a_i1*x₁^(k) + ... + a_i,i-1*x_{i-1}^(k) + a_i,i+1*x_{i+1}^(k) + ... + a_in*x_n^(k)) ]

重要点：

x_i^(k+1) 是 x_i 的第 k+1 次迭代值。
x_j^(k) (j ≠ i) 是 x_j 的第 k 次迭代值。
在计算 x_i 的新值时，必须使用上一轮迭代中所有其他变量的旧值，这意味着我们不能一边计算新值一边覆盖旧值，需要为解向量准备两个副本（一个存放旧值，一个存放新值）。

迭代步骤

初始化：为解向量 x 提供一个初始猜测值，通常可以设为零向量 x = [0, 0, ..., 0]ᵀ。
迭代计算：根据上述迭代公式，用 x 的旧值计算出所有分量 x_i 的新值，并存入一个新的向量 x_new 中。
更新解：将计算出的新向量 x_new 复制回解向量 x 中,为下一次迭代做准备。
检查收敛：计算新向量 x_new 和旧向量 x 之间的某种“距离”（通常是范数），||x_new - x||，如果这个距离小于一个预设的极小值（epsilon），则认为解已经收敛,迭代结束。
重复：如果未收敛，则返回步骤2,继续下一次迭代。

C语言编程实现

下面是一个完整的C语言程序,用于演示如何使用雅可比迭代法求解线性方程组。

我们将求解以下方程组作为示例：

10x₁ - x₂ + 2x₃ = 6
-x₁ + 11x₂ - x₃ + 3x₄ = 25
2x₁ - x₂ + 10x₃ - x₄ = -11
3x₂ - x₃ + 8x₄ = 15

这个方程组的系数矩阵是严格对角占优的,因此雅可比迭代法保证收敛。

C语言代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// 定义最大迭代次数，防止不收敛导致死循环
#define MAX_ITERATIONS 100
// 定义收敛精度
#define EPSILON 0.00001
// 使用雅可比迭代法求解线性方程组 Ax = b
// n: 方程组的阶数
// A: 系数矩阵 (n x n)
// b: 常数项向量 (n x 1)
// x: 解向量 (n x 1)，初始值会被使用
void jacobi_iteration(int n, double A[n][n], double b[n], double x[n]) {
    double x_new[n]; // 用于存放迭代后的新解
    int iteration;
    double error; // 用于记录误差
    for (iteration = 0; iteration < MAX_ITERATIONS; iteration++) {
        // 1. 根据旧值计算所有分量的新值
        for (int i = 0; i < n; i++) {
            double sum = 0.0;
            // sum = Σ (A[i][j] * x[j]) for j != i
            for (int j = 0; j < n; j++) {
                if (i != j) {
                    sum += A[i][j] * x[j];
                }
            }
            // 计算新解 x_i^(k+1) = (b_i - sum) / A[i][i]
            x_new[i] = (b[i] - sum) / A[i][i];
        }
        // 2. 计算误差，判断是否收敛
        error = 0.0;
        for (int i = 0; i < n; i++) {
            // 使用无穷范数 (max norm)，即最大分量的变化量
            double diff = fabs(x_new[i] - x[i]);
            if (diff > error) {
                error = diff;
            }
        }
        // 3. 更新解向量
        for (int i = 0; i < n; i++) {
            x[i] = x_new[i];
        }
        // 打印当前迭代次数和误差（可选）
        printf("Iteration %d, Error: %e\n", iteration + 1, error);
        // 4. 检查收敛条件
        if (error < EPSILON) {
            printf("\nConverged after %d iterations.\n", iteration + 1);
            return; // 成功收敛，退出函数
        }
    }
    // 如果达到最大迭代次数仍未收敛
    printf("\nWarning: Did not converge after %d iterations.\n", MAX_ITERATIONS);
}
int main() {
    // 定义方程组的阶数
    const int n = 4;
    // 1. 定义系数矩阵 A (必须是严格对角占优或不可约对角占优)
    double A[n][n] = {
        {10.0, -1.0,  2.0,  0.0},
        {-1.0, 11.0, -1.0,  3.0},
        { 2.0, -1.0, 10.0, -1.0},
        { 0.0,  3.0, -1.0,  8.0}
    };
    // 2. 定义常数项向量 b
    double b[n] = {6.0, 25.0, -11.0, 15.0};
    // 3. 定义并初始化解向量 x
    // 初始值可以设为0，或者其他合理的猜测值
    double x[n] = {0.0, 0.0, 0.0, 0.0};
    printf("Starting Jacobi iteration...\n");
    printf("Initial guess: [");
    for(int i = 0; i < n; i++) printf("%f ", x[i]);
    printf("]\n\n");
    // 4. 调用雅可比迭代法函数
    jacobi_iteration(n, A, b, x);
    // 5. 打印最终结果
    printf("\nSolution vector x:\n");
    for (int i = 0; i < n; i++) {
        printf("x[%d] = %f\n", i + 1, x[i]);
    }
    return 0;
}

代码讲解

宏定义：
- MAX_ITERATIONS：设置了一个上限,防止因为不收敛而导致程序无限循环。
- EPSILON (ε)：收敛的阈值，当解向量的最大变化量小于这个值时,我们认为解已经足够精确。
jacobi_iteration 函数：
- 参数：接收方程组的阶数 n、系数矩阵 A、常数项向量 b 和解向量 x。
- x_new 数组：这是雅可比迭代法的关键，我们用它来存储每一轮计算出的所有新解，而不会覆盖掉本轮迭代中还需要用到的旧解 x。
- 主循环 for (iteration = ...)：控制迭代次数。
- 内层循环 for (int i = 0; i < n; i++)：计算每个变量 x_i 的新值。
  - sum 变量计算了 Σ (A[i][j] * x[j]) (j ≠ i) 的部分。
  - x_new[i] = (b[i] - sum) / A[i][i] 实现了雅可比迭代的核心公式。
- 误差计算：
  - 这里我们使用了无穷范数（Infinity Norm），即 max(|x_new[i] - x[i]|)。
  - 这种方式计算简单,且能有效反映解的变化情况。
- 更新解：在计算完所有 x_new 后，一次性将 x_new 的内容复制到 x 中,为下一次迭代做准备。
- 收敛判断：如果误差 error 小于 EPSILON，打印成功信息并返回，如果循环结束仍未收敛,则打印警告信息。
main 函数：
- 定义了示例方程组的系数矩阵 A 和常数项 b。
- 初始化解向量 x 为零向量。
- 调用 jacobi_iteration 函数进行求解。
- 最后打印出最终的解向量。

预期输出

运行上述程序，你将得到类似下面的输出（具体数值可能因浮点数精度略有差异）：

Starting Jacobi iteration...
Initial guess: [0.000000 0.000000 0.000000 0.000000 ]
Iteration 1, Error: 3.125000e+00
Iteration 2, Error: 2.734375e-01
Iteration 3, Error: 8.789062e-02
Iteration 4, Error: 2.713535e-02
Iteration 5, Error: 8.271011e-03
Iteration 6, Error: 2.521286e-03
Iteration 7, Error: 7.692565e-04
Iteration 8, Error: 2.346964e-04
Iteration 9, Error: 7.157481e-05
Iteration 10, Error: 2.183079e-05
Iteration 11, Error: 6.659263e-06
Converged after 11 iterations.
Solution vector x:
x[1] = 1.000000
x[2] = 2.000000
x[3] = -1.000000
x[4] = 1.000000

可以看到，程序在11次迭代后成功收敛，并得到了精确的解 [1, 2, -1, 1]。

Jacobi迭代法思想及C编程如何实现？

雅可比迭代法的思想

适用条件

核心思想

迭代步骤

C语言编程实现

C语言代码

代码讲解

预期输出

相关文章

目录[+]