求解線性方程組的三種基本迭代法

前言

　　在實際項目的一些矩陣運算模塊中，往往需要對線性方程組進行求解以得到最終結果。

　　然而，你無法讓計算機去使用克萊默法則或者高斯消元法這樣的純數學方法來進行求解。

　　計算機解決這個問題的方法是迭代法。本文將介紹三種最為經典的迭代法並用經典C++源代碼實現之。

迭代法簡介

　　從解的某個近似值出發，通過構造一個無窮序列去逼近精確解的方法。

雅克比迭代法

　　計算流程：

　　　　1. 初始化系數矩陣等計算環境

　　　　2. 設置精度控制和迭代次數控制變量

　　　　3. 采用如下式子進行迭代計算：

　　　　　　

　　　　4. 循環執行 3，若(條件a)當前滿足迭代精度控制的解分量數等於解向量維數或者(條件b)迭代次數達到最大次數則跳出循環，進入到 5。

　　　　5. 若是因為條件a跳出循環則迭代成功，打印解向量；若是因為條件b退出循環則表示在指定的迭代次數內方程組未求解完。

　　說明：最經典原始的一種解方程組的迭代法。

高斯迭代法

　　計算流程：

　　　　1. 初始化系數矩陣等計算環境

　　　　2. 設置精度控制和迭代次數控制變量

　　　　3. 采用如下式子進行迭代計算：

　　　　　　

　　　　4. 循環執行 3，若(條件a)當前滿足迭代精度控制的解分量數等於解向量維數或者(條件b)迭代次數達到最大次數則跳出循環，進入到 5。

　　　　5. 若是因為條件a跳出循環則迭代成功，打印解向量；若是因為條件b退出循環則表示在指定的迭代次數內方程組未求解完。

　　說明：相對於雅克比迭代法，該方法不需要一個額外的輔助向量保存上一次迭代計算的結果，節省了空間。

超松弛迭代法

　　計算流程：

　　　　1. 初始化系數矩陣等計算環境

　　　　2. 設置精度控制和迭代次數控制變量以及松弛因子omiga

　　　　3. 采用如下式子進行迭代計算：

　　　　　　

　　　　4. 循環執行 3，若(條件a)當前滿足迭代精度控制的解分量數等於解向量維數或者(條件b)迭代次數達到最大次數則跳出循環，進入到 5。

　　　　5. 若是因為條件a跳出循環則迭代成功，打印解向量；若是因為條件b退出循環則表示在指定的迭代次數內方程組未求解完。

　　說明：

　　　　1. 該方法同樣不需要一個額外的輔助向量保存上一次迭代計算的結果。

　　　　2. 需要設置一個w因子參數，一般取0-2。當小於1時該方法為低松弛迭代法，高於1時為超松弛迭代法。經驗上來看一般取1.4-1.6來實現超松弛迭代。

源代碼實現

　　第一部分：迭代類聲明 (Iterator.h)

// 頭文件保護符
#ifndef Iterator_H
#define Iterator_H

// 迭代計算類
class CIterator
{
public:
    // 構造函數
    CIterator ();
    // 析構函數
    ~CIterator ();
    // 設置計算環境 (如系數矩陣等)
    bool SetEnvironment ();
    // 雅克比迭代法計算方程解
    bool JacobianCalculate ();
    // 高斯迭代法計算方程解
    bool GaussCalculate ();
    // 超松弛迭代法計算方程解
    bool RelaxationCalculate (double omiga);
    // 獲取系數矩陣
    double ** GetMatrixA ();
    // 獲取系數矩陣的階
    int GetOrder ();
    // 獲取方程組右值向量
    double * GetVectorB ();
    // 獲取方程解向量
    double * GetSolution ();

private:
    double **matrixA;    // 系數矩陣
    int order;            // 系數矩陣的階
    double *vectorB;    // 方程組右值向量
    double *solution;    // 解向量
};

#endif

　　第二部分：迭代類實現 (Iterator.cpp)

#include "Iterator.h"
#include <iostream>
#include <iomanip>

using namespace std;

// 構造函數
CIterator :: CIterator ()
{
    matrixA = NULL;
    order = 0;
    vectorB = NULL;
    solution = NULL;
}

// 析構函數
CIterator :: ~CIterator ()
{
    // 釋放內存空間
    if (!vectorB) {
        delete [] vectorB;
        vectorB = NULL;
    }
    if (!solution) {
        delete [] solution;
        solution = NULL;
    }
    if (matrixA!=NULL) {
        for (int i=0; i<order; i++) {
            delete [] matrixA[i];
            matrixA[i] = NULL;
        }
        delete [] matrixA;
        matrixA = NULL;
    }
}

// 設置計算環境 (如系數矩陣等)
bool CIterator :: SetEnvironment ()
{
    cout << "系數矩陣階數: ";
    cin >> order;
    cout << endl;

    matrixA = new double *[order];
    for (int i=0; i<order; i++)
        matrixA[i] = new double [order];

    for (int i=0; i<order; i++) {
        cout << "第 " << i << " 行系數: ";
        for (int j=0; j<order; j++)
            cin >> matrixA[i][j];
    }
    cout << endl;

    vectorB = new double [order];
    cout << "b 向量: ";
    for (int i=0; i<order; i++)
            cin >> vectorB[i];
    cout << endl;

    solution = new double [order];

    cout << "運算環境搭建完畢" << endl << endl;

    return true;
}

// 雅克比迭代法計算方程解
bool CIterator :: JacobianCalculate ()
{
    cout << "下面使用雅克比迭代法計算該線性方程組" << endl << endl;

    // 設置迭代精度控制值
    cout << "迭代精度控制值: ";
    double bias;
    cin >> bias;
    cout << endl;

    // 設置迭代次數控制值
    cout << "迭代次數控制值: ";
    int itMax;
    cin >> itMax;
    cout << endl;

    // 當前滿足迭代精度控制的解分量數
    int meetPrecisionCount = 0;

    // 輔助向量，存放上一次迭代的解向量。
    double * auxiliary = new double [order];

    // 初始化解向量
    for (int i=0; i<order; i++) {
        auxiliary[i] = 0;
        solution[i] = 1;
    }

    // 當前迭代次數
    int itCur = 0;

    // 若當前滿足迭代精度控制的解分量數等於解向量維數或者迭代次數達到最大次數則跳出循環
    while (meetPrecisionCount<order && itCur<itMax) {
        
        // 當前滿足迭代精度控制的解分量數清零
        meetPrecisionCount = 0;

        // 將當前解向量存入輔助向量
        for (int i=0; i<order; i++)
            auxiliary[i] = solution[i];

        // 計算新的解向量
        for (int i=0; i<order; i++) {

            // 暫存當前解分量
            double temp = solution[i];

            // 清零當前解分量
            solution[i] = 0;

            // 雅克比迭代計算
            for (int j=0; j<i; j++) {
                solution[i] += matrixA[i][j]*auxiliary[j];
            }
            for (int j=i+1; j<order; j++) {
                solution[i] += matrixA[i][j]*auxiliary[j];
            }
            solution[i] = (vectorB[i]-solution[i])/matrixA[i][i];

            // 更新當前滿足迭代精度控制的解分量數
            if (fabs(temp-solution[i])<bias)
                meetPrecisionCount++;
        }

        // 當前迭代次數遞增
        itCur++;
    }

    // 若在規定的迭代次數內未完成計算任務，則返回錯誤。
    if (itCur == itMax) return false;

    return true;
}

// 高斯迭代法計算方程解
bool CIterator :: GaussCalculate ()
{
    cout << "下面使用高斯迭代法計算該線性方程組" << endl << endl;

    // 設置迭代精度控制值
    cout << "迭代精度控制值: ";
    double bias;
    cin >> bias;
    cout << endl;

    // 設置迭代次數控制值
    cout << "迭代次數控制值: ";
    int itMax;
    cin >> itMax;
    cout << endl;

    // 當前滿足迭代精度控制的解分量數
    int meetPrecisionCount = 0;

    // 初始化解向量
    for (int i=0; i<order; i++) {
        solution[i] = 0;
    }

    // 當前迭代次數
    int itCur = 0;

    // 若當前滿足迭代精度控制的解分量數等於解向量維數或者迭代次數達到最大次數則跳出循環
    while (meetPrecisionCount<order && itCur<itMax) {
        
        // 當前滿足迭代精度控制的解分量數清零
        meetPrecisionCount = 0;

        // 計算新的解向量
        for (int i=0; i<order; i++) {

            // 暫存當前解分量
            double temp = solution[i];

            // 清零當前解分量
            solution[i] = 0;

            // 高斯迭代計算
            for (int j=0; j<i; j++) {
                solution[i] += matrixA[i][j]*solution[j];
            }
            for (int j=i+1; j<order; j++) {
                solution[i] += matrixA[i][j]*solution[j];
            }
            solution[i] = (vectorB[i]-solution[i])/matrixA[i][i];
            
            // 更新當前滿足迭代精度控制的解分量數
            if (fabs(temp-solution[i])<bias)
                meetPrecisionCount++;
        }

        // 當前迭代次數遞增
        itCur++;
    }

    // 若在規定的迭代次數內未完成計算任務，則返回錯誤。
    if (itCur == itMax) return false;

    return true;
}

// 超松弛迭代法計算方程解
bool CIterator :: RelaxationCalculate (double omiga)
{
    cout << "下面使用超松弛迭代法計算該線性方程組" << endl << endl;

    // 設置迭代精度控制值
    cout << "迭代精度控制值: ";
    double bias;
    cin >> bias;
    cout << endl;

    // 設置迭代次數控制值
    cout << "迭代次數控制值: ";
    int itMax;
    cin >> itMax;
    cout << endl;

    // 當前滿足迭代精度控制的解分量數
    int meetPrecisionCount = 0;

    // 當前滿足迭代精度控制的解分量數
    for (int i=0; i<order; i++) {
        solution[i] = 0;
    }

    // 當前迭代次數
    int itCur = 0;

    // 若當前滿足迭代精度控制的解分量數等於解向量維數或者迭代次數達到最大次數則跳出循環
    while (meetPrecisionCount<order && itCur<itMax) {
        
        // 當前滿足迭代精度控制的解分量數清零
        meetPrecisionCount = 0;

        // 計算新的解向量
        for (int i=0; i<order; i++) {

            // 暫存當前解分量
            double temp = solution[i];

            // 清零當前解分量
            solution[i] = 0;

            // 超松弛迭代計算
            for (int j=0; j<i; j++) {
                solution[i] += matrixA[i][j]*solution[j];
            }
            for (int j=i+1; j<order; j++) {
                solution[i] += matrixA[i][j]*solution[j];
            }
            solution[i] = omiga*(vectorB[i]-solution[i])/matrixA[i][i] + (1-omiga)*temp;
            
            // 更新當前滿足迭代精度控制的解分量數
            if (fabs(temp-solution[i])<bias)
                meetPrecisionCount++;
        }

        // 當前迭代次數遞增
        itCur++;
    }

    // 若在規定的迭代次數內未完成計算任務，則返回錯誤。
    if (itCur == itMax) return false;

    return true;
}

// 獲取系數矩陣
double ** CIterator :: GetMatrixA ()
{
    return this->matrixA;
}

// 獲取系數矩陣的階
int CIterator :: GetOrder ()
{
    return this->order;
}

// 獲取方程組右值向量
double * CIterator :: GetVectorB ()
{
    return this->vectorB;
}

// 獲取方程解向量
double * CIterator :: GetSolution ()
{
    return this->solution;
}

　　第三部分：主函數 (main.cpp)

// 迭代計算類
#include "Iterator.h"
#include <iostream>

using namespace std;

// 打印方程解
void printResults (CIterator iterator);

int main()
{
    // 定義迭代類對象
    CIterator iterator;

    // 設置計算環境 (如系數矩陣等)
    iterator.SetEnvironment();
    
    // 雅克比迭代法計算方程解
    if (!iterator.JacobianCalculate()) {
        cout << "規定次數內未完成迭代計算" << endl;
        return EXIT_FAILURE;
    }

    // 高斯迭代法計算方程解
    /*
    if (!iterator.GaussCalculate()) {
        cout << "規定次數內未完成迭代計算" << endl;
        return EXIT_FAILURE;
    }
    */

    // 超松弛迭代法計算方程解
    /*
    double omiga = 1.5;        // 松弛迭代系數
    if (!iterator.RelaxationCalculate(omiga)) {
        cout << "規定次數內未完成迭代計算" << endl;
        return EXIT_FAILURE;
    }
    */

    // 打印方程解
    printResults (iterator);

    return EXIT_SUCCESS;
}

// 打印方程解
void printResults (CIterator iterator)
{
    cout << "計算成功，打印方程解:  " << endl;
    for (int i=0; i<iterator.GetOrder(); i++)
        cout << "x" << i << " = " << iterator.GetSolution()[i] << endl;

    cout << endl;
}

程序演示

　　使用超松弛迭代法解如下方程組：

　　

　　將項目主函數中雅克比和高斯迭代計算的部分注釋掉，運行：

　　

說明

　　不是每個方程組都能迭代得到解，我們通常將系數矩陣轉換為對角占優矩陣(右向量部分也要跟着轉)。滿足此條件的方程組的解向量大都能用這幾種迭代法算出來。