求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列时，不考虑第一个和第二个值；类推）；

>>>第二层循环，循环变量 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 lastRowN = n;
 43                 int colCount = 2 * n;
 44                 //把EB中索引为j的列的值化为1
 45                 for (int i = firstRowN; i < lastRowN; i++)
 46                 {
 47                     //如果EBijNum值为0，就把0所在的行与此刻最后一行调换位置
 48                     //并且循环变量i的终止值减去1,直到EBijNum值不为0
 49                     //最多调换到0所在的行的下一行
 50                     double EBijNum = EandB_normal[i, j];
 51                     while(EBijNum == 0 && lastRowN > i + 1)
 52                     {
 53                         for (int k = 0; k < colCount; k++)
 54                         {
 55                             rowHaveZero[k] = EandB_normal[i, k];
 56                             EandB_normal[i, k] = EandB_normal[lastRowN - 1, k];
 57                             EandB_normal[lastRowN - 1, k] = rowHaveZero[k];
 58                         }
 59                         lastRowN -= 1;
 60                         EBijNum = EandB_normal[i, j];
 61                     }
 62                     //如果while循环是由第二个判断跳出
 63                     //即EBijNum始终为0
 64                     if (EBijNum == 0)
 65                     {
 66                         //循环变量i的终止值再减去1，然后跳出循环
 67                         lastRowN -= 1;
 68                         break;
 69                     }
 70                     //如果为负数，该行变号
 71                     if(EBijNum < 0)
 72                     {
 73                         for(int k = 0; k < colCount; k++)
 74                         {
 75                             EandB_normal[i, k] = (-1) * EandB_normal[i, k];
 76                         }
 77                         EBijNum = EandB_normal[i, j];
 78                     }
 79                     //将该值变为1，则其余值都除以EBijNum
 80                     for (int k = 0; k < colCount; k++)
 81                     {
 82                         EandB_normal[i, k] = EandB_normal[i, k] / EBijNum;
 83                     }
 84                 }
 85 
 86                 //自n列起，每列只保留一个1，呈阶梯状
 87                 int secondRowN = firstRowN + 1;
 88                 for(int i = secondRowN; i < lastRowN; i++)
 89                 {
 90                     for(int k = 0; k < colCount; k++)
 91                     {
 92                         EandB_normal[i, k] = EandB_normal[i, k] 
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 - EandB_normal[firstRowN, k];
 93                     }
 94                 }
 95             
 96                 if(lastRowN == firstRowN)
 97                 {
 98                     //矩阵不满秩
 99                     IsFullRank = false;
100                     break;
101                 }
102             }
103             //不满秩，结束运算
104             if (!IsFullRank)
105             {
106                 double[,] error = new double[n, n];
107                 for(int i = 0; i < n; i++)
108                 {
109                     for(int j = 0; j < n; j++)
110                     {
111                         error[i, j] = 0;
112                     }
113                 }
114                 //返还值均为0的矩阵
115                 return error;
116                 
117             }
118            
119             //将上三角矩阵变为单位矩阵
120             for(int j = 2*n - 1; j > n; j--)
121             {
122                 //firstRowN为参考行
123                 //secondRowN为运算行
124                 int firstRowN = j - n;
125                 int secondRowN = firstRowN - 1;
126                 int colCount = j + 1;
127                 //从最后一列起，每列只保留一个1，其余减为0
128                 for(int i = secondRowN; i > -1; i--)
129                 {
130                     double EBijNum = EandB_normal[i, j];
131                     for(int k = 0; k < colCount; k++)
132                     {
133                         EandB_normal[i, k] = EandB_normal[i, k] 
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 - EandB_normal[firstRowN, k] * EBijNum;
134                     }
135                 }
136 
137             }
138 
139 
140 
141             //####提取逆矩阵####
142             for(int i = 0; i < n; i++)
143             {
144                 for(int j = 0; j < n; j++)
145                 {
146                     B_inverse[i, j] = EandB_normal[i, j];
147                 }
148             }
149 
150             return B_inverse;
151         }

 

 

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或者为负的情况，所以需要理清思路之后进行编写。

 

posted @ 2020-04-11 16:20  Ganders  阅读(1190)  评论(0)    收藏  举报