MyMathLib系列(行列式计算3)

到今天，行列式和线性方程组部分就完成了。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MyMathLib
{

    /// <summary>
    /// 行列式计算,本程序属于MyMathLib的一部分，欢迎使用，参考，提意见。
    /// 有时间用函数语言改写，做自己得MathLib,里面的算法经过验证，但没经过
    /// 严格测试，如需参考，请慎重.
    /// </summary>
    public static partial class LinearAlgebra
    {#region 线性方程组
        /// <summary>
        /// 根据拉普拉斯定理计算行列式值。
        /// </summary>
        /// <param name="Determinants">N阶行列式</param>
        /// <returns>计算结果</returns>
        public static decimal CalcDeterByLaplaceLaw(decimal[,] Determinants)
        {
            var n = Determinants.GetLength(0);
            //如果阶数小于3，则没必要采用拉普拉斯展开
            if (n <= 3)
            {
                return CalcDeterminantAij(Determinants, false);
            }
            var theRows = GetLaplaceRowsOdd(n);
            return CalcDeterByLaplaceLaw(Determinants, theRows);
        }
        /// <summary>
        /// 求解线性方程组,这里要求N
        /// </summary>
        /// <param name="CoefficientDeterminant">线性方程组系数行列式</param>
        /// <param name="ConstantTerms">常数项</param>
        /// <returns></returns>
        public static decimal[] LinearEquations(int UnknownElements,decimal[,] CoefficientDeterminant, decimal[] ConstantTerms)
        {
            var theRowCount = CoefficientDeterminant.GetLength(0);
            var theColCount = CoefficientDeterminant.GetLength(1);
            if (UnknownElements == theRowCount && theColCount == UnknownElements)
            {
                var theD = CalcDeterByLaplaceLaw(CoefficientDeterminant);
                if(theD==0)
                {
                    return null;
                }
                decimal[] theResults = new decimal[UnknownElements];

                for (int i = 1; i <= UnknownElements; i++)
                {
                    //置换第i列，注意保存原来的值，下次计算时恢复
                    var theTemp = new decimal[UnknownElements];
                    for (int j = 1; j <= UnknownElements; j++)
                    {
                        theTemp[j-1] = CoefficientDeterminant[j-1, i-1];
                        CoefficientDeterminant[j - 1, i - 1] = ConstantTerms[j - 1];
                    }
                    var theDi = CalcDeterByLaplaceLaw(CoefficientDeterminant) / theD;
                    theResults[i-1] = theDi;
                    //复原系数行列式.
                    for (int j = 1; j <= UnknownElements; j++)
                    {
                        CoefficientDeterminant[j - 1, i - 1] = theTemp[j - 1];
                    }
                }
                return theResults;

            }
            else
            {
                throw new Exception("参数格式不正确!");
            }
        }
        /// <summary>
        /// 求解线性方程组(消元法),这里与化三角求行列式的方法类似，这里要求方程个数和元个数相同。
        /// 如果方程个数小于元的个数，消元没问题，但涉及到一般解，会牵扯到符号运算，这里暂不考虑.
        /// </summary>
        /// <param name="CoefficientDeterminant">线性方程组系数行列式,最右边N+1列是常数项</param>
        /// <returns></returns>
        public static decimal[] LinearEquationsEM(int UnknownElements,decimal[,] CoefficientDeterminant)
        {
            var theRowCount = CoefficientDeterminant.GetLength(0);
            var theColCount = CoefficientDeterminant.GetLength(1);
            if (UnknownElements == theRowCount && theColCount == UnknownElements + 1)
            {
                decimal[] theResults = new decimal[UnknownElements];
                int theN = UnknownElements;
                //从第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 (CoefficientDeterminant[j, i] == 0)
                        {
                            continue;
                        }

                        //如果[j,i]的上一行[j-1, i]的值为0则交换
                        if (CoefficientDeterminant[j - 1, i] == 0)
                        {
                            for (int k = 0; k <= theN; k++)//这里要交换常数项，所以是k <= theN
                            {
                                decimal 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 <= theN; k++)//这里要计算常数项，所以是k <= theN
                            {
                                CoefficientDeterminant[j, k] = CoefficientDeterminant[j, k] - CoefficientDeterminant[j - 1, k] * theRate;
                            }
                        }
                    }
                }
                //处理结果
                if (CoefficientDeterminant[UnknownElements - 1, UnknownElements - 1] == 0)
                {
                    if (CoefficientDeterminant[UnknownElements - 1, UnknownElements] == 0)
                    {
                        throw new Exception("无效方程,有无穷个解！");
                    }
                    else
                    {
                        throw new Exception("方程无解！");
                    }
                }
                //结果处理，回代
                for (int i = UnknownElements - 1; i >= 0; i--)
                {
                    //计算已求项
                    decimal theTempDec = 0;
                    for (int j = i + 1; j < theN; j++)
                    {
                        theTempDec += CoefficientDeterminant[i, j] * theResults[j];
                    }
                    //计算结果,如果系数为0，则为无效方程
                    if (CoefficientDeterminant[i, i] == 0)
                    {
                        throw new Exception("无效方程");
                    }
                    theResults[i] = (CoefficientDeterminant[i, UnknownElements] - theTempDec) / CoefficientDeterminant[i, i];
                }
                return theResults;

            }
            else
            {
                throw new Exception("参数格式不正确!");
            }
        }

        #endregion

    }
}