解线性方程组迭代法+p190+1+8+C语言版

/*
 * Solving a system of linear equations using iterative methods:
 *   1. Jacobi Iterative Method (雅可比迭代法)
 *   2. Gauss-Seidel Iterative Method (高斯 - 赛德尔迭代法)
 *   3. SOR (Successive Over-Relaxation) Method (逐次超松弛迭代法)
 *
 * Target System (系数矩阵 A 需满足对角占优以保证收敛):
 *   5x1 +  2x2 + x3 =  -12
 *   -x1 +  4x2 +  2x3  =  20
 *   2x1 -  3x2 + 10x3 = 3
 *
 * Theoretical Exact Solution: x1 = -4, x2 = 3, x3 = 2
 */

#include <stdio.h>
#include <math.h>

#define N       3           /* 未知数个数 (number of unknowns) */
#define MAX_IT  1000        /* 最大迭代次数，防止不收敛时死循环 */
#define TOL     1e-8        /* 收敛容差 (convergence tolerance)，使用无穷范数判断 */

/* ------------------------------------------------------------------ */
/*  Print a solution vector
 *  功能：格式化打印解向量 x 的各个分量
 * ------------------------------------------------------------------ */
static void print_solution(const double x[N])
{
    for (int i = 0; i < N; i++)
        // 输出格式：x1 = value, 保留 8 位小数
        printf("  x%d = %12.8f\n", i + 1, x[i]);
}

/* ------------------------------------------------------------------ */
/*  Infinity-norm of the difference between two vectors
 *  功能：计算两个向量差的无穷范数 (max|a_i - b_i|)
 *  用途：作为迭代收敛的判断依据。当误差小于 TOL 时，认为已收敛。
 * ------------------------------------------------------------------ */
static double inf_norm_diff(const double a[N], const double b[N])
{
    double max = 0.0;
    for (int i = 0; i < N; i++) {
        double d = fabs(a[i] - b[i]); // 计算当前分量的绝对误差
        if (d > max) max = d;         // 更新最大误差
    }
    return max;
}

/* ================================================================== */
/*  1. JACOBI METHOD (雅可比迭代法)
 *  核心思想：在第 k+1 次迭代计算 x_i 时，全部使用第 k 次迭代的旧值。
 *  公式：x_i^(k+1) = (b_i - Σ_{j!=i} A_ij * x_j^(k)) / A_ii
 *  特点：易于并行化，但收敛速度通常慢于 Gauss-Seidel。
 * ================================================================== */
void jacobi(void)
{
    printf("============================================================\n");
    printf(" Jacobi Iterative Method\n");
    printf("============================================================\n");

    /* 系数矩阵 A (Coefficient Matrix) 和 右端向量 b (RHS Vector) */
    double A[N][N] = {
            { 5, 2,  1},
            {-1,  4, 2},
            { 2, -3,  10}
    };
    double b[N] = {-12, 20, 3};

    double x[N]    = {0, 0, 0};   /* 当前迭代值 (current iterate), 初始化为 0 */
    double x_new[N];              /* 下一次迭代值 (next iterate), 暂存新结果 */
    int    iter    = 0;           /* 迭代计数器 */
    double err;                   /* 当前误差 */

    // 打印表头
    printf(" k  |   x1            x2            x3          | error\n");
    printf("----+------------------------------------------------+----------\n");

    do {
        /* Jacobi 更新核心循环：
         * 注意：计算 x_new[i] 时，右侧全部使用 x[j] (旧值)，
         * 即使 j < i 已经计算出新的 x_new[j]，在本轮也不使用。
         */
        for (int i = 0; i < N; i++) {
            double sigma = 0.0;
            // 计算非对角线元素的加权和：Σ_{j!=i} A_ij * x_j
            for (int j = 0; j < N; j++)
                if (j != i) sigma += A[i][j] * x[j];

            // 应用迭代公式
            x_new[i] = (b[i] - sigma) / A[i][i];
        }

        // 计算新旧向量之间的无穷范数误差
        err = inf_norm_diff(x_new, x);
        iter++;

        // 将新值覆盖旧值，准备下一轮迭代
        for (int i = 0; i < N; i++) x[i] = x_new[i];

        // 控制输出频率：前 10 次每次打印，之后每 10 次打印一次
        if (iter <= 10 || iter % 10 == 0)
            printf(" %-3d | %12.8f  %12.8f  %12.8f | %.2e\n",
                   iter, x[0], x[1], x[2], err);

    } while (err > TOL && iter < MAX_IT); // 收敛条件：误差<TOL 或 达到最大迭代次数

    printf("\n Converged after %d iterations\n", iter);
    printf(" Solution:\n");
    print_solution(x);
    printf("\n");
}

/* ================================================================== */
/*  2. GAUSS-SEIDEL METHOD (高斯 - 赛德尔迭代法)
 *  核心思想：在第 k+1 次迭代计算 x_i 时，对于 j < i 的分量，立即使用
 *  刚刚更新过的 x_j^(k+1) 新值；对于 j > i 的分量，仍使用旧值 x_j^(k)。
 *  公式：x_i^(k+1) = (b_i - Σ_{j<i} A_ij*x_j^(k+1) - Σ_{j>i} A_ij*x_j^(k)) / A_ii
 *  特点：通常比 Jacobi 收敛更快，因为利用了最新信息。
 * ================================================================== */
void gauss_seidel(void)
{
    printf("============================================================\n");
    printf(" Gauss-Seidel Iterative Method\n");
    printf("============================================================\n");

    double A[N][N] = {
            { 5, 2,  1},
            {-1,  4, 2},
            { 2, -3,  10}
    };
    double b[N] = {-12, 20, 3};

    double x[N]     = {0, 0, 0};
    double x_old[N]; // 用于保存上一轮迭代结果，以便计算误差
    int    iter     = 0;
    double err;

    printf(" k  |   x1            x2            x3          | error\n");
    printf("----+------------------------------------------------+----------\n");

    do {
        // 保存旧值用于误差计算 (因为 x 数组会在循环中被就地更新)
        for (int i = 0; i < N; i++) x_old[i] = x[i];

        /* Gauss-Seidel 更新核心循环：
         * 注意：计算 x[i] 时，如果 j < i，使用的 x[j] 已经是本轮更新过的新值。
         * 这体现了"即时更新"的特性。
         */
        for (int i = 0; i < N; i++) {
            double sigma = 0.0;
            for (int j = 0; j < N; j++)
                if (j != i) sigma += A[i][j] * x[j]; // x[j] 可能是新值也可能是旧值

            x[i] = (b[i] - sigma) / A[i][i]; // 就地更新
        }

        // 计算当前解与上一轮解的误差
        err = inf_norm_diff(x, x_old);
        iter++;

        if (iter <= 10 || iter % 10 == 0)
            printf(" %-3d | %12.8f  %12.8f  %12.8f | %.2e\n",
                   iter, x[0], x[1], x[2], err);

    } while (err > TOL && iter < MAX_IT);

    printf("\n Converged after %d iterations\n", iter);
    printf(" Solution:\n");
    print_solution(x);
    printf("\n");
}

/* ================================================================== */
/*  3. SOR METHOD (逐次超松弛迭代法)
 *  核心思想：在 Gauss-Seidel 的基础上引入松弛因子 omega (ω)。
 *  公式：x_i^(k+1) = (1-ω)*x_i^(k) + ω * (GS_step)
 *  其中 GS_step 是高斯 - 赛德尔的一步计算结果。
 *  - 若 ω = 1.0，退化为 Gauss-Seidel 方法。
 *  - 若 1.0 < ω < 2.0，称为超松弛，旨在加速收敛。
 *  - 若 0.0 < ω < 1.0，称为低松弛，有时用于改善稳定性。
 * ================================================================== */
void sor(double omega)
{
    printf("============================================================\n");
    printf(" SOR Method  (omega = %.2f)\n", omega);
    printf("============================================================\n");

    double A[N][N] = {
            { 5, 2,  1},
            {-1,  4, 2},
            { 2, -3,  10}
    };
    double b[N] = {-12, 20, 3};

    double x[N]     = {0, 0, 0};
    double x_old[N];
    int    iter     = 0;
    double err;

    printf(" k  |   x1            x2            x3          | error\n");
    printf("----+------------------------------------------------+----------\n");

    do {
        for (int i = 0; i < N; i++) x_old[i] = x[i];

        for (int i = 0; i < N; i++) {
            double sigma = 0.0;
            for (int j = 0; j < N; j++)
                if (j != i) sigma += A[i][j] * x[j];

            /* 第一步：计算标准的 Gauss-Seidel 预测值 */
            double gs = (b[i] - sigma) / A[i][i];

            /* 第二步：应用松弛因子进行加权平均 */
            // x_new = (1-ω)*x_old + ω*gs
            x[i] = (1.0 - omega) * x[i] + omega * gs;
        }

        err = inf_norm_diff(x, x_old);
        iter++;

        if (iter <= 10 || iter % 10 == 0)
            printf(" %-3d | %12.8f  %12.8f  %12.8f | %.2e\n",
                   iter, x[0], x[1], x[2], err);

    } while (err > TOL && iter < MAX_IT);

    printf("\n Converged after %d iterations\n", iter);
    printf(" Solution:\n");
    print_solution(x);
    printf("\n");
}

/* ================================================================== */
/*  MAIN FUNCTION
 *  程序入口：初始化参数并依次调用三种算法进行求解
 * ================================================================== */
int main(void)
{
    printf("\n");
    printf("  System of linear equations:\n");

    // 【已修复】此处打印的方程组系数现在与代码中定义的矩阵 A 完全一致
    printf("    5x1 +  2x2 + 1x3 =  -12\n");
    printf("   -x1 +  4x2 +  2x3  =  20\n");
    printf("   2x1 -  3x2 + 10x3 = 3\n");

    // 此处打印的精确解现在与代码中矩阵 A 和 b 的真实解一致
    printf("  Exact solution: x1=-4, x2=3, x3=2\n");
    printf("  Tolerance: %.0e,  Max iterations: %d\n\n", TOL, MAX_IT);

    // 执行 Jacobi 迭代
    jacobi();

    // 执行 Gauss-Seidel 迭代
    gauss_seidel();

    // 执行 SOR 迭代，测试不同的松弛因子
    // 循环条件改为动态计算数组长度，确保测试所有 4 个 omega 值
    double omegas[] = {0.9, 1.03, 1.0, 1.1};
    int num_omegas = sizeof(omegas) / sizeof(omegas[0]);

    for (int k = 0; k < num_omegas; k++)
        sor(omegas[k]);

    return 0;
}