C#，码海拾贝（14）——实矩阵与复数矩阵“求逆”的高斯·约当算法，《C#数值计算算法编程》源代码升级改进版
2023-09-11 14:15:48 时间
1、矩阵求逆
矩阵求逆，即求矩阵的逆矩阵。矩阵是线性代数的主要内容，很多实际问题用矩阵的思想去解既简单又快捷。逆矩阵又是矩阵理论的很重要的内容，逆矩阵的求法自然也就成为线性代数研究的主要内容之一。
设A是数域上的一个n阶方阵，若在相同数域上存在另一个n阶矩B，使得： AB=BA=E。 则我们称B是A的逆矩阵，而A则被称为可逆矩阵。其中，E为单位矩阵。
典型的矩阵求逆方法有：利用定义求逆矩阵、初等变换法、伴随阵法、恒等变形法等。
2、矩阵求逆的计算过程简述
求元素为具体数字的矩阵的逆矩阵，常用初等变换法‘如果A可逆，则A’可通过初等变换，化为单位矩阵 I ，即存在初等矩阵使 P1,P2,P3,...Ps：
（1）P1 x P2 x P3 x ...x Ps x A = I；
（2）用 A-1 右乘上式两端，得：P1 x P2 x P3 x ...x Ps x I = A-1；
比较（1）、（2）两式，可以看到当A通过初等变换化为单位处阵的同时，对单位矩阵I作同样的初等变换，就化为A的逆矩阵 。
这就是求逆矩阵的初等行变换法，它是实际应用中比较简单的一种方法。需要注意的是，在作初等变换时只允许作行初等变换。同样，只用列初等变换也可以求逆矩阵。
用此方法求逆矩阵，对于小型矩阵，特别是二阶方阵求逆既方便、快捷，又有规律可循。因为二阶可逆矩阵的伴随矩阵，只需要将主对角线元素的位置互换，次对角线的元素变号即可。
若可逆矩阵是二阶或二阶以上矩阵，在求逆矩阵的过程中，需要求9个或9个以上代数余子式，还要计算一个三阶或三阶以上行列式，工作量大且中途难免出现符号及计算的差错。对于求出的逆矩阵是否正确，一般要通过 来检验。一旦发现错误，必须对每一计算逐一排查。
using System;
namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 矩阵类
/// 作者：周长发
/// 改进：深度混淆
/// https://blog.csdn.net/beijinghorn
/// </summary>
public partial class Matrix
{
/// <summary>
/// 实矩阵求逆的全选主元高斯－约当法
/// </summary>
/// <param name="src">源矩阵</param>
/// <returns>求逆是否成功</returns>
public static bool InvertGaussJordan(Matrix src)
{
// 分配内存
int[] pnRow = new int[src.Columns];
int[] pnCol = new int[src.Columns];
// 消元
for (int k = 0; k <= src.Columns - 1; k++)
{
double d = 0.0;
for (int i = k; i <= src.Columns - 1; i++)
{
for (int j = k; j <= src.Columns - 1; j++)
{
int w = i * src.Columns + j;
double p = Math.Abs(src[w]);
if (p > d)
{
d = p;
pnRow[k] = i;
pnCol[k] = j;
}
}
}
// 失败
if (Math.Abs(d) < float.Epsilon)
{
return false;
}
if (pnRow[k] != k)
{
for (int j = 0; j <= src.Columns - 1; j++)
{
int u = k * src.Columns + j;
int v = pnRow[k] * src.Columns + j;
double p = src[u];
src[u] = src[v];
src[v] = p;
}
}
if (pnCol[k] != k)
{
for (int i = 0; i <= src.Columns - 1; i++)
{
int u = i * src.Columns + k;
int v = i * src.Columns + pnCol[k];
double p = src[u];
src[u] = src[v];
src[v] = p;
}
}
{
int w = k * src.Columns + k;
src[w] = 1.0 / src[w];
for (int j = 0; j <= src.Columns - 1; j++)
{
if (j != k)
{
int u = k * src.Columns + j;
src[u] = src[u] * src[w];
}
}
for (int i = 0; i <= src.Columns - 1; i++)
{
if (i != k)
{
for (int j = 0; j <= src.Columns - 1; j++)
{
if (j != k)
{
int u = i * src.Columns + j;
src[u] = src[u] - src[i * src.Columns + k] * src[k * src.Columns + j];
}
}
}
}
for (int i = 0; i <= src.Columns - 1; i++)
{
if (i != k)
{
int u = i * src.Columns + k;
src[u] = -src[u] * src[w];
}
}
}
}
// 调整恢复行列次序
for (int k = src.Columns - 1; k >= 0; k--)
{
if (pnCol[k] != k)
{
for (int j = 0; j <= src.Columns - 1; j++)
{
int u = k * src.Columns + j;
int v = pnCol[k] * src.Columns + j;
double p = src[u];
src[u] = src[v];
src[v] = p;
}
}
if (pnRow[k] != k)
{
for (int i = 0; i <= src.Columns - 1; i++)
{
int u = i * src.Columns + k;
int v = i * src.Columns + pnRow[k];
double p = src[u];
src[u] = src[v];
src[v] = p;
}
}
}
// 成功返回
return true;
}
/// <summary>
/// 复矩阵求逆的全选主元高斯－约当法
/// </summary>
/// <param name="src">源矩阵</param>
/// <param name="mtxImag">复矩阵的虚部矩阵，当前矩阵为复矩阵的实部</param>
/// <returns>求逆是否成功</returns>
public static bool InvertGaussJordan(Matrix src, Matrix mtxImag)
{
int i, j, k, z, u, v, w;
double p, q, s, t, d, b;
// 分配内存
int[] pnRow = new int[src.Columns];
int[] pnCol = new int[src.Columns];
// 消元
for (k = 0; k <= src.Columns - 1; k++)
{
d = 0.0;
for (i = k; i <= src.Columns - 1; i++)
{
for (j = k; j <= src.Columns - 1; j++)
{
u = i * src.Columns + j;
p = src[u] * src[u] + mtxImag[u] * mtxImag[u];
if (p > d)
{
d = p;
pnRow[k] = i;
pnCol[k] = j;
}
}
}
// 失败
if (Math.Abs(d) < float.Epsilon)
{
return false;
}
if (pnRow[k] != k)
{
for (j = 0; j <= src.Columns - 1; j++)
{
u = k * src.Columns + j;
v = pnRow[k] * src.Columns + j;
t = src[u];
src[u] = src[v];
src[v] = t;
t = mtxImag[u];
mtxImag[u] = mtxImag[v];
mtxImag[v] = t;
}
}
if (pnCol[k] != k)
{
for (i = 0; i <= src.Columns - 1; i++)
{
u = i * src.Columns + k;
v = i * src.Columns + pnCol[k];
t = src[u];
src[u] = src[v];
src[v] = t;
t = mtxImag[u];
mtxImag[u] = mtxImag[v];
mtxImag[v] = t;
}
}
z = k * src.Columns + k;
src[z] = src[z] / d;
mtxImag[z] = -mtxImag[z] / d;
for (j = 0; j <= src.Columns - 1; j++)
{
if (j != k)
{
u = k * src.Columns + j;
p = src[u] * src[z];
q = mtxImag[u] * mtxImag[z];
s = (src[u] + mtxImag[u]) * (src[z] + mtxImag[z]);
src[u] = p - q;
mtxImag[u] = s - p - q;
}
}
for (i = 0; i <= src.Columns - 1; i++)
{
if (i != k)
{
v = i * src.Columns + k;
for (j = 0; j <= src.Columns - 1; j++)
{
if (j != k)
{
u = k * src.Columns + j;
w = i * src.Columns + j;
p = src[u] * src[v];
q = mtxImag[u] * mtxImag[v];
s = (src[u] + mtxImag[u]) * (src[v] + mtxImag[v]);
t = p - q;
b = s - p - q;
src[w] = src[w] - t;
mtxImag[w] = mtxImag[w] - b;
}
}
}
}
for (i = 0; i <= src.Columns - 1; i++)
{
if (i != k)
{
u = i * src.Columns + k;
p = src[u] * src[z];
q = mtxImag[u] * mtxImag[z];
s = (src[u] + mtxImag[u]) * (src[z] + mtxImag[z]);
src[u] = q - p;
mtxImag[u] = p + q - s;
}
}
}
// 调整恢复行列次序
for (k = src.Columns - 1; k >= 0; k--)
{
if (pnCol[k] != k)
{
for (j = 0; j <= src.Columns - 1; j++)
{
u = k * src.Columns + j;
v = pnCol[k] * src.Columns + j;
t = src[u];
src[u] = src[v];
src[v] = t;
t = mtxImag[u];
mtxImag[u] = mtxImag[v];
mtxImag[v] = t;
}
}
if (pnRow[k] != k)
{
for (i = 0; i <= src.Columns - 1; i++)
{
u = i * src.Columns + k;
v = i * src.Columns + pnRow[k];
t = src[u];
src[u] = src[v];
src[v] = t;
t = mtxImag[u];
mtxImag[u] = mtxImag[v];
mtxImag[v] = t;
}
}
}
// 成功返回
return true;
}
}
}
相关文章
- Word控件Spire.Doc 【超链接】教程(8)：在 C#/VB.NET 中链接到 Word 文档中的书签
- C#【必备技能篇】委托Action、 Action＜T＞、Func＜T＞、Predicate＜T＞（精简易懂，赞！）
- C#，图像二值化（10）——全局阈值的灰度平均值算法（Gray-Average Thresholding）及其源代码
- C#，核心基础算法——大数计算类Skyiv.BigInterger和任意精度算术运算类BigArithmetic源代码
- C#，斯特林数（Stirling Number）的算法与源代码
- C#，动态规划的集合划分问题（DP Partition problem）算法与源代码
- C#，双向链表（Doubly Linked List）归并排序（Merge Sort）算法与源代码
- C#，计算几何，计算机图形学（Computer Graphics）洪水填充算法（Flood Fill Algorithm）与源代码
- C#，笛卡尔树（Cartesian Tree）的构造、遍历算法与源代码
- C#，计算几何，随机点集之三角剖分的德劳内（Delaunay）算法的源代码
- C#，最小生成树（MST）博鲁夫卡（Boruvka）算法的源代码
- C#，数值计算，求平方根之巴比伦算法（Babylonian algorithm）的源代码
- C#，简单选择排序算法（Simple Select Sort）的源代码与数据可视化
- C#，码海拾贝（19）——一般实矩阵的QR分解（QR Decomposition）方法之C#源代码，《C#数值计算算法编程》源代码升级改进版
- C#，码海拾贝（15）——“对称正定矩阵”的求逆和“托伯利兹矩阵”求逆的“埃兰特”方法之C#源代码，《C#数值计算算法编程》源代码升级改进版
- C#，码海拾贝（11）——拉格朗日（Lagrange）三点式曲面插值算法，《C#数值计算算法编程》源代码升级改进版
- C#，码海拾贝（09）——阿克玛（Akima）曲线光滑插值算法，《C#数值计算算法编程》源代码升级改进版
- C#，码海拾贝（06）——连分式（Continued Fraction）曲线插值算法，《C#数值计算算法编程》源代码升级改进版
- C#，码海拾贝（05）——拉格朗日（Lagrange）三点式曲线插值算法，《C#数值计算算法编程》源代码升级改进版
- C#，码海拾贝（04）——拉格朗日（Lagrange）曲线插值算法及其拓展，《C#数值计算算法编程》源代码升级改进版
- C#，数值计算（Numerical Recipes in C#），矩阵的奇异值分解（SVD，Singular Value Decomposition）算法与源代码
- C#，初学琼林（06）——组合数的算法、数据溢出问题的解决方法及相关C#源代码
- C#，码海拾贝（16）——求“矩阵秩”的全选主元“高斯消去法（Gauss Elimination）”C#源代码，《C#数值计算算法编程》源代码升级改进版
- C# 贪心算法 -- 钱币找零问题
- C# checked和unchecked运算符
- C# 数据处理 相关算法