2-2、线性代数方程组的求解

内容来自王晓华老师

这块内容有点硬核，先做了解，主要学习如何使用迭代解决问题的步骤

两种常用的迭代法求解方程组，分别是雅可比迭代法（Jacobi）和高斯—赛德尔迭代法（Gauss—Seidel）

完整的雅可比迭代算法

const double PRECISION = 0.000000001;
const int VEC_N = 16;  //实际方程组的个数 n 必须小于VEC_N

void jacobi_eqn(double a[][VEC_N], double b[], int n, double x[])
{
    double xt[VEC_N];
    double max_diff = 0.0;

    do
    {
        for (int i = 0; i < n; i++)
        {
            double sum = 0.0;
            for (int k = 0; k < n; k++)
            {
                if (i != k)  //对角线元素不计算
                {
                    sum += a[i][k] * x[k];
                }
            }
            xt[i] = (b[i] - sum) / a[i][i];   //根据关系计算 xi 的新值
        }

        max_diff = calc_max_diff(xt, x, n);

        for (int j = 0; j < n; j++)
        {
            x[j] = xt[j]; //更新 x，准备下一次迭代
        }
    } while (max_diff > PRECISION);
}

高斯-赛德尔迭代法
雅可比迭代法每次迭代计算时，将上一次的迭代变量整体带入到迭代关系式中，计算新的迭代变量值，也就是所谓的整体迭代。在迭代收敛的前提下，如果迭代变量中的每个分量 x 在计算出新的迭代值后，直接带入迭代，参与其他迭代分量的计算，则能显著地提高迭代效果，这种改进的方法就是高斯-赛德尔迭代法。
从算法实现的角度理解，这种改进思想就是每计算出一个迭代分量的新迭代值，就立即让它参与到其他迭代分量的计算中，其迭代关系可以理解为：

const double PRECISION = 0.000000001;
const int VEC_N = 16;  //实际方程组的个数 n 必须小于 VEC_N

void gauss_seidel_eqn(double a[][VEC_N], double b[], int n, double x[])
{
    double max_diff = 0.0;

    do
    {
        max_diff = 0.0;
        for (int i = 0; i < n; i++)
        {
            double sum = 0.0;
            for (int k = 0; k < n; k++)
            {
                if (i != k)  //对角线元素不计算
                {
                    sum += a[i][k] * x[k];
                }
            }
            double xt = (b[i] - sum) / a[i][i];   //根据关系计算 xi 的新值
            if (std::fabs(xt - x[i]) > max_diff) //max_diff 只保留差值最大的
            {
                max_diff = std::fabs(xt - x[i]);
            }
            x[i] = xt;
        }
    } while (max_diff > PRECISION);
}