迭代法求解线性方程组（C++实现）

本系列是数值分析相关算法的文章，这次用迭代法求解线性方程组，不同于上次用高斯消元法之类的求解。迭代法对于稀疏矩阵方程组的运算，会大大提高。而如果用高斯相关的算法求解，会浪费大量资源计算无用的东西，所以有必要研究此算法。
本文章主要使用了3个算法，分别是雅可比迭代法、高斯-塞德尔迭代法、超松弛迭代法。每一个算法都比前一次的迭代次数少。下面贴代码！

#include <iostream>
#include<cmath>
using namespace std;
//数组的大小
const int n = 4;

//子程序，用来向量的无穷范数，该定义可以参看《数值分析》（第5版——清华大学出版社）P164页
template<class T>
T MAX(T a[n])
{
    T max = 0;
    for (int i = 0; i<n; i++)
    {
        if (fabs(a[i])>max)
            max = fabs(a[i]);
    }
    return max;
}
//求解线性方程组的雅可比迭代法
template<class T>
void Jacobi(T a[n][n], T b[n])
{
    int count = 0;
    T aaa;
    T bbb;
    int i, j;
    T  ccc;
    T xk[n] = { 0 };
    double x0[n];
    double xxx;
    bool ddd;
    //求解的过程，递归求解。
    do {
        count++;
        for (j = 0; j<n; j++)
            x0[j] = xk[j];

        for (i = 0; i<n; i++)
        {
            T sum = 0;
            for (j = 0; j<n; j++)
            {
                if (j == i)
                    continue;
                else
                    sum += a[i][j] * x0[j];
            }
            xk[i] = (b[i] - sum) / a[i][i];
    //      cout << "xk[" << i + 1 << "]=" << xk[i] << endl;
        }
        aaa = MAX(xk);                           //求出xk的无穷范数
    //  cout << aaa << endl;  
        bbb = MAX(x0);                           //求出x0的无穷范数
    //  cout << bbb << endl;
        ccc = fabs(bbb - aaa);                     //做差
    //  cout << ccc << endl;
        xxx = fabs(ccc - 0.000001);                       //与要满足的误差精度作比较，我设定的是迭代精度为0.000001
    //  cout << xxx << endl;
        if (xxx<0.000001)
            ddd = 0;
        else
            ddd = 1;
    } while (ddd);                           //判断是否在迭代精度外

    cout << "this is JACOBI method  " << endl;                    //输出最后结果
    for (i = 0; i<n; i++)
        cout << "x[" << i + 1 << "]= " << xk[i] << endl;
    cout << "recursive times:  " << count << endl;
}

//高斯-塞德尔迭代法
template<class T>
void Gausssin_Jacobi(T a[n][n], T b[n])
{
    int count = 0;
    T aaa;
    T bbb;
    int i, j;
    T  ccc;
    double xxx;
    bool ddd;
    T x0[n];
    T xk[n] = { 0 };
    //迭代计算过程
    do {
        count++;          //记录迭代次数
        //先复制数组，可以与后面计算后的迭代结果做差
        for (j = 0; j<n; j++)
            x0[j] = xk[j];
        //计算过程
        for (i = 0; i<n; i++)
        {
            T sum1 = 0;
            T sum2 = 0;
            for (j = 0; j<=i - 1; j++)
            {
                sum1 += a[i][j] * xk[j];
            }
            for (j = i + 1; j<n; j++)
                sum2 += a[i][j] * x0[j];

            xk[i] = (b[i] - sum1 - sum2) / a[i][i];
        //  cout << "this is the XK[" << i + 1 << "]=" << xk[i] << endl;
        }
        aaa = MAX(xk);                                         //求出xk的范数
    //  cout << "this is XK de fanshu" << aaa << endl;
        bbb = MAX(x0);                                         //求出x0的范数
    //  cout << "this is x0 de fanshu" << bbb << endl;
        ccc = fabs(bbb - aaa);                                 //范数之差
    //  cout << "cha zhi de daxiao" << ccc << endl;
        xxx = fabs(ccc - 0.000001);                           //范数之差与精度的比较
    //  cout << "this is panduan de daxiao" << xxx << endl << endl;
        if (xxx<0.000001)
            ddd = 0;
        else
            ddd = 1;
    } while (ddd);                                 

    //输出结果
    cout << "this is Gausssin_Jacobi method  " << endl;
    for (i = 0; i<n; i++)
        cout << "x[" << i + 1 << "]= " << xk[i] << endl;
     cout << "recursive times:" << count << endl;
}


//超松弛迭代法（successive over relaxation method）（SOR方法）
template<class T>
void SOR(T a[n][n], T b[n],double w)
{
    int count = 0;
    T aaa;
    T bbb;
    int i, j;
    T  ccc;
    double xxx;
    bool ddd;
    T x0[n];
    T xk[n] = { 0 };
    T xk_[n] = { 0 };
    //迭代计算过程
    do {
        count++;               //记录迭代次数
        //先复制数组，可以与后面计算后的迭代结果做差
        for (j = 0; j<n; j++)
            x0[j] = xk[j];
        //计算过程
        for (i = 0; i<n; i++)
        {
            T sum1 = 0;
            T sum2 = 0;
            for (j = 0; j <= i - 1; j++)
            {
                sum1 += a[i][j] * xk[j];
            }
            for (j = i+1 ; j<n; j++)
                sum2 += a[i][j] * x0[j];

            xk_[i] = (b[i] - sum1 - sum2) / a[i][i];
            xk[i] = (1 - w)*x0[i] + w*xk_[i];
  //    cout << "this is the XK[" << i + 1 << "]=" << xk[i] << endl;
        }
        aaa = MAX(xk);                                         //求出xk的范数
//      cout << "this is XK de fanshu" << aaa << endl;
        bbb = MAX(x0);                                         //求出x0的范数
   //   cout << "this is x0 de fanshu" << bbb << endl;
        ccc = fabs(bbb - aaa);                                 //范数之差
   //   cout << "cha zhi de daxiao" << ccc << endl;
        xxx = fabs(ccc - 0.000001);                           //范数之差与精度的比较
  //    cout << "this is panduan de daxiao" << xxx << endl << endl;
        if (xxx<0.000001)
            ddd = 0;
        else
            ddd = 1;
    } while (ddd);

    //输出结果
    cout << "this is successive over relaxation  method  " << endl;
    for (i = 0; i<n; i++)
        cout << "x[" << i + 1 << "]= " << xk[i] << endl;
    cout << "recursive times: "<< count << endl;
    cout << endl;
}

下面是测试和输出：

int main()
{
    double a[4][4] = { -4,1,1,1,1,-4,1,1,1,1,-4,1,1,1,1,-4 };
    double x[4] = { 1,1,1,1 };
    Jacobi(a, x);
    Gausssin_Jacobi(a, x);
    cout<<endl;

    SOR(a, x,1.3);
    return 0;
}                                                                                                                                    this is JACOBI method
x[1]= -0.999994
x[2]= -0.999994
x[3]= -0.999994
x[4]= -0.999994
recursive times:  42
this is Gausssin_Jacobi method
x[1]= -0.999997
x[2]= -0.999997
x[3]= -0.999997
x[4]= -0.999998
recursive times:23

this is successive over relaxation  method
x[1]= -1
x[2]= -0.999999
x[3]= -1
x[4]= -1
recursive times: 12

对于输出，可以看到，SOR方法满足同样精度的话，迭代次数更少。在从与精确解的精度比较，显然SOR方法的精度更高，但是SOR的松弛因子的选取是个很大的问题，选好了对于计算有很大的便利。这是该算法性能的关键。