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

（一）矩阵求逆的计算机算法

此算法称为高斯约旦（Gauss-Jordan）消去法，也称全选主元高斯约旦（Gauss-Jordan）消去法，主要用于求解线性方程组和矩阵的逆，与之相对的是高斯（Gauss）消去法，主要用于求解线性方程组和矩阵的秩。高斯约旦（Gauss-Jordan）消去法是对Jordan对高斯（Gauss）消去法的补充。

#include "stdafx.h"#include<iostream>using namespace std;#include<fstream>#include<math.h>#define size 10//源文件采用const int size，会报错#define size2 20//源文件采用const int size2，会报错typedef double matrix4 [size][size];//数组类型定义，规定matrix4为具有[10,10]个元素的二维双精度型数组类型typedef double matrix5 [size][size2];//数组类型定义，规定matrix5为具有[10,20]个元素的二维双精度型数组类型matrix4 a, b, c;int i, j, flag;int matmul(matrix4 a, matrix4 b, int ira, int ica, int icb, double(*pc)[size])//matmul矩阵相乘函数{int i, j, k;for(i=1;i<=ira;i++)for (j = 1; j <= icb; j++){*(*(pc + i) + j) = 0.0;//初始化pc[i][j]for (k = 1; k <= ica; k++)*(*(pc + i) + j) = *(*(pc + i) + j) + a[i][k] * b[k][j];}return (0);}/*matmul*/int matinv(matrix4 a, int n, double(*pb)[size], int &pflag)//martinv是矩阵求逆函数{int kswap, i, j, k;double pivot, temp;matrix5 c;//c是增广矩阵pflag = 1;for (i = 1; i <= n; i++)for (j = 1; j <= 2 * n; j++)c[i][j] = 0.0;//对增广矩阵清零for (i = 1; i <= n; i++)for (j = 1; j <= n; j++)c[i][j] = a[i][j];//对c的左半部分赋予A的元素for (i = 1; i <= n; i++)c[i][i + n] = 1.0;//将c的右半部分做成单位阵for (k = 1; k <= n; k++){pivot = c[k][k]; kswap = k + 1;while(fabs(pivot)<1.0e-8){if (kswap > n){pflag = 0;//矩阵求逆失败cout << "Matrix inversion fails"; break;}for (j = 1; j <= n * 2; j++){temp = c[k][j]; c[k][j] = c[kswap][j]; c[kswap][j] = temp;}/*for j,交换两行*/pivot = c[k][k];kswap = kswap + 1;}/*while*/for (j = k; j <= 2 * n; j++)c[k][j] = c[k][j] / pivot;//主元素所在行除以主元素，使主元素成为1for(i=k+1;i<=n;i++)//主元素是c[k][k]{temp = c[i][k];//使主元素（除主元素外）所在列为0for (j = k; j <= 2 * n; j++)c[i][j] = c[i][j] - temp*c[k][j];}for(i=k-1;i>=1;i--){temp = c[i][k];for (j = k; j <= 2 * n; j++)c[i][j] = c[i][j] - temp*c[k][j];}/*for i*/}//for kfor (k = 1; k <= n; k++)for (j = 1; j <= 2 * n; j++)*(*(pb + k) + j) = c[k][j + n];//传递到指针数组*pb,以便输出return (0);}/*martinv*/int main(){fstream fs; int n;fs.open("1-2-5in.txt", ios::in);if (!fs) { cout << "open fail"; return (1); }fs >> n;//输入方阵阶数for (i = 1; i <= n; i++)for (j = 1; j <= n; j++)fs >> a[i][j];//输入方阵的各元素matinv(a, n, b, flag);//调用矩阵求逆函数for (i = 1; i <= n; i++){for (j = 1; j <= n; j++)cout << " " << b[i][j];cout<< "\n";}matmul(a, b, n, n, n, c);//调用矩阵相乘函数for (i = 1; i <= n; i++){for (j = 1; j <= n; j++)cout<< " " <<c[i][j];cout<< "\n";} return 0;}

（二）代码理解

矩阵求逆的思路是建立一个增广矩阵C，C的左半部分是矩阵A，右半部分是单位阵I，通过矩阵初等运算，把C的左半部分A化成单位阵，则右半部分就是矩阵A的逆矩阵。初等运算主要包括：

• 交换矩阵的两行，关键是找到主元pivot，通过交换使得C[k][j]不等于0
• 用一个常数乘矩阵的某行，a、是将每行的主元pivot化为1，以便向单位阵转化，b、为下一步相加减铺垫
• 将矩阵的两行相加减，为了使主元pivot所在列（除主元外）的元素变为0

a)、交换矩阵的两行

for (k = 1; k <= n; k++){pivot = c[k][k]; kswap = k + 1;while(fabs(pivot)<1.0e-8){if (kswap > n){pflag = 0;//矩阵求逆失败cout << "Matrix inversion fails"; break;}for (j = 1; j <= n * 2; j++){temp = c[k][j]; c[k][j] = c[kswap][j]; c[kswap][j] = temp;//交换两行                        }       pivot = c[k][k];kswap = kswap + 1;}/*while*/

交换两行的作用是找到不为零的主元pivot，始终按顺序使第1行至第n行的pivot不为0，while为真说明c[k][k]为0（fabs函数为取绝对值函数），假定c[1][1]为0，则调用while循环语句，与第2行进行交换，检查与第二行交换后的pivot，若仍为0，则继续执行while，与第3行交换，一直交换循环直至找到不为0的pivot，即与第1列中最先不为0的元素所在的行进行交换，kswap > n 的作用是检查c[n][n]主元，如果（n-1）行的所有交换、相乘、相加减都处理完毕，主元c[n][n]为0，相当于行列式值为0，则矩阵的逆不存在

b)、用一个常数去乘以矩阵的某行，本例为主元素所在行除以主元素，使主元素为1

for (j = k; j <= 2 * n; j++)c[k][j] = c[k][j] / pivot;//主元素所在行除以主元素，使主元素成为1for (k = 1; k <= n; k++)//大的for循环

主元素所在行都除以主元素，使pivot c[k][k]=1，此外主元素所在行的其他元素都要进行相应的除法运算，需要说明的是因为最外面的for大循环存在，算法的思路是一行一行处理

c)、将矩阵的两行相加减，为了使主元pivot所在列（除主元外）的元素变为0

for(i=k+1;i<=n;i++)//主元素是c[k][k],以k为界，对k+1行进行变换{temp = c[i][k];//使主元素（除主元素外）所在列为0for (j = k; j <= 2 * n; j++)c[i][j] = c[i][j] - temp*c[k][j];}for(i=k-1;i>=1;i--)//主元素是c[k][k],以k为界，对k-1行进行变换{temp = c[i][k];for (j = k; j <= 2 * n; j++)c[i][j] = c[i][j] - temp*c[k][j];}/*for i*/

因为主元素为1，通过for循环使主元素所在列的其他元素变为0，以pivot c[1][1]为例，k=1，k-1之前的变换不存在，所以相当于从第2行到第n行变换，使第2行到第n行的首元素变为0，至此c[1][1]的变换暂时完成，再进行pivot c[2][2]变换时，k-1=1的变换存在，所以除使第3行到第n行的第2列元素为0外，还要进行第1行的消0变换，这就是for(i=k-1)的作用，至此通过大循环就能把C的左半部分化为单位阵，则右半部分就是A的逆阵

（三）问题

#define size 10//源文件采用const int size，在visual studio 2017下会报错#define size2 20//源文件采用const int size2，在visual studio 2017下会报错

a)、此外定义函数时，需要指定函数返回值类型，int,double等
b)、当需要输入的数据过多时，可以采用读取.txt文件的方式进行数据的写入。本例中“1-2-5in.txt”文件与.cpp放在同一个文件夹下面，debug运行.cpp文件时会查找到此txt文件

fstream fs; int n;fs.open("1-2-5in.txt", ios::in);if (!fs) { cout << "open fail"; return (1); }