在行列式消元中的判断条件有问题,这里给出订正后的代码:
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; } } } } } }
原文:http://blog.csdn.net/hawksoft/article/details/42222403