求解线性方程组的三种基本迭代法

前言

　　在实际项目的一些矩阵运算模块中，往往需要对线性方程组进行求解以得到最终结果。

　　然而，你无法让计算机去使用克莱默法则或者高斯消元法这样的纯数学方法来进行求解。

　　计算机解决这个问题的方法是迭代法。本文将介绍三种最为经典的迭代法并用经典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;
}

程序演示

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

　　


　　将项目主函数中雅克比和高斯迭代计算的部分注释掉，运行：

　　


说明

　　不是每个方程组都能迭代得到解，我们通常将系数矩阵转换为对角占优矩阵(右向量部分也要跟着转)。满足此条件的方程组的解向量大都能用这几种迭代法算出来。