降階法計算行列式方法有個地方有Bug(原文也已更正,此為更正後部分)
今天用此函數做方程求解時發現有誤,特此更正:
///
/// 降階法計算行列式
///
/// N階行列式
/// 是否0優化
/// 計算結果
public static decimal CalcDeterminantAij(decimal[,] Determinants, bool ZeroOptimization = false)
{
var theN = Determinants.GetLength(0);
//如果為2階,直接計算
if (theN == 2)
{
return Determinants[0, 0] * Determinants[1, 1] - Determinants[0, 1] * Determinants[1, 0];
}
if (theN == 1)
{
return Determinants[0, 0];
}
if (theN == 0)
{
throw new Exception("參數錯誤!");
}
if (ZeroOptimization)
{
//找0最多的行
int theRowIndex = 0;
int theMaxZeroCountR = -1;
for (int i = 0; i < theN; i++)
{
int theZeroNum = 0;
for (int j = 0; j < theN; j++)
{
if (Determinants[i, j] == 0)
{
theZeroNum++;
}
}
if (theZeroNum > theMaxZeroCountR)
{
theRowIndex = i;
theMaxZeroCountR = theZeroNum;
}
}
//找0最多的列
int theColIndex = 0;
int theMaxZeroCountC = -1;
for (int i = 0; i < theN; i++)
{
int theZeroNum = 0;
for (int j = 0; j < theN; j++)
{
if (Determinants[j, i] == 0)
{
theZeroNum++;
}
}
if (theZeroNum > theMaxZeroCountC)
{
theColIndex = i;
theMaxZeroCountC = theZeroNum;
}
}
if (theMaxZeroCountR >= theMaxZeroCountC)
{
decimal theRetDec = 0;
//第i=theRowIndex+1行展開
int i = theRowIndex + 1;
for (int j = 1; j <= theN; j++)
{
var theSign = CalcDeterMijSign(i, j);
var theNewMij = GetDeterminantMij(Determinants, i, j);
theRetDec += theSign * Determinants[i - 1, j - 1] * CalcDeterminantAij(theNewMij, ZeroOptimization);
}
return theRetDec;
}
else
{
decimal theRetDec = 0;
//第j=theColIndex+1列展開
int j = theColIndex + 1;
for (int i = 1; i <= theN; i++)
{
var theSign = CalcDeterMijSign(i, j);
var theNewMij = GetDeterminantMij(Determinants, i, j);
theRetDec += theSign * Determinants[i, j] * CalcDeterminantAij(theNewMij, ZeroOptimization);
}
return theRetDec;
}
}
else
{
//采用隨機法展開一行
var i = new Random().Next(1, theN);
decimal theRetDec = 0;
for (int j = 1; j <= theN; j++)
{
var theSign = CalcDeterMijSign(i, j);
var theNewMij = GetDeterminantMij(Determinants, i, j);
//此處修改theRetDec += theSign * Determinants[i, j] * CalcDeterminantAij(theNewMij, ZeroOptimization);
theRetDec += theSign * Determinants[i-1, j-1] * CalcDeterminantAij(theNewMij, ZeroOptimization);
}
return theRetDec;
}
}