1、概述
给定任意n阶矩阵B,如果满秩,利用初等行变换解算矩阵B的逆矩阵B_inverse。
2、思路
1)利用初等行变换,那么要将单位矩阵E和n阶矩阵B合并(规定为EandB_normal[ n, 2 * n])
2)将EandB_normal[ n, 2 * n]转为右半部分为上三角的矩阵
>>>这一步转换比较复杂一点,具体实现就是:
>>>第一层循环,循环变量 j 从第n列开始到第2 * n - 1列结束,目的就是将该列值都转为1,方便后边变为上三角矩阵(需要注意的是,对于第n列,应该考虑把每个值都变为1;但是到第n + 1列时,就不考虑第一个值了;第n + 2列时,不考虑第一个和第二个值;类推);
>>>第二层循环,循环变量 i 从第j - n行开始到第n - 1行结束,目的是对每一行都进行除以EandB_normal[ i, j]值的运算,这样EandB_normal[ i, j]的值就变为了1(需要注意的是,如果EandB_normal[ i, j]的值为0的话,我们考虑将该行与最后一行调换,同时循环变量 i 到第n - 2行结束;如果调换之后,EandB_normal[ i, j]的值仍然为0,那么再将该行与此时的最后一行调换,类推;但是如果一直调换,直到发现始终为0,就说明矩阵B不满秩,退出计算;如果EandB_normal[ i, j]值为负数,该行同时变号);
>>>第三层循环,循环变量 k 从第0列开始到第2 * n - 1列结束,目的是将上一步中循环到的行中的每一个值都除以EandB_normal[ i, j]的值;
>>>循环全部完成之后,矩阵EandB_normal[ n, 2 * n]就变成了右半部分为上三角的矩阵。
3)接上一步,将该矩阵转为右半部分为单位矩阵的矩阵,此时即为矩阵B的逆矩阵与单位矩阵的合并(规定为B_inverse_andE[ n, 2 * n])
>>>这一步中的循环变量是递减的,具体实现就是:
>>>第一层循环,循环变量 j 从第2 * n - 1列开始到第n列结束,目的是将该列值只保留一个1,其余变为0;
>>>第二层循环,循环变量 i 从第 j - n行开始到第0行结束;
>>>第三层循环,循环变量 k 从第0列开始到第2 * n - 1列结束;拿 j = 2 * n - 1, i = n - 1举例,此时,我们希望第n - 2行的值都加上该行最后一个值的相反数与第n - 1行乘积的对应值,第n - 3行的值都加上该行最后一个值得相反数与第n - 1行乘积的对应值,类推;(需要注意的是,j = 2 * n - 2时,i从第n - 2行开始循环,j = 2 * n - 3时,i从第n - 2行开始循环,类推);
>>>当循环全部完成之后,B_inverse_andE[ n, 2 * n]的右半部分就变为了单位矩阵,左半部分为矩阵B的逆矩阵。
4)接上一步,将B的逆矩阵取出来(规定为B_inverse[n, n])
3、代码
1 //求逆 2 public static double[,] Inverse(double [,] matrixB,int orderNum) 3 { 4 //判断是否满秩 5 bool IsFullRank = true; 6 //n为阶级 7 int n = orderNum; 8 9 10 //####赋值#### 11 //矩阵B 12 //矩阵B的逆矩阵 13 //单位矩阵E和矩阵B组成的矩阵 14 double[,] B_normal = matrixB; 15 double[,] B_inverse = new double[n, n]; 16 double[,] EandB_normal = new double[n, 2 * n]; 17 for(int i = 0; i < n; i++) 18 { 19 for(int j = 0; j < n; j++) 20 { 21 if(i==j) 22 EandB_normal[i, j] = 1; 23 else 24 EandB_normal[i, j] = 0; 25 26 } 27 for(int k = n; k < 2*n; k++) 28 { 29 EandB_normal[i, k] = B_normal[i, k - n]; 30 } 31 } 32 33 34 35 //####计算#### 36 //中间变量数组,与最后一行值进行调换的时候用到 37 double[] rowHaveZero = new double[2 * n]; 38 //EB矩阵右边的n*n变为上三角矩阵 39 for(int j = n; j < 2*n; j++) 40 { 41 int firstRowN = j - n; 42 int rowCount = n; 43 int colCount = 2 * n; 44 //把EB中索引为j的列的值化为1 45 for (int i = firstRowN; i < rowCount; i++) 46 { 47 //如果为0,就把0所在的行与最后一行调换位置 48 //并且总行数减去1,直到不为0 49 double EBijNum = EandB_normal[i, j]; 50 while(EBijNum == 0 && rowCount > i + 1) 51 { 52 for (int k = 0; k < colCount; k++) 53 { 54 rowHaveZero[k] = EandB_normal[i, k]; 55 EandB_normal[i, k] = EandB_normal[rowCount - 1, k]; 56 EandB_normal[rowCount - 1, k] = rowHaveZero[k]; 57 } 58 rowCount -= 1; 59 EBijNum = EandB_normal[i, j]; 60 } 61 //如果while循环是由第二个判断跳出 62 //即EBijNum始终为0 63 if (EBijNum == 0) 64 { 65 //总行数再减去1,然后跳出循环 66 rowCount -= 1; 67 break; 68 } 69 //如果为负数,该行变号 70 if(EBijNum < 0) 71 { 72 for(int k = 0; k < colCount; k++) 73 { 74 EandB_normal[i, k] = (-1) * EandB_normal[i, k]; 75 } 76 EBijNum = EandB_normal[i, j]; 77 } 78 //将该值变为1,则其余值都除以EBijNum 79 for(int k = 0; k < colCount; k++) 80 { 81 EandB_normal[i, k] = EandB_normal[i, k] / EBijNum; 82 } 83 } 84 85 //自n列起,每列只保留一个1,呈阶梯状 86 int secondRowN = firstRowN + 1; 87 for(int i = secondRowN; i < rowCount; i++) 88 { 89 for(int k = 0; k < colCount; k++) 90 { 91 EandB_normal[i, k] = EandB_normal[i, k]
- EandB_normal[firstRowN, k]; 92 } 93 } 94 95 if(rowCount == firstRowN) 96 { 97 //矩阵不满秩 98 IsFullRank = false; 99 break; 100 } 101 } 102 //不满秩,结束运算 103 if (!IsFullRank) 104 { 105 double[,] error = new double[n, n]; 106 for(int i = 0; i < n; i++) 107 { 108 for(int j = 0; j < n; j++) 109 { 110 error[i, j] = 0; 111 } 112 } 113 //返还值均为0的矩阵 114 return error; 115 116 } 117 118 //将上三角矩阵变为单位矩阵 119 for(int j = 2*n - 1; j > n; j--) 120 { 121 int firstRowN = j - n; 122 int secondRowN = firstRowN - 1; 123 int colCount = j + 1; 124 //从最后一列起,每列只保留一个1,其余减为0 125 for(int i = secondRowN; i > -1; i--) 126 { 127 double EBijNum = EandB_normal[i, j]; 128 for(int k = 0; k < colCount; k++) 129 { 130 EandB_normal[i, k] = EandB_normal[i, k]
- EandB_normal[firstRowN, k] * EBijNum; 131 } 132 } 133 134 } 135 136 137 138 //####提取逆矩阵#### 139 for(int i = 0; i < n; i++) 140 { 141 for(int j = 0; j < n; j++) 142 { 143 B_inverse[i, j] = EandB_normal[i, j]; 144 } 145 } 146 147 return B_inverse; 148 }
4、测试运行
1)测试5阶矩阵
1 int orderNum = 5; 2 double[,] B = new double[orderNum, orderNum]; 3 B[0, 0] = 0; B[0, 1] = 0; B[0, 2] = 0; B[0, 3] = 1; B[0, 4] = 3; 4 B[1, 0] = 0; B[1, 1] = 0; B[1, 2] = 0; B[1, 3] = -1; B[1, 4] = 2; 5 B[2, 0] = 1; B[2, 1] = 1; B[2, 2] = 1; B[2, 3] = 0; B[2, 4] = 0; 6 B[3, 0] = 0; B[3, 1] = 1; B[3, 2] = 1; B[3, 3] = 0; B[3, 4] = 0; 7 B[4, 0] = 0; B[4, 1] = 0; B[4, 2] = 1; B[4, 3] = 0; B[4, 4] = 0; 8 9 double[,] B_inverse = ClassMatrixCalculate.Inverse(B, orderNum);
2)运行结果
5、总结
矩阵求逆的原理并不复杂,就是使用初等行变换,但是实现起来需要牵扯到多层循环,并且循环变量的起始值和结束值需要根据情况变化,
还要判断所处位置的值是否为0或者为负的情况,所以需要理清思路之后进行编写。