求n阶矩阵B的逆矩阵(C#)

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列时,不考虑第一个和第二个值;类推);

扫描二维码关注公众号,回复: 10695005 查看本文章

>>>第二层循环,循环变量 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或者为负的情况,所以需要理清思路之后进行编写。

猜你喜欢

转载自www.cnblogs.com/Ganders/p/12680620.html