首页 > 代码库 > 求解线性方程组的三种基本迭代法
求解线性方程组的三种基本迭代法
前言
在实际项目的一些矩阵运算模块中，往往需要对线性方程组进行求解以得到最终结果。
然而，你无法让计算机去使用克莱默法则或者高斯消元法这样的纯数学方法来进行求解。
计算机解决这个问题的方法是迭代法。本文将介绍三种最为经典的迭代法并用经典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. 设置精度控制和迭代次数控制变量以及松弛因子w
3. 采用如下式子进行迭代计算：
4. 循环执行 3，若(条件a)当前满足迭代精度控制的解分量数等于解向量维数或者(b)迭代次数达到最大次数则跳出循环，进入到 5。
5. 若是因为条件a跳出循环则迭代成功，打印解向量；若是因为条件b退出循环则表示在指定的迭代次数内方程组未求解完。
说明：
1. 该方法同样不需要一个额外的辅助向量保存上一次迭代计算的结果。
2. 需要设置一个w因子参数，一般取0-2。当小于1时该方法为低松弛迭代法，高于1时为超松弛迭代法。经验上来看一般取1.4-1.6来实现超松弛迭代。
源代码实现
第一部分：迭代类声明 (Iterator.h)
1 // 头文件保护符 2 #ifndef Iterator_H 3 #define Iterator_H 4 5 // 迭代计算类 6 class CIterator 7 { 8 public: 9 // 构造函数10 CIterator ();11 // 析构函数12 ~CIterator ();13 // 设置计算环境 (如系数矩阵等)14 bool SetEnvironment ();15 // 雅克比迭代法计算方程解16 bool JacobianCalculate ();17 // 高斯迭代法计算方程解18 bool GaussCalculate ();19 // 超松弛迭代法计算方程解20 bool RelaxationCalculate (double omiga);21 // 获取系数矩阵22 double ** GetMatrixA ();23 // 获取系数矩阵的阶24 int GetOrder ();25 // 获取方程组右值向量26 double * GetVectorB ();27 // 获取方程解向量28 double * GetSolution ();29 30 private:31 double **matrixA; // 系数矩阵32 int order; // 系数矩阵的阶33 double *vectorB; // 方程组右值向量34 double *solution; // 解向量35 };36 37 #endif
第二部分：迭代类实现 (Iterator.cpp)
1 #include "Iterator.h" 2 #include <iostream> 3 #include <iomanip> 4 5 using namespace std; 6 7 // 构造函数 8 CIterator :: CIterator () 9 { 10 matrixA = NULL; 11 order = 0; 12 vectorB = NULL; 13 solution = NULL; 14 } 15 16 // 析构函数 17 CIterator :: ~CIterator () 18 { 19 // 释放内存空间 20 if (!vectorB) { 21 delete [] vectorB; 22 vectorB = NULL; 23 } 24 if (!solution) { 25 delete [] solution; 26 solution = NULL; 27 } 28 if (matrixA!=NULL) { 29 for (int i=0; i<order; i++) { 30 delete [] matrixA[i]; 31 matrixA[i] = NULL; 32 } 33 delete [] matrixA; 34 matrixA = NULL; 35 } 36 } 37 38 // 设置计算环境 (如系数矩阵等) 39 bool CIterator :: SetEnvironment () 40 { 41 cout << "系数矩阵阶数: "; 42 cin >> order; 43 cout << endl; 44 45 matrixA = new double *[order]; 46 for (int i=0; i<order; i++) 47 matrixA[i] = new double [order]; 48 49 for (int i=0; i<order; i++) { 50 cout << "第 " << i << " 行系数: "; 51 for (int j=0; j<order; j++) 52 cin >> matrixA[i][j]; 53 } 54 cout << endl; 55 56 vectorB = new double [order]; 57 cout << "b 向量: "; 58 for (int i=0; i<order; i++) 59 cin >> vectorB[i]; 60 cout << endl; 61 62 solution = new double [order]; 63 64 cout << "运算环境搭建完毕" << endl << endl; 65 66 return true; 67 } 68 69 // 雅克比迭代法计算方程解 70 bool CIterator :: JacobianCalculate () 71 { 72 cout << "下面使用雅克比迭代法计算该线性方程组" << endl << endl; 73 74 // 设置迭代精度控制值 75 cout << "迭代精度控制值: "; 76 double bias; 77 cin >> bias; 78 cout << endl; 79 80 // 设置迭代次数控制值 81 cout << "迭代次数控制值: "; 82 int itMax; 83 cin >> itMax; 84 cout << endl; 85 86 // 当前满足迭代精度控制的解分量数 87 int meetPrecisionCount = 0; 88 89 // 辅助向量，存放上一次迭代的解向量。 90 double * auxiliary = new double [order]; 91 92 // 初始化解向量 93 for (int i=0; i<order; i++) { 94 auxiliary[i] = 0; 95 solution[i] = 1; 96 } 97 98 // 当前迭代次数 99 int itCur = 0;100 101 // 若当前满足迭代精度控制的解分量数等于解向量维数或者迭代次数达到最大次数则跳出循环102 while (meetPrecisionCount<order && itCur<itMax) {103 104 // 当前满足迭代精度控制的解分量数清零105 meetPrecisionCount = 0;106 107 // 将当前解向量存入辅助向量108 for (int i=0; i<order; i++)109 auxiliary[i] = solution[i];110 111 // 计算新的解向量112 for (int i=0; i<order; i++) {113 114 // 暂存当前解分量115 double temp = solution[i];116 117 // 清零当前解分量118 solution[i] = 0;119 120 // 雅克比迭代计算121 for (int j=0; j<i; j++) {122 solution[i] += matrixA[i][j]*auxiliary[j];123 }124 for (int j=i+1; j<order; j++) {125 solution[i] += matrixA[i][j]*auxiliary[j];126 }127 solution[i] = (vectorB[i]-solution[i])/matrixA[i][i];128 129 // 更新当前满足迭代精度控制的解分量数130 if (fabs(temp-solution[i])<bias)131 meetPrecisionCount++;132 }133 134 // 当前迭代次数递增135 itCur++;136 }137 138 // 若在规定的迭代次数内未完成计算任务，则返回错误。139 if (itCur == itMax) return false;140 141 return true;142 }143 144 // 高斯迭代法计算方程解145 bool CIterator :: GaussCalculate ()146 {147 cout << "下面使用高斯迭代法计算该线性方程组" << endl << endl;148 149 // 设置迭代精度控制值150 cout << "迭代精度控制值: ";151 double bias;152 cin >> bias;153 cout << endl;154 155 // 设置迭代次数控制值156 cout << "迭代次数控制值: ";157 int itMax;158 cin >> itMax;159 cout << endl;160 161 // 当前满足迭代精度控制的解分量数162 int meetPrecisionCount = 0;163 164 // 初始化解向量165 for (int i=0; i<order; i++) {166 solution[i] = 0;167 }168 169 // 当前迭代次数170 int itCur = 0;171 172 // 若当前满足迭代精度控制的解分量数等于解向量维数或者迭代次数达到最大次数则跳出循环173 while (meetPrecisionCount<order && itCur<itMax) {174 175 // 当前满足迭代精度控制的解分量数清零176 meetPrecisionCount = 0;177 178 // 计算新的解向量179 for (int i=0; i<order; i++) {180 181 // 暂存当前解分量182 double temp = solution[i];183 184 // 清零当前解分量185 solution[i] = 0;186 187 // 高斯迭代计算188 for (int j=0; j<i; j++) {189 solution[i] += matrixA[i][j]*solution[j];190 }191 for (int j=i+1; j<order; j++) {192 solution[i] += matrixA[i][j]*solution[j];193 }194 solution[i] = (vectorB[i]-solution[i])/matrixA[i][i];195 196 // 更新当前满足迭代精度控制的解分量数197 if (fabs(temp-solution[i])<bias)198 meetPrecisionCount++;199 }200 201 // 当前迭代次数递增202 itCur++;203 }204 205 // 若在规定的迭代次数内未完成计算任务，则返回错误。206 if (itCur == itMax) return false;207 208 return true;209 }210 211 // 超松弛迭代法计算方程解212 bool CIterator :: RelaxationCalculate (double omiga)213 {214 cout << "下面使用超松弛迭代法计算该线性方程组" << endl << endl;215 216 // 设置迭代精度控制值217 cout << "迭代精度控制值: ";218 double bias;219 cin >> bias;220 cout << endl;221 222 // 设置迭代次数控制值223 cout << "迭代次数控制值: ";224 int itMax;225 cin >> itMax;226 cout << endl;227 228 // 当前满足迭代精度控制的解分量数229 int meetPrecisionCount = 0;230 231 // 当前满足迭代精度控制的解分量数232 for (int i=0; i<order; i++) {233 solution[i] = 0;234 }235 236 // 当前迭代次数237 int itCur = 0;238 239 // 若当前满足迭代精度控制的解分量数等于解向量维数或者迭代次数达到最大次数则跳出循环240 while (meetPrecisionCount<order && itCur<itMax) {241 242 // 当前满足迭代精度控制的解分量数清零243 meetPrecisionCount = 0;244 245 // 计算新的解向量246 for (int i=0; i<order; i++) {247 248 // 暂存当前解分量249 double temp = solution[i];250 251 // 清零当前解分量252 solution[i] = 0;253 254 // 超松弛迭代计算255 for (int j=0; j<i; j++) {256 solution[i] += matrixA[i][j]*solution[j];257 }258 for (int j=i+1; j<order; j++) {259 solution[i] += matrixA[i][j]*solution[j];260 }261 solution[i] = omiga*(vectorB[i]-solution[i])/matrixA[i][i] + (1-omiga)*temp;262 263 // 更新当前满足迭代精度控制的解分量数264 if (fabs(temp-solution[i])<bias)265 meetPrecisionCount++;266 }267 268 // 当前迭代次数递增269 itCur++;270 }271 272 // 若在规定的迭代次数内未完成计算任务，则返回错误。273 if (itCur == itMax) return false;274 275 return true;276 }277 278 // 获取系数矩阵279 double ** CIterator :: GetMatrixA ()280 {281 return this->matrixA;282 }283 284 // 获取系数矩阵的阶285 int CIterator :: GetOrder ()286 {287 return this->order;288 }289 290 // 获取方程组右值向量291 double * CIterator :: GetVectorB ()292 {293 return this->vectorB;294 }295 296 // 获取方程解向量297 double * CIterator :: GetSolution ()298 {299 return this->solution;300 }
第三部分：主函数 (main.cpp)
1 // 迭代计算类 2 #include "Iterator.h" 3 #include <iostream> 4 5 using namespace std; 6 7 // 打印方程解 8 void printResults (CIterator iterator); 9 10 int main()11 {12 // 定义迭代类对象13 CIterator iterator;14 15 // 设置计算环境 (如系数矩阵等)16 iterator.SetEnvironment();17 18 // 雅克比迭代法计算方程解19 if (!iterator.JacobianCalculate()) {20 cout << "规定次数内未完成迭代计算" << endl;21 return EXIT_FAILURE;22 }23 24 // 高斯迭代法计算方程解25 /*26 if (!iterator.GaussCalculate()) {27 cout << "规定次数内未完成迭代计算" << endl;28 return EXIT_FAILURE;29 }30 */31 32 // 超松弛迭代法计算方程解33 /*34 double omiga = 1.5; // 松弛迭代系数35 if (!iterator.RelaxationCalculate(omiga)) {36 cout << "规定次数内未完成迭代计算" << endl;37 return EXIT_FAILURE;38 }39 */40 41 // 打印方程解42 printResults (iterator);43 44 return EXIT_SUCCESS;45 }46 47 // 打印方程解48 void printResults (CIterator iterator)49 {50 cout << "计算成功，打印方程解: " << endl;51 for (int i=0; i<iterator.GetOrder(); i++)52 cout << "x" << i << " = " << iterator.GetSolution()[i] << endl;53 54 cout << endl;55 }
程序演示
使用超松弛迭代法解如下方程组：
将主函数中雅克比和高斯迭代计算的部分注释掉，运行：
说明
不是每个方程组都能迭代得到解，我们通常将系数矩阵转换为对角占优矩阵(右向量部分也要跟着转)。满足这种条件的方程组的解向量大都能用这几种迭代法算出来。
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。