C语言求矩阵的逆（高斯法） |
您所在的位置：网站首页 › 矩阵求逆计算器怎么用 › C语言求矩阵的逆（高斯法） |
C语言求矩阵的逆（高斯法）
|
初等变换法是常用的矩阵求逆方法之一 相对于伴随法，初等行变换法有着较低的时间复杂度，可以进行相对高维的矩阵运算，但同时也会损失一点点精度。 伴随法可参考之前的博客：C语言求矩阵的逆（伴随法） LU矩阵分解法求逆可参考：LU分解法求矩阵逆 更新： 原代码在进行增广矩阵进行内存拷贝的时，使用的是n（2*col），应改为col(原矩阵的列)，此处已做更改，并在原代码中指出。最后释放内存时需要循环释放，已做修改。_msize函数在Mac/Linux中无法使用，这里给出了另一种解决思路，见方法二。更改选主元时的bug。目录 数学原理 矩阵的初等行变换 利用增广矩阵求逆 选择主元 程序设计 法一 整体代码 测试 法二 整体代码 测试 数学原理 矩阵的初等行变换矩阵的初等变换又分为矩阵的初等行变换和矩阵的初等列变换。矩阵的初等行变换和初等列变换统称为初等变换。另外：分块矩阵也可以定义初等变换。 定义：如果B可以由A经过一系列初等变换得到，则称矩阵A与B称为等价。 数域V上的矩阵具有以下三种行变换： 1）以V中一个非零的数，乘以矩阵的某一行 2）把矩阵的某一行的a倍加到另一行(a∈V) 3）互换矩阵中两行的位置 矩阵A通过以上三种行变换变成矩阵B，称为A→B （注意：不要与行列式的行变换性质混淆） 利用增广矩阵求逆如果想求矩阵A的逆，可以通过拼上一个单位E矩阵，形成A|E的增广矩阵，通过初等行变换使其成为E|B，B就是A的逆矩阵。 例如： (1)形成A|E (2)进行初等行变换，A|E→E|B 逆矩阵B则是 主元就是在矩阵消去过程中，每列的要保留的非零元素，用它可以把该列其他消去。在阶梯型矩阵中，主元就是每个非零行第一个非零元素就是主元。 高斯消去法在消元过程中可能出现零主元，即 选择主元的方法： 1）找到主对角线以下每列最大的数Max所在的行数k 2）利用初等行变换——换行变换，将k行与当前主元行互换（记录总共换行次数n） 3）以当前主元行为基，利用初等行变换——消法变换，将主对角线下方消0 4）行列式每次换行需变号，行列式最后的符号为 5）每次进行高斯消去前都必须选择主元，计算n维的行列式，则需要进行n-1次主元选择 程序设计 法一整体流程： 1）判断传入指针是否为空，判断矩阵是否为方针 2）为增广矩阵、输出的逆矩阵开辟空间，并初始化为0 3）将原矩阵中的数据拷贝到增广矩阵中，并将增广矩阵右侧化为单位阵 4）逐列开始，选择主元，将其他位置化为0，每列阶梯位置化为1 6）将增广矩阵右侧的逆矩阵拷贝到输入矩阵中，并释放增广矩阵内存 整体代码 #include #include #include #include double** Matrix_inver(double** src) { //step 1 //判断指针是否为空 if (src == NULL)exit(-1); int i, j, k, row, col, n，principal; double** res, ** res2,tmp;//res为增广矩阵，res为输出的逆矩阵 double Max; //判断矩阵维数 row = (double)_msize(src) / (double)sizeof(double*); col = (double)_msize(*src) / (double)sizeof(double); if (row != col)exit(-1); //step 2 res = (double**)malloc(sizeof(double*) * row); res2 = (double**)malloc(sizeof(double*) * row); n = 2 * row; for (i = 0; i < row; i++) { res[i] = (double*)malloc(sizeof(double) * n); res2[i] = (double*)malloc(sizeof(double) * col); memset(res[i], 0, sizeof(res[0][0]) * n);//初始化 memset(res2[i], 0, sizeof(res2[0][0]) * col); } //step 3 //进行数据拷贝 for (i = 0; i < row; i++) { //此处源代码中的n已改为col，当然方阵的row和col相等，这里用col为了整体逻辑更加清晰 memcpy(res[i], src[i], sizeof(res[0][0]) * col); } //将增广矩阵右侧变为单位阵 for (i = 0; i < row; i++) { for (j = col; j < n; j++) { if (i == (j - row)) res[i][j] = 1.0; } } for (j = 0; j < col; j++) { //step 4 //整理增广矩阵，选主元 principal = j; Max = fabs(res.[principal][j]); // 用绝对值比较 // 默认第一行的数最大 // 主元只选主对角线下方 for (i = j; i < row; i++) { if (fabs(res.[i][j]) > Max) { principal = i; Max = fabs(res.[i][j]); } } if (j != principal) { for (k = 0; k < n; k++) { tmp = res.[principal][k]; res.[principal][k] = res.[j][k]; res.[j][k] = tmp; } } //step 5 //将每列其他元素化0 for (i = 0; i < row; i++) { if (i == j || res[i][j] == 0)continue; double b = res[i][j] / res[j][j]; for (k = 0; k < n; k++) { res[i][k] += b * res[j][k] * (-1); } } //阶梯处化成1 double a = 1.0 / res[j][j]; for (i = 0; i < n; i++) { res[j][i] *= a; } } //step 6 //将逆矩阵部分拷贝到res2中 for (i = 0; i < row; i++) { memcpy(res2[i], res[i] + row, sizeof(res[0][0]) * row); } //必须释放res内存！ for(i = 0; i < row; i++) { free(res[i]); } free(res); return res2; }上述代码中： 函数传入的参数需是以malloc开辟的动态矩阵，如固定为二维数组，需自行更改 step 1： row = (double)_msize(src) / (double)sizeof(double*); col = (double)_msize(*src) / (double)sizeof(double); 为判断矩阵的维数；其中_msize为库函数，需要包含头文件 返回指针指向内存的小心，单位为字节。step 2： 增广矩阵的行数不变，列数是原来的2倍；对照该方法开辟矩阵空间，便可以理解step1中如何计算矩阵维数。step 3： memset和memcpy都需要包含头文件需要注意memset是以字节对内存进行赋值的，对非字符（char类型）的类型进行赋值时，只使用赋0或-1。step 4： 选择主元：该步骤具有两层意义：1）相对程度上提高计算精度2）保证阶梯处的数都非0step 5： j代表从每列开始进行消0，i和k代表行和列；不要混淆j和k：j代表消零的大循环，k代表赋值的小循环。step 6： res[i] + row代表只拷贝增广矩阵中逆矩阵部分，函数结束前必须释放增广矩阵的内存，否则将造成内存泄漏。关于上述函数不太理解的地方可以参考之前发的判断矩阵维数和伴随法求逆矩阵 C语言动态内存部分可参考C语言动态内存管理 测试相关测试函数（创建矩阵，初始化矩阵，打印矩阵） double** MakeMat(int n) { int i = 0; if (n Max) { principal = i; Max = fabs(res.data[i][j]); } } if (j != principal) { for (k = 0; k < n; k++) { tmp = res.data[principal][k]; res.data[principal][k] = res.data[j][k]; res.data[j][k] = tmp; } } // step 5 // 将每列其他元素化0 for (i = 0; i < src.row; i++) { if (i == j || res.data[i][j] == 0) continue; double b = res.data[i][j] / res.data[j][j]; for (k = 0; k < n; k++) { res.data[i][k] += b * res.data[j][k] * (-1); } } // 阶梯处化成1 double a = 1.0 / res.data[j][j]; for (i = 0; i < n; i++) { res.data[j][i] *= a; } } // step 6 // 将逆矩阵部分拷贝到res2中 for (i = 0; i < src.row; i++) { memcpy(res2.data[i], res.data[i] + src.row, sizeof(res.data[0][0]) * src.row); } // 释放增广矩阵内存 free_Matrix(res); return res2; } 测试 void test() { Matrix a = MakeMatrix(3, 3); a.data[0][0] = 1; a.data[1][1] = 2; a.data[0][1] = 1; a.data[2][0] = 1; a.data[2][2] = 2; Matrix res = Matrix_Guass_inver(a); } |
【本文地址】