在行列式消元中的判断条件有问题,这里给出订正后的代码:
1)MyMathLib.LinearAlgebra.CalcDeterminant方法
/// <summary>
/// 化三角法行列式计算,
/// </summary>
/// <param name="Determinants">N阶行列式</param>
/// <returns>计算结果</returns>
public static double CalcDeterminant(double[,] Determinants)
{
int theSign = 1;//正负符号,如果需要变换,将记录变换后的符号.
int theN = Determinants.GetLength(0);
//从第1列到第theN-1列
for (int i = 0; i < theN - 1; i++)
{
//从第theN-1行到第i+1行,将D[j,i]依次变为0
for (int j = theN - 1; j > i; j--)
{
//如果为当前值为0,则不处理,继续处理上一行
if (Determinants[j, i] == 0)
{
continue;
}
//修订部分:
//如果左上邻元素[j-1, i-1]以及其左边的元素都为0方可交换
//因为当前元素的左边元素已经全部是零,因此如果要交换不能使本行左边产生非零数,
//则需要左上邻及其所有元素皆为0.
var theCanDo = true;
for (int s = i - 1; s >= 0; s--)
{
if (Determinants[j - 1, s] != 0)
{
theCanDo = false;
break;
}
}
if (theCanDo)
{
//如果[j,i]的上一行[j-1, i]的值为0则交换
if (Determinants[j - 1, i] == 0)
{
//每次交换,行列式的值大小不变,符号取反
theSign = 0 - theSign;
for (int k = 0; k < theN; k++)
{
double theTmpDec = Determinants[j, k];
Determinants[j, k] = Determinants[j - 1, k];
Determinants[j - 1, k] = theTmpDec;
}
}
else
{
//将当前行减去上一行与theRate的积。
double theRate = Math.Round(Determinants[j, i] / Determinants[j - 1, i], ConstDef.Decimals);
for (int k = 0; k < theN; k++)
{
Determinants[j, k] = Math.Round(Determinants[j, k] - Determinants[j - 1, k] * theRate, ConstDef.Decimals);
}
}
}
}
}
//结果为对角线上的元素的乘积,注意符号的处理。
double theRetDec = theSign;
for (int i = 0; i < theN; i++)
{
theRetDec *= Determinants[i, i];
}
return theRetDec;
}
2.MyMathLib.LinearAlgebra.EquationsElimination
/// <summary>
/// 方程组消元,最后一列为系数,结果就在CoefficientDeterminant里.
/// 本算法也可以用来求矩阵的秩.
/// </summary>
/// <param name="CoefficientDeterminant">方程组系数数组</param>
public static void EquationsElimination(double[,] CoefficientDeterminant)
{
var theRowCount = CoefficientDeterminant.GetLength(0);
var theColCount = CoefficientDeterminant.GetLength(1);
int theN = theRowCount;
int theE = theColCount - 1;
//从第1列到第theE-1列,最后一列不用处理.
for (int i = 0; i < theE; i++)
{
//从第theN-1行到第1行,将D[j,i]依次变为0,需要注意的是:
//如果第j-1行,的左元素全部为0,才能继续交换.
for (int j = theN - 1; j > 0; j--)
{
//如果为当前值为0,则不处理,继续处理上一行
if (CoefficientDeterminant[j, i] == 0)
{
continue;
}
//******这里增加修改了判断是否继续交换消元的条件:
//如果左上邻元素[j-1, i-1]以及其左边的元素都为0方可交换
//因为当前元素的左边元素已经全部是零,因此如果要交换不能使本行左边产生非零数,
//则需要左上邻及其所有元素皆为0.
var theCanDo = true;
for (int s = i - 1; s >= 0; s--)
{
if (CoefficientDeterminant[j - 1, s] != 0)
{
theCanDo = false;
break;
}
}
if (theCanDo)
{
//如果[j,i]的上一行[j-1, i]的值为0则交换
if (CoefficientDeterminant[j - 1, i] == 0)
{
for (int k = 0; k <= theE; k++)//这里要交换常数项,所以是k <= theE
{
double theTmpDec = CoefficientDeterminant[j, k];
CoefficientDeterminant[j, k] = CoefficientDeterminant[j - 1, k];
CoefficientDeterminant[j - 1, k] = theTmpDec;
}
}
else
{
//将当前行减去上一行与theRate的积。
//var theRate = CoefficientDeterminant[j, i] / CoefficientDeterminant[j - 1, i];
//for (int k = 0; k <= theE; k++)//这里要计算常数项,所以是k <= theE
//{
// CoefficientDeterminant[j, k] = CoefficientDeterminant[j, k] - CoefficientDeterminant[j - 1, k] * theRate;
//}
//改进:做乘法,可以避免小数换算带来的误差
var theRate2 = CoefficientDeterminant[j, i];
var theRate1 = CoefficientDeterminant[j - 1, i];
for (int k = 0; k <= theE; k++)//这里要计算常数项,所以是k <= theE
{
CoefficientDeterminant[j, k] = CoefficientDeterminant[j, k] * theRate1 - CoefficientDeterminant[j - 1, k] * theRate2;
}
}
}
}
}
}