迭代法求解线性方程组(C++)

尽管所有的线性方程组都可以通过Gauss消去法和LU分解法求解, 但是对于大型的稀疏线性代数方程组, 往往求解效率较低. 这时候常使用的求解算法是迭代法, 其基本形式为:

$${\mathbf{x}}^{(k+1)}={\mathbf{F}}_k({\mathbf{x}}^{(k)},{\mathbf{x}}^{(k?1)},?,{\mathbf{x}}^{(0)}),k=0,1,?$$

若${\mathbf{x}}^{(k+1)}$只与${\mathbf{x}}^{(k)}$有关, 且$F_k$是线性的, 即

$${\mathbf{x}}^{(k+1)}={\mathbf{B}}_k{\mathbf{x}}^{(k)}+{\mathbf{f}}_k$$

其中, ${\mathbf{B}}_k\in {\mathbb{R}}^{n\times n}$, 称为单步线性迭代法, ${\mathbf{B}}_k$称为迭代矩阵. 若${\mathbf{B}}_k$与${\mathbf{f}}_k$都与$k$无关, 即

$${\mathbf{x}}^{(k+1)}=\mathbf{B}{\mathbf{x}}^{(k)}+\mathbf{f}$$

称为单步定常线性迭代法. 本文主要讨论这种类型的迭代法.

求解算法

Jacobi迭代法

设有线性方程组

$$\mathbf{A}\mathbf{x}=\mathbf{b}$$

其中, $\mathbf{A}=(a_{ij})$非奇异, 我们将$\mathbf{A}$的对角线$\mathbf{D}$, 上三角部分$\mathbf{U}$, 下三角部分$\mathbf{L}$分离, 即

$$\mathbf{A}=\mathbf{L}+\mathbf{D}+\mathbf{U}$$

若$\mathbf{D}$非奇异, 我们可以将$\mathbf{L}\mathbf{x}$和$\mathbf{U}\mathbf{x}$移到方程的右端, 再将方程两端同时左乘${\mathbf{D}}^{?1}$, 得

$$\mathbf{x}=?{\mathbf{D}}^{?1}(\mathbf{L}+\mathbf{U})\mathbf{x}+{\mathbf{D}}^{?1}\mathbf{b}$$

若我们取${\mathbf{B}}_J:=?{\mathbf{D}}^{?1}(\mathbf{L}+\mathbf{U}),{\mathbf{f}}_J:={\mathbf{D}}^{?1}\mathbf{b}$, 则有如下迭代公式

$${\mathbf{x}}^{(k+1)}={\mathbf{B}}_J\mathbf{x}+{\mathbf{f}}_J$$

上式被称为Jacobi迭代法. 因此, 我们可以将算法表示为如下伪代码:

Gauss-Seidel迭代法

若我们把Jacobi迭代法还原为方程组的形式, 我们有

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i?\sum _{j=1}^{i?1}{a}_{ij}{x}_j^{(k)}?\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}),i=1,2,?,n$$

观察上式可知, 在实际计算$x_i^{(k+1)}$的时候, 分量$x_1^{(k+1)},?,x_{i?1}^{(k+1)}$已经算出, 可以将$x_1^{(k)},?,x_{i?1}^{(k)}$替换为$x_1^{(k+1)},?,x_{i?1}^{(k+1)}$, 这样每次迭代都可以用到当前最新的计算结果, 这就是Gauss-Seidel迭代法. 相应地, Gauss-Seidel迭代法也有矩阵形式:

$${\mathbf{x}}^{(k+1)}={\mathbf{B}}_{GS}\mathbf{x}+{\mathbf{f}}_{GS}$$

其中, ${\mathbf{B}}_{GS}:=?(\mathbf{D}+\mathbf{L}{)}^{?1}\mathbf{U},{\mathbf{f}}_{GS}:=(\mathbf{D}+\mathbf{L}{)}^{?1}\mathbf{b}$. 同样地, Gauss-Seidel迭代法可表示为伪代码如下:

超松弛迭代法

将Gauss-Seidel迭代法和Jacobi迭代法作加权平均, 我们就得到了超松弛迭代法(又叫SOR法). 其中Gauss-Seidel迭代法的权重$\omega$称为松弛因子.

同以上两种迭代法, 超松弛迭代法也可以表示为如下的矩阵形式

$${\mathbf{x}}^{(k+1)}={\mathbf{B}}_{\omega }\mathbf{x}+{\mathbf{f}}_{\omega }$$

其中, ${\mathbf{B}}_{\omega }:=(\mathbf{D}+\omega \mathbf{L}{)}^{?1}((1?\omega )\mathbf{D}?\omega \mathbf{U}),{\mathbf{f}}_{\omega }:=\omega (\mathbf{D}+\omega \mathbf{L}{)}^{?1}\mathbf{b}$. 同样地, 超松弛迭代法可表示为伪代码如下:

算法实现

准备工作

首先进行预处理如下:

#include <armadillo>
#include <vector>
#include <QFile>
#include <QString>
#include <string>
using namespace arma;
using namespace std;

为了方便写入向量, 编写一个向量导出函数.¹

void write_vec(QFile &file, const vec &x)
{
    string s;
    for (auto &i : x)
        s += (to_string(i)) += ',';
    *(s.end() - 1) = '\n';
    file.write(s.c_str());
}

Jacobi迭代法

矩阵法

/*
 * Jacobi迭代法(矩阵法)
 * A:系数矩阵
 * b:常数向量
 * x:迭代初始向量
 * u:迭代过程
 * e:迭代精度
 */
void Jacobi_mat(const mat &A, const vec &b, vec &x, vector<vec> &u, const double &e = 1e-6)
{
    vec x1(x), b1(b.n_elem, 1, fill::none);
    u.push_back(x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        A1.at(i, i) = 0;
        x.at(i) = b1.at(i) = b.at(i) / A.at(i, i);
        for (unsigned j(0); j != A.n_cols; ++j)
            if (i != j)
                x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x1.at(j);
    }
    u.push_back(x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        u.push_back(x = b1 - A1 * x);
    }
}

分量法

/*
 * Jacobi迭代法(分量法)
 * A:系数矩阵
 * b:常数向量
 * x:迭代初始向量
 * u:迭代过程
 * e:迭代精度
 */
void Jacobi_ele(const mat &A, const vec &b, vec &x, vector<vec> &u, const double &e = 1e-6)
{
    vec x1(x), b1(b.n_elem, 1, fill::none);
    u.push_back(x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        x.at(i) = b1.at(i) = b.at(i) / A.at(i, i);
        for (unsigned j(0); j != A.n_cols; ++j)
            if (i != j)
                x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x1.at(j);
    }
    u.push_back(x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        for (unsigned i(0); i != x.n_rows; ++i)
        {
            x.at(i) = b1.at(i);
            for (unsigned j(0); j != A.n_cols; ++j)
                if (i != j)
                    x.at(i) -= A1.at(i, j) * x1.at(j);
        }
        u.push_back(x);
    }
}

/*
 * Jacobi迭代法(分量法)
 * A:系数矩阵
 * b:常数向量
 * x:迭代初始向量
 * p:迭代过程保存路径
 * e:迭代精度
 */
void Jacobi(const mat &A, const vec &b, vec &x, const QString &p, const double &e = 1e-6)
{
    QFile file(p);
    if (!file.open(QIODevice::WriteOnly))
        throw "文件保存失败！";
    vec x1(x), b1(b.n_elem, 1, fill::none);
    write_vec(file, x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        x.at(i) = b1.at(i) = b.at(i) / A.at(i, i);
        for (unsigned j(0); j != A.n_cols; ++j)
            if (i != j)
                x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x1.at(j);
    }
    write_vec(file, x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        for (unsigned i(0); i != x.n_rows; ++i)
        {
            x.at(i) = b1.at(i);
            for (unsigned j(0); j != A.n_cols; ++j)
                if (i != j)
                    x.at(i) -= A1.at(i, j) * x1.at(j);
        }
        write_vec(file, x);
    }
    file.close();
}

Gauss-Seidel迭代法

矩阵法

void Gauss_Seidel_mat(const mat &A, const vec &b, vec &x, vector<vec> &u, const double &e = 1e-6)
{
    mat T(inv(trimatl(A))), A1(T * trimatu(A, 1));
    vec b1(T * b), x1;
    u.clear();
    u.push_back(x);
    do
    {
        x1 = x;
        u.push_back(x = b1 - A1 * x);
    } while (norm(x - x1) > e);
}

分量法

/*
 * Gauss-Seidel迭代法(分量法)
 * A:系数矩阵
 * b:常数向量
 * x:迭代初始向量
 * u:迭代过程
 * e:迭代精度
 */
void Gauss_Seidel_ele(const mat &A, const vec &b, vec &x, vector<vec> &u, const double &e = 1e-6)
{
    vec x1(x), b1(b.n_elem, 1, fill::none);
    u.push_back(x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        x.at(i) = b1.at(i) = b.at(i) / A.at(i, i);
        unsigned j(-1);
        while (++j != i)
            x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x.at(j);
        while (++j != A.n_cols)
            x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x.at(j);
    }
    u.push_back(x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        for (unsigned i(0); i != x.n_rows; ++i)
        {
            x.at(i) = b1.at(i);
            unsigned j(-1);
            while (++j != i)
                x.at(i) -= A1.at(i, j) * x.at(j);
            while (++j != A.n_cols)
                x.at(i) -= A1.at(i, j) * x.at(j);
        }
        u.push_back(x);
    }
}

/*
 * Gauss-Seidel迭代法(分量法)
 * A:系数矩阵
 * b:常数向量
 * x:迭代初始向量
 * p:迭代过程保存路径
 * e:迭代精度
 */
void Gauss_Seidel(const mat &A, const vec &b, vec &x, const QString &p, const double &e = 1e-6)
{
    QFile file(p);
    if (!file.open(QIODevice::WriteOnly))
        throw "文件保存失败！";
    vec x1(x), b1(b.n_elem, 1, fill::none);
    write_vec(file, x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        x.at(i) = b1.at(i) = b.at(i) / A.at(i, i);
        unsigned j(-1);
        while (++j != i)
            x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x.at(j);
        while (++j != A.n_cols)
            x.at(i) -= (A1.at(i, j) = A.at(i, j) / A.at(i, i)) * x.at(j);
    }
    write_vec(file, x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        for (unsigned i(0); i != x.n_rows; ++i)
        {
            x.at(i) = b1.at(i);
            unsigned j(-1);
            while (++j != i)
                x.at(i) -= A1.at(i, j) * x.at(j);
            while (++j != A.n_cols)
                x.at(i) -= A1.at(i, j) * x.at(j);
        }
        write_vec(file, x);
    }
    file.close();
}

超松弛迭代法

矩阵法

/*
 * SOR法(矩阵法)
 * A:系数矩阵
 * b:常数向量
 * w:松弛因子
 * x:迭代初始向量
 * u:迭代过程
 * e:迭代精度
 */
void SOR_mat(const mat &A, const vec &b, const double &w, vec &x, vector<vec> &u, const double &e = 1e-6)
{
    mat L(trimatl(A, -1)), D(diagmat(A)), U(trimatu(A, 1)), T(inv(trimatl(w * L + D))), A1(T * ((1 - w) * D - w * U));
    vec b1(T * b * w), x1;
    u.clear();
    u.push_back(x);
    do
    {
        x1 = x;
        u.push_back(x = b1 + A1 * x);
    } while (norm(x - x1) > e);
}

分量法

/*
 * SOR法(分量法)
 * A:系数矩阵
 * b:常数向量
 * w:松弛因子
 * x:迭代初始向量
 * u:迭代过程
 * e:迭代精度
 */
void SOR_ele(const mat &A, const vec &b, const double &w, vec &x, vector<vec> &u, const double &e = 1e-6)
{
    double w1(1 - w);
    vec x1(x), b1(b.n_elem, 1, fill::none);
    u.push_back(x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        (x.at(i) *= w1) += b1.at(i) = w * b.at(i) / A.at(i, i);
        unsigned j(-1);
        while (++j != i)
            x.at(i) -= (A1.at(i, j) = w * A.at(i, j) / A.at(i, i)) * x.at(j);
        while (++j != A.n_cols)
            x.at(i) -= (A1.at(i, j) = w * A.at(i, j) / A.at(i, i)) * x.at(j);
    }
    u.push_back(x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        for (unsigned i(0); i != x.n_rows; ++i)
        {
            (x.at(i) *= w1) += b1.at(i);
            unsigned j(-1);
            while (++j != i)
                x.at(i) -= A1.at(i, j) * x.at(j);
            while (++j != A.n_cols)
                x.at(i) -= A1.at(i, j) * x.at(j);
        }
        u.push_back(x);
    }
}

/*
 * SOR法(分量法)
 * A:系数矩阵
 * b:常数向量
 * w:松弛因子
 * x:迭代初始向量
 * p:迭代过程保存路径
 * e:迭代精度
 */
void SOR(const mat &A, const vec &b, const double &w, vec &x, const QString &p, const double &e = 1e-6)
{
    QFile file(p);
    if (!file.open(QIODevice::WriteOnly))
        throw "文件保存失败！";
    double w1(1 - w);
    vec x1(x), b1(b.n_elem, 1, fill::none);
    write_vec(file, x);
    mat A1(A.n_rows, A.n_cols, fill::none);
    for (unsigned i(0); i != x.n_rows; ++i)
    {
        if (abs(A.at(i, i)) < e)
            throw "对角线元素为零！";
        (x.at(i) *= w1) += b1.at(i) = w * b.at(i) / A.at(i, i);
        unsigned j(-1);
        while (++j != i)
            x.at(i) -= (A1.at(i, j) = w * A.at(i, j) / A.at(i, i)) * x.at(j);
        while (++j != A.n_cols)
            x.at(i) -= (A1.at(i, j) = w * A.at(i, j) / A.at(i, i)) * x.at(j);
    }
    write_vec(file, x);
    while (norm(x1 - x) > e)
    {
        x1 = x;
        for (unsigned i(0); i != x.n_rows; ++i)
        {
            (x.at(i) *= w1) += b1.at(i);
            unsigned j(-1);
            while (++j != i)
                x.at(i) -= A1.at(i, j) * x.at(j);
            while (++j != A.n_cols)
                x.at(i) -= A1.at(i, j) * x.at(j);
        }
        write_vec(file, x);
    }
    file.close();
}

实例分析

例1

用迭代法求解线性方程组

$$\left(\begin{array}{ccc}8& 1& 2\\ 8& 7& 2\\ 4& 9& 9\end{array}\right)\mathbf{x}=\left(\begin{array}{c}10\\ 18\\ 17\end{array}\right)$$

选取初值${\mathbf{x}}_0=(0,0,0{)}^T$, 代入程序并测量时间得解向量为$\mathbf{x}=(1.0625,1.3333,0.0833{)}^T$, 运行时间见下表:

算法实现方式	耗时($\mu$s)
Jacobi迭代法(分量法)	0.0393829
Jacobi迭代法(矩阵法)	0.0456848
Gauss-Seidel迭代法(分量法)	0.00231934
Gauss-Seidel迭代法(矩阵法)	0.00382996
SOR法(分量法)	0.00256348
SOR法(矩阵法)	0.00396729

其中, 超松弛迭代法的松弛因子取1.2. 由表可知, 矩阵法实现比分量法实现较慢, 但可读性更高, 这是因为矩阵法实现涉及到矩阵求逆运算, 耗时较多; 在迭代方法方面, Gauss-Seidel迭代法和超松弛迭代法较快, 而Jacobi迭代较慢.

例2

用迭代法求解线性方程组

$$?\begin{array}{l}7x_1+x_2+2x_3=10\\ x_1+8x_2+2x_3=8\\ 2x_1+2x_2+9x_3=6\end{array}$$

取初值为$(0,0,0{)}^T$, 利用Jacobi迭代法有迭代过程如下表:

$k$	$x_1^{(k)}$	$x_2^{(k)}$	$x_3^{(k)}$
0	0.000000	0.000000	0.000000
1	1.428571	1.000000	0.666667
2	1.095238	0.654762	0.126984
3	1.298753	0.831349	0.277778
4	1.230442	0.768211	0.193311
5	1.263595	0.797867	0.222521
6	1.251013	0.786420	0.208564
7	1.256636	0.791482	0.213904
8	1.254387	0.789445	0.211529
9	1.255357	0.790319	0.212482
10	1.254960	0.789960	0.212072
11	1.255128	0.790112	0.212240
12	1.255058	0.790049	0.212169
13	1.255088	0.790076	0.212198
14	1.255075	0.790064	0.212186
15	1.255081	0.790069	0.212191
16	1.255078	0.790067	0.212189
17	1.255079	0.790068	0.212190
18	1.255079	0.790068	0.212190

而利用Gauss-Seidel迭代法有迭代过程如下表:

$k$	$x_1^{(k)}$	$x_2^{(k)}$	$x_3^{(k)}$
0	0.000000	0.000000	0.000000
1	1.428571	0.821429	0.166667
2	1.263605	0.800383	0.208003
3	1.254802	0.791149	0.212011
4	1.254976	0.790125	0.212200
5	1.255068	0.790067	0.212192
6	1.255078	0.790067	0.212190
7	1.255079	0.790068	0.212190

由此可见, Gauss-Seidel迭代法在本例中收敛速度较快.

---

1. 这里用了Qt的QFile, 也可以方便地改写为使用标准库的FILE或文件输出流. ??