《电路计算C++与MATLAB》学习笔记（五）

（一）Gauss（高斯）消去法解线性代数方程代码与理解

Gauss消去法，又叫全选主元高斯（Gauss）消去法，一般用于求解线性方程组和矩阵的秩，与之相对的是Gauss-Jordan消元法，也称全选主元高斯约旦（Gauss-Jordan）消去法，是Jordan对Gauss消去法的一种改进，一般用于求解线性方程组和矩阵的逆。

#include "stdafx.h"#include<stdlib.h>#include<iostream>using namespace std;#include<fstream>#include<math.h>#define size 10typedef double matrix[size][size];void gauss(matrix& a, int m, int&flag1)//a是实型二维数组{int i, j, k, n, rowindex; double pivot, inter;n = m + 1;for (i = 1; i <= m; i++){pivot = fabs(a[i][i]);rowindex = i;for(j=i+1;j<=m;j++)if(fabs(a[j][i])>pivot){pivot = fabs(a[j][i]); rowindex = j;}if (pivot < 1.0e-6) { flag1 = 0; break; }if(rowindex!=i)//互换两行for(k=i;k<=n;k++){inter = a[i][k]; a[i][k] = a[rowindex][k]; a[rowindex][k] = inter;}//正向消去for (j = i + 1; j <= m; j++)for (k = i + 1; k <= n; k++)a[j][k] = a[j][k] - a[j][i] * a[i][k] / a[i][i];}/*for i*///反向回代a[m][n] = a[m][n] / a[m][m];for(i=m-1;i>=1;i--){for (k = i + 1; k <= m; k++)a[i][n] = a[i][n] - a[i][k]*a[k][n];a[i][n] = a[i][n] / a[i][i];}//for i}//for gaussint main(){int m, ma, i, j, flag; matrix r, r1; fstream fs;flag = 1;fs.open("1-3-1in.txt", ios::in);//打开输入文件1-3-in.txtif (!fs) { cout << "open fail"; return (1); }fs >> m;//输入方程阶数ma = m + 1;//ma是增广矩阵的列for(i=1;i<=m;i++)for(j=1;j<=ma;j++){fs >> r[i][j];  r1[i][j] = r[i][j];//输入方程的各系数和右方项}gauss(r, m, flag);//调用gauss函数if (flag == 0) { cout << "Determinant equals zero" << endl; exit(1); }for (i = 1; i <= m; i++)cout << " " << r[i][ma];//输出解答cout << endl;double c[10];for(i=1;i<=m;i++){c[i] = 0.0;for (j = 1; j <= m; j++)c[i] = c[i] + r1[i][j] * r[j][ma];}for (i = 1; i <= m; i++)cout << " " << c[i];//输出校验结果    return 0;}

（二）问题与总结

a)、书本的代码较为混乱，主要体现在以下几点

• 采用#include"matrix.cpp",此用法错误
• 二维矩阵的表示方法有误，书本采用a(i,j)方式，因此改为a[i][j]二维数组形式

b)、Gauss 消元法介绍，可以求矩阵的秩和进行矩阵的逆变换

高斯消元法主要分为三个步骤，这个程序是建立增广矩阵C=[A|B]，将矩阵A化为上三角矩阵，再进行未知数回代的方法，从而得出方程的解

①交换两行，作用是找到主元素pivot,并使pivot尽可能大，只要不为0即可

for(j=i+1;j<=m;j++)if(fabs(a[j][i])>pivot){pivot = fabs(a[j][i]); rowindex = j;}if (pivot < 1.0e-6) { flag1 = 0; break; }if(rowindex!=i)//互换两行for(k=i;k<=n;k++){inter = a[i][k]; a[i][k] = a[rowindex][k]; a[rowindex][k] = inter;}

②正向消去，这里代码省去了主元素所在列的下面的元素消为0，但因为行变换在程序中是独立的（需要执行循环），在进行行变换时，实际就已经消去了，最后就转化为上三角矩阵

for (j = i + 1; j <= m; j++)for (k = i + 1; k <= n; k++)a[j][k] = a[j][k] - a[j][i] * a[i][k] / a[i][i];}/*for i*/

③反向回代，在化为上三角矩阵后，最后一列只有一个未知数，故可以求出来

a[m][n] = a[m][n] / a[m][m];

求出最后一列中的未知数后，再往上一列回代，又可以求出另一未知数，依次循环回代，就可以求解出列向量的所有未知数

for(i=m-1;i>=1;i--){for (k = i + 1; k <= m; k++)a[i][n] = a[i][n] - a[i][k]*a[k][n];a[i][n] = a[i][n] / a[i][i];}//for i

④补充与总结

在反向回代的过程中，有些算法是再把上三角矩阵化为单位阵，这样矩阵右边就是方程的解，实际这就是高斯约旦（Gauss-Jordan）消去法的方法，这也就是高斯（Gauss）消元法与高斯约旦（Gauss-Jordan）消去法的差异，因此高斯（Gauss）消元法也可用于求矩阵的秩，高斯约旦（Gauss-Jordan）消去法也可用于求矩阵的逆