MyMathLib系列(行列式計算4--向量部分)

MyMathLib系列(行列式計算4--向量部分)

1)將向量組進行消元，變換成階梯矩陣，這是求向量組的極大線性無關組的基本算法。這個方法在前面曾經給出過，但在這裡做了改進，目的是為了可以判斷是否線性相關：

        /// 

        /// 方程組消元，最後一列為系數,結果就在CoefficientDeterminant裡.
        /// 本算法也可以用來求矩陣的秩.
        /// 
        /// 方程組系數數組
        public static void EquationsElimination(decimal[,] 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.
                    for (int s = i - 1; s >= 0; s--)
                    {
                        if (CoefficientDeterminant[j - 1, s] != 0)
                        {
                            break;
                        }
                    }
                    //如果[j,i]的上一行[j-1, i]的值為0則交換
                    if (CoefficientDeterminant[j - 1, i] == 0)
                    {
                        for (int k = 0; k <= theE; k++)//這裡要交換常數項，所以是k <= theE
                        {
                            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 <= 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;
                        }
                    }
                }
            }
        }

後面的矩陣求秩的過程也會用到這種消元法，不過後面的矩陣求秩在這個算法上略有修改，采用了另外一種方式。下面是向量相關的算法：

/*Vector.cs
 * Albert.Tian on 20141225 
 */
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MyMathLib
{
    /// 

    /// 向量相關運算，這裡沒有給出向量的基本變換。因為
    /// 向量組可以用矩陣表達，因此很多向量的運算最終都歸於矩陣運算。
    /// 本算法中的求極大線性無關組，是根據課本算法得出，不是最簡單的一種方法，
    /// 在矩陣表達中，可以根據矩陣的初等變換，很容易求到極大線性無關組.
    /// 
    public class Vectors
    {
        /// 

        /// 判斷從0到指定結束位的元素是否全是0.
        /// 
        /// 當前向量
        /// 結束索引
        /// 
        public static bool IsZeroVector(decimal[] CurrVector, int EndIndex = 0)
        {
            int theEndIndex = EndIndex;
            if (theEndIndex == 0)
            {
                theEndIndex = CurrVector.Length;
            }

            var theIsZeroVector = true;
            for (int i = 0; i < theEndIndex; i++)
            {
                if (CurrVector[i] != 0)
                {
                    theIsZeroVector = false;
                    break;
                }
            }
            return theIsZeroVector;
        }
        /// 

        /// 判斷當前向量是否與前面的向量線性相關，采用的是消元化
        /// 
        /// 當前向量
        /// 已有的線性無關的向量組
        /// 
        private static bool IsLinearCorrelation(decimal[] CurrVector, List PrevMaxLIVectors)
        {
            //零向量必是線性相關
            var theIsZeroVector = IsZeroVector(CurrVector);
            if (theIsZeroVector)
            {
                return true;
            }
            //如果前面沒有向量，則當前非零向量必是線性無關的.
            if (PrevMaxLIVectors.Count <= 0)
            {
                return false;
            }
            //構造方程式，判斷是否線性相關
            var theECount = CurrVector.Length;
            var theVCount = PrevMaxLIVectors.Count;
            
            //加入當前向量為常量向量,因此方程組的列為theVCount+1.
            var theEquealGroup = new decimal[theECount, theVCount + 1];
            //導入PrevMaxLIVectors
            for (int i = 0; i < theVCount; i++)
            {
                for (int j = 0; j < theECount; j++)
                {
                    theEquealGroup[j, i] = PrevMaxLIVectors[i][j];
                }
            }

            //加入當前向量作為常量向量
            for (int j = 0; j < theECount; j++)
            {
                theEquealGroup[j, theVCount] = CurrVector[j];
            }
            //消元,這裡的消元法變成了求矩陣秩的基本算法.
            LinearAlgebra.EquationsElimination(theEquealGroup);
            //如果theRA=theRA1，則至少有一個非零解。
            //這裡如此判斷，主要是因為消元的方法。當然，這在矩陣有證。
            var theRA = 0;
            var theRA1 = 0;
            for (int i = 0; i < theECount; i++)
            {
                if (!IsZeroVector(theEquealGroup.GetVector(i), theVCount))
                {
                    theRA++;
                }
                if (!IsZeroVector(theEquealGroup.GetVector(i)))
                {
                    theRA1++;
                }
            }
            return (theRA >= theRA1);
        }
        /// 

        /// 這個函數暫時無用，判斷兩個向量的等價性
        /// 
        /// 
        /// 
        /// 
        private static bool IsEquivalent(decimal[] V1, decimal[] V2)
        {
            var theV1IsZeroVector = IsZeroVector(V1);
            var theV2IsZeroVector = IsZeroVector(V2);
            //兩個都是0向量
            if (theV1IsZeroVector && theV2IsZeroVector)
            {
                return true;
            }
            //其中一個零向量，另一個是非零向量
            if (theV1IsZeroVector || theV2IsZeroVector)
            {
                return false;
            }
            //一個分量肯定成比例
            if (V1.Length <= 1)
            {
                return true;
            }

            decimal thePreRate = 0;
            var theCount = V1.Length;
            var theIndex = 0;
            for (int i = 0; i < theCount; i++)
            {
                //如果都等於0，則不比較
                if (V1[i] == V2[i] && V2[i] == 0)
                {
                    continue;
                }
                //如果其中一個是0，則必然不成比例,不等價
                if (V1[i] == 0 || V2[i] == 0)
                {
                    return false;
                }
                if (theIndex == 0)
                {
                    thePreRate = V1[i] / V2[i];
                    theIndex++;
                }
                else
                {

                    if (thePreRate != V1[i] / V2[i])
                    {
                        return false;
                    }

                }
            }
            return true;
        }
        /// 

        /// 獲取向量組的極大線性無關組，方法是依照教科書寫的。
        /// 
        /// 向量組
        /// 
        public static List GetMaxLinearIndependentGroup(decimal[,] Vectors)
        {
            //向量個數
            var theVCount = Vectors.GetLength(0);
            //元數個數
            var theECount = Vectors.GetLength(1);
            List theTempVectors = new List();
            for (int i = 0; i < theVCount; i++)
            {
                var theTempVector = Vectors.GetVector(i);
                //如果當前向量不能被前面的向量線性表出，則加入到極大組
                //能不能線性表出，實際上就是求當前向量加入前面的向量組後，
                //還是不是線性相關
                if (!IsLinearCorrelation(theTempVector, theTempVectors))
                {
                    theTempVectors.Add(theTempVector);
                }
            }
            return theTempVectors;
        }
        
    }
}