本系列是數值分析相關演算法的文章，這次用迭代法求解線性方程組，不同於上次用高斯消元法之類的求解。迭代法對於稀疏矩陣方程組的運算，會大大提高。而如果用高斯相關的演算法求解，會浪費大量資源計算無用的東西，所以有必要研究此演算法。
本文章主要使用了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的鬆弛因子的選取是個很大的問題，選好了對於計算有很大的便利。這是該演算法效能的關鍵。