基于C语言的Jacobi迭代和Gauss-Seidel迭代的方程组求解实现

Jacobi迭代方法介绍

Jacobi迭代法是一种简单的迭代求解方法，适用于严格对角占优矩阵。其基本思想是利用当前迭代步的已知解来更新下一个迭代步的解。在C语言实现中，我们首先需要定义系数矩阵A、常数向量b以及迭代向量x。然后，通过循环迭代，不断更新x的值，直到满足收敛条件或达到最大迭代次数。

Jacobi迭代法的优点是简单直观，但收敛速度相对较慢。特别是对于非严格对角占优的矩阵，Jacobi迭代法可能不收敛。因此，在实际应用中，我们需要根据问题的具体情况选择合适的迭代方法。

Gauss-Seidel迭代方法介绍

Gauss-Seidel迭代法是对Jacobi迭代法的一种改进。在Gauss-Seidel迭代法中，我们使用当前迭代步已更新的解来更新下一个解。这种策略可以加速收敛速度，特别是在对角线元素较大的情况下。在C语言实现中，我们需要对系数矩阵A进行适当的分解，以便在迭代过程中方便地获取已更新的解。

与Jacobi迭代法相比，Gauss-Seidel迭代法具有更快的收敛速度。但是，Gauss-Seidel迭代法并不总是收敛的，其收敛性取决于系数矩阵A的性质。因此，在使用Gauss-Seidel迭代法时，我们需要确保系数矩阵A满足一定的条件（如严格对角占优）。

具体代码实现

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "time.h"

#define N 3   //要求解的数组大小N x N
#define ERROR 1e-6  //误差要求
#define IterMaxNum 50000   //最大迭代次数

float current_error(float* last, float* now)
{
    float max_error = 0;
    float temp;
    for (int i = 0; i < N; i++)
    {
        temp = *(now + i) - *(last + i);
        if (temp < 0)
            temp = -temp;
        if (temp > max_error)
            max_error = temp;
    }
    return max_error;
}

//雅可比迭代函数
float* Jacobi(float* A, float* b)
{
    clock_t start, end;
    double cpu_time_used;

    // 记录开始时间
    start = clock();
    int i, j;
    int iter_num = 0;   //目前的迭代次数
    float error_now = 100;   //目前迭代下的误差

    //生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
    float* x = (float*)malloc(N * sizeof(float));
    memset(x, 0, N * sizeof(float));
    float* x0 = (float*)malloc(N * sizeof(float));
    memset(x0, 0, N * sizeof(float));

    while(error_now > ERROR)
    {
        //将上一次迭代的结果保存到x0中去
        memcpy(x0, x, N * sizeof(float));

        for (i = 0; i < N; i++)
        {
            float sum = 0;
            for (j = 0; j < N; j++)
            {
                //printf("%f\t", A[i * N + j]);
                if (i != j)
                    sum += A[i * N + j] * (x0[j]);
            }
            x[i] = (b[i] - sum) / A[i * N + i];
        }

        error_now = current_error(x, x0);
        iter_num++;

        printf("当前迭代次数为:【%d】，迭代误差为:【%f】\n", iter_num, error_now);
        printf("迭代结果为：%f,%f,%f\n", x[0], x[1], x[2]);
        if (iter_num > IterMaxNum)
        {
            printf("迭代次数超过最大值要求，返回最后一次迭代结果");
            break;
        }
    }
    // 记录结束时间
    end = clock();

    // 计算经过的CPU时间
    cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

    // 打印Jacobi迭代的运行时间
    printf("\n【Jacobi迭代运行了 %f 秒。】\n", cpu_time_used);

    free(x0);
    return x;
}


float* GS(float* A, float* b)
{
    int i, j;
    int iter_num = 0;   //目前的迭代次数
    float error_now = 100;   //目前迭代下的误差

    //生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
    float* x = (float*)malloc(N * sizeof(float));
    memset(x, 0, N * sizeof(float));
    float* x0 = (float*)malloc(N * sizeof(float));
    memset(x0, 0, N * sizeof(float));

    clock_t start, end;
    double cpu_time_used;

    // 记录开始时间
    start = clock();
    while (error_now > ERROR)
    {
        memcpy(x0, x, N * sizeof(float));

        for (i = 0; i < N; i++)
        {
            float sum = 0;
            for (j = 0; j < N; j++)
            {
                //printf("%f\t", A[i * N + j]);
                if (i != j)
                    sum += A[i * N + j] * (x[j]);
            }
            x[i] = (b[i] - sum) / A[i * N + i];
        }

        error_now = current_error(x, x0);  //获得当前迭代的误差结果
        iter_num++;   //迭代次数加1

        printf("当前迭代次数为:【%d】，迭代误差为:【%f】\n", iter_num, error_now);
        printf("迭代结果为：%f,%f,%f\n", x[0], x[1], x[2]);
        if (iter_num > IterMaxNum)
        {
            printf("迭代次数超过最大值要求，返回最后一次迭代结果");
            break;
        }
    }
    // 记录结束时间
    end = clock();

    // 计算经过的CPU时间
    cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

    // 打印Gauss-Seidel迭代的运行时间
    printf("\n【Gauss-Seidel迭代运行了 %f 秒。】\n", cpu_time_used);
    free(x0);
    return x;
}


int main()
{
    // 例子
    float A[N][N] = {
            {10,-1,-2},
            {-1,10,-2},
            {-1, -1, 5},
    };   //要求解的矩阵
    float b[N] = {72, 83, 42};

    float* x;

//    x = Jacobi((float*)A, b);

    x = GS((float*)A, b);

    printf("\n【最终结果为】：%f, %f, %f\n", x[0], x[1], x[2]);

    free(x);   //释放x
}

示例题目

以如下方程组为例，进行迭代方法的测试

实现效果

Jacobi迭代实现效果，迭代了【17】达到了收敛

Gauss-Seidel迭代实现效果，而Gauss-Seidel只迭代了【10】次即可达到收敛。

注意，有些示例中，Jacobi迭代是比Gauss-Seidel收敛的快的，但更多的时候Gauss-Seidel是收敛更快的，并且存在有时候Jacobi收敛，而Gauss-Seidel却发散，要根据情况而定。