[计算机数值分析]高斯-塞德尔迭代公式解线性方程组

Spring-_-Bear 的 CSDN 博客导航

---

在 雅可比迭代公式 的基础上，设已得到近似值（$x_1^{(k)},x_2^{(k)},x_3^{(k)}$），则由方程 (1)₁

$$\{\begin{array}{l}{x}_1=0.1x_2+0.2x_3+0.72\\ x_2=0.1x_1+0.2x_3+0.83\\ x_3=0.2x_1+0.2x_2+0.84\end{array}\text{\hspace{1em}(1)}$$

得

$$x_1^{(k+1)}=0.1x_2^{(k)}+0.2x_3^{(k)}+0.72$$

对于收敛的迭代过程，所求出的 “新值” $x_1^{(k+1)}$ 常比 “老值” $x_1^{(k)}$ 更准确些，因此可用它替代老值 $x_1^{(k)}$ 作进一步的计算，将 $x_1^{(k+1)}$ 与 $x_3^{(k)}$ 一起代入 (1)₂ 得

$$x_2^{(k+1)}=0.1x_1^{(k+1)}+0.2x_3^{(k)}+0.83$$

处于同样得考虑，我们再用新值 $x_2^{(k+1)}$ 替代老值 $x_2^{(1)}$，由方程 (1)₃ 得

$$x_3^{(k+1)}=0.2x_1^{(k+1)}+0.2x_2^{(k+1)}+0.84$$

这种充分利用新值建立起来的迭代公式称作解方程组（1）的高斯 - 塞德尔公式，如（2）所示。

$$\{\begin{array}{l}{x}_1^{(k+1)}=0.1x_2^{(k)}+0.2x_3^{(k)}+0.72\\ x_2^{(k+1)}=0.1x_1^{(k+1)}+0.2x_3^{(k)}+0.83\\ x_3^{(k+1)}=0.2x_1^{(k+1)}+0.2x_2^{(k+1)}+0.84\end{array}\text{\hspace{1em}(2)}$$

仍取初值 $x_1^{(0)}=x_2^{(0)}=x_3^{(0)}=0$，按式（2）的计算结果如下表，与 雅可比迭代公式 相比效果更好。

k	x1(k)	x2(k)	x3(k)
0	0.00000	0.00000	0.00000
1	0.72000	0.90200	1.16440
2	1.04308	1.16719	1.28205
3	1.09313	1.19572	1.29777
4	1.09913	1.19947	1.29972
5	1.09989	1.19993	1.29997
6	1.09999	1.19999	1.30000

运行示例：

程序源码：

#include <iostream>
#include <iomanip>

#define MAX 10

using namespace std;

int main(void)
{
    int row;
    cout << "请输入增广矩阵行数:";
    cin >> row;

    int col;
    cout << "请输入增广矩阵列数:";
    cin >> col;

    double a[MAX][MAX];
    for (int i = 1; i <= row; i++)
    {
        cout << "请输入增广矩阵第 " << i << " 行的元素：";
        for (int j = 1; j <= col; j++)
        {
            cin >> a[i][j];
        }
    }

    // 对增广矩阵元素进行处理化简
    for (int i = 1; i <= row; i++)
    {
        for (int j = 1; j <= col; j++)
        {
            if (i != j && j != col)
            {
                // 非行列坐标相等的元素除以负行列坐标相等的元素
                a[i][j] /= -a[i][i];
            }
            if (j == col)
            {
                // 常数项除以行列相等的元素
                a[i][j] /= a[i][i];
            }
        }
    }

    // 将矩阵主对角线上的元素赋值为 0
    for (int i = 1; i <= row; i++)
    {
        a[i][i] = 0;
    }

    double b[MAX];
    cout << "请输入迭代初值:";
    for (int i = 1; i < col; i++)
    {
        cin >> b[i];
    }

    int N;
    cout << "请输入最大迭代次数:";
    cin >> N;

    double accuracy;
    cout << "请输入精度：";
    cin >> accuracy;

    double error = 0;
    for (int i = 1; i <= N; i++)
    {
        cout << "第 " << i << " 次迭代:";

        // 用 count 控制迭代值按每行 col-1 个输出，即自变量个数输出
        int count = 1;
        // 每一次迭代 (i)，求取每一行的行列坐标相等的值即 a[j][j]
        for (int j = 1; j <= row; j++)
        {
            // j 循环 col 次，求得 a[j][j]
            double item;
            for (int k = 1; k <= col; k++)
            {
                if (j == k)
                {
                    continue;
                }

                if (k != col)
                {
                    // 如果 j 不等于 k 且不等于常数项的下标，则每一项等于系数与对应迭代初值的乘积
                    item = a[j][k] * b[k];
                }
                else
                {
                    // j=col 即常数项下标时，直接加上常数项，不对常数项处理
                    item = a[j][col];
                }

                // 累加 item，求得 a[1][1], a[2][2]......
                a[j][j] += item;
            }
            // 输出 a[1][1], a[2][2]......
            cout << "\t\t" << a[j][j] << "\t\t";

            // 相邻两次迭代结果的偏差
            error = abs(a[count][count] - b[count]);
            // 更新 b[] 里的值为本次迭代结束的值
            b[count] = a[count][count];
            // 赋 a[1][1], a[2][2]......值为 0，避免拥有初值影响下一次迭代结果
            a[count][count] = 0;

            if (count % (col - 1) == 0)
                cout << endl;
            count++;
        }

        // 解符合精度要求，直接结束迭代过程
        if (error <= accuracy)
        {
            break;
        }

        if (i > N)
        {
            cout << "达到允许的最大迭代次数，未找到符合精度要求的根！" << endl;
            break;
        }
    }

    return 0;
}