标签：
利用TExp类的运算来求矩阵的特征值，初等因子等：using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace MyMathLib
{
/// <summary>
/// 一元多项式计算
/// </summary>
public class PolynomialOfOneBasic
{
/// <summary>
/// 化成对阶梯矩阵
/// </summary>
/// <param name="CoefficientDeterminant">系数矩阵</param>
public static void TransToEchelonMatrix(TExp[,] CoefficientDeterminant)
{
var theRowCount = CoefficientDeterminant.GetLength(0);
var theColCount = CoefficientDeterminant.GetLength(1);
int theN = theRowCount;
int theE = theColCount;
//从第1列到第theE列.
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] == (TExp)0)
{
continue;
}
//******这里增加修改了判断是否继续交换消元的条件:
//如果左上邻元素[j-1, i-1]以及其左边的元素都为0方可交换
//因为当前元素的左边元素已经全部是零，因此如果要交换不能使本行左边产生非零数，
//则需要左上邻及其所有元素皆为0.
var theCanDo = true;
for (int s = i - 1; s >= 0; s--)
{
if (CoefficientDeterminant[j - 1, s] != (TExp)0)
{
theCanDo = false;
break;
}
}
if (theCanDo)
{
//如果[j,i]的上一行[j-1, i]的值为0则交换
if (CoefficientDeterminant[j - 1, i] == (TExp)0)
{
for (int k = 0; k < theE; k++)//这里要交换常数项
{
TExp theTmpDec = CoefficientDeterminant[j, k];
CoefficientDeterminant[j, k] = CoefficientDeterminant[j - 1, k];
CoefficientDeterminant[j - 1, k] = theTmpDec;
}
}
else
{
var theRate2 = CoefficientDeterminant[j, i];
var theRate1 = CoefficientDeterminant[j - 1, i];
for (int k = i; k < theE; k++)//这里要计算常数项
{
CoefficientDeterminant[j, k] =
CoefficientDeterminant[j, k] * theRate1
- CoefficientDeterminant[j - 1, k] * theRate2;
}
//将第i行，以最高次数系数来消除大数项.
double theMaxRation = 0;
uint theMaxPower = 0;
for (int k = i+1; k < theE; k++)//这里要计算常数项
{
var theMaxItem = CoefficientDeterminant[j, k].GetMaxPowerExp();
if (theMaxItem.Power > theMaxPower)
{
theMaxPower = theMaxItem.Power;
theMaxRation = theMaxItem.Ratio;
}
else if (theMaxItem.Power == theMaxPower)
{
if (Math.Abs(theMaxItem.Ratio) > Math.Abs(theMaxRation))
{
theMaxRation = theMaxItem.Ratio;
}
}
}
//
if (theMaxRation != 0)
{
theMaxRation = 1 / theMaxRation;
for (int k = i + 1; k < theE; k++)//这里要计算常数项
{
CoefficientDeterminant[j, k] = CoefficientDeterminant[j, k] * theMaxRation;
}
}
}
}
}
}
}
/// <summary>
/// 变换成对角矩阵,这个算法的问题在于会增加方程解。
/// </summary>
/// <returns>变换过程记录.</returns>
public static void TransToStandardForm(TExp[,] CoefficientDeterminant)
{
//先把矩阵转换成上三角矩阵。
TransToEchelonMatrix(CoefficientDeterminant);
//然后从最后一列开始，第1行开始变换为0.
var theColCount = CoefficientDeterminant.GetLength(0);
var theRowCount = CoefficientDeterminant.GetLength(1);
for (int j = theColCount - 1; j >= 0; j--)
{
//从下到上找到第1个非0行，作为基准行(减少行)
//因为矩阵的下半部分全为0,则开始找的位置在对角线上开始.
int theR = -1;
int theStartIndex = Math.Min(j, theRowCount - 1);
for (int i = theStartIndex; i >= 0; i--)
{
if (CoefficientDeterminant[i, j] != (TExp)0)
{
theR = i;
break;
}
}
//如果找到基准行，则开始变换，利用减去基准行*一个系数来消除第0到thR-1行的元素
if (theR >= 0)
{
for (int i = 0; i < theR; i++)
{
var theRate2 = CoefficientDeterminant[theR, j];
var theRate1 = CoefficientDeterminant[i, j];
for (int k = 0; k < theColCount; k++)//这里要计算常数项
{
CoefficientDeterminant[i, k] =
(CoefficientDeterminant[i, k] * theRate2)
- (CoefficientDeterminant[theR, k] * theRate1);
}
//将第i行，以最高次数系数来消除大数项.
double theMaxRation = 0;
uint theMaxPower = 0;
for (int k = 0; k < theColCount; k++)//这里要计算常数项
{
var theMaxItem = CoefficientDeterminant[i, k].GetMaxPowerExp();
if (theMaxItem.Power > theMaxPower)
{
theMaxPower = theMaxItem.Power;
theMaxRation = theMaxItem.Ratio;
}
else if (theMaxItem.Power == theMaxPower)
{
if (Math.Abs(theMaxItem.Ratio) > Math.Abs(theMaxRation))
{
theMaxRation = theMaxItem.Ratio;
}
}
}
//
if (theMaxRation != 0)
{
theMaxRation = 1 / theMaxRation;
for (int k = 0; k < theColCount; k++)//这里要计算常数项
{
CoefficientDeterminant[i, k] = CoefficientDeterminant[i, k] * theMaxRation;
}
}
}
}
}
}
/// <summary>
/// 转换到对角矩阵，不会增加矩阵乘积的次数，主要用于初等因子的寻找.
/// </summary>
/// <param name="Elements">Lambda矩阵.</param>
public static void TransToDiagMatrix(TExp[,] Elements)
{
var theRowCount = Elements.GetLength(0);
var theColCount = Elements.GetLength(1);
if (theRowCount != theColCount || theColCount<=0)
{
throw new Exception("本方法仅支持方阵,阶数n大于等于0！");
}
//主要用于Lambda矩阵的处理，这里从第theColCount开始处理
//相当于从右上角开始沿对角线进行处理[theRowCount-j-1,j]
for (int j = theRowCount-1; j >= 0; j--)
{
//寻找i行i列中的常数元素 0 表示没找到 1 表示遍历行找到，-1 表示列遍历找到.
//如果当前的元素为0，如果没找到，就把找到的非常数项交换到当前位置.
//非常数项的次数应该尽量的低.
int theNotZeroR = -1;
int theNotZeroC = -1;
uint theNotZeroPower = uint.MaxValue;
var theHasFind = 0;
var theIndex = -1;
for (int k = j; k >= 0; k--)
{
if (Elements[theRowCount - j - 1, k].GetMaxPowerExp().Power < theNotZeroPower)
{
if (!Elements[theRowCount - j - 1, k].IsZero)
{
theNotZeroPower = Elements[theRowCount - j - 1, k].GetMaxPowerExp().Power;
theNotZeroR = theRowCount - j - 1;
theNotZeroC = k;
}
}
if (Elements[theRowCount - k - 1, j].GetMaxPowerExp().Power < theNotZeroPower)
{
if (!Elements[theRowCount - k - 1, j].IsZero)
{
theNotZeroPower = Elements[theRowCount - k - 1, j].GetMaxPowerExp().Power;
theNotZeroR = theRowCount - k - 1;
theNotZeroC = j;
}
}
if (Elements[theRowCount - j - 1, k].IsNotZeroConst)
{
theHasFind = 1;
theIndex = k;
break;
}
if (Elements[theRowCount - k - 1, j].IsNotZeroConst)
{
theHasFind = -1;
theIndex = theRowCount - k - 1;
break;
}
}
if (theHasFind != 0)
{
switch (theHasFind)
{
case -1:
Elements.SwapRow(theIndex, j);
theNotZeroR = -1;
theNotZeroC = -1;
break;
case 1:
Elements.SwapCol(theIndex, j);
theNotZeroC =-1;
theNotZeroR = -1;
break;
}
//处理[theRowCount-j-1,j]所在行列的元素
RemoveNoZeroElement(Elements, theRowCount - j - 1, j);
}
else
{
//如果没找到常数项，且当前项为0项，则交换
if (Elements[theRowCount - j - 1, j].IsZero)
{
//行相同，交换列
if (theRowCount - j - 1 == theNotZeroR)
{
Elements.SwapCol(theNotZeroC, j);
}
else
{
if(j== theNotZeroC)
{
Elements.SwapRow(theNotZeroR, theRowCount - j - 1);
}
}
}
//继续找,这次是找两个相加后为常量的项
//以元素[theRowCount-j-1,j]为准，先向下找.
#region 在同一列上向下找
var theHasFind2 = false;
var theR1 = -1;
var theR2 = -1;
double theSign = 1;
//这里需要考虑到对应成比例的问题.
for (int r1 = theRowCount - j - 1; r1 < theRowCount; r1++)
{
if (Elements[r1, j].IsZero)
{
continue;
}
else
{
for (int r2 = r1 + 1; r2 < theRowCount; r2++)
{
//零元素可以不考虑.
if (Elements[r2, j].IsZero)
{
continue;
}
else
{
var theTempRatio1 = Elements[r1, j].GetMaxPowerItemRatio();
var theTempRatio2 = Elements[r2, j].GetMaxPowerItemRatio();
var theTempRate = theTempRatio1 / theTempRatio2;
var theTemp = Elements[r1, j] + (theTempRate * Elements[r2, j]);
if (theTemp.IsNotZeroConst || theTemp.IsZero)
{
theHasFind2 = true;
theR1 = r1;
theR2 = r2;
theSign = theTempRate;
break;
}
else
{
theTemp = Elements[r1, j] - (theTempRate * Elements[r2, j]);
if (theTemp.IsNotZeroConst || theTemp.IsZero)
{
theHasFind2 = true;
theR1 = r1;
theR2 = r2;
theSign = -theTempRate;
break;
}
}
}
}
if (theHasFind2)
{
break;
}
}
}
if (theHasFind2)
{
var theContinue2 = true;
if (theR1 == theRowCount - j - 1)
{
var theTempExp = Elements[theRowCount - j - 1, j] + Elements[theR2, j] * theSign;
if (theTempExp.IsZero)
{
theContinue2 = false;
}
}
if (theContinue2)
{
//1、先将[r,j]变为非零常量项
for (int c = 0; c <= j; c++)
{
Elements[theR1, c] = Elements[theR1, c] + (Elements[theR2, c] * theSign);
}
//2、交换theRowCount-j-1 行和r1(theR1)
Elements.SwapRow(theRowCount - j - 1, theR1);
}
//3、以[theRowCount-j-1,j]为准，变其所在行列的非零元素为零元素.
//处理[theRowCount-j-1,j]所在行列的元素
//如果当前元素为0，则继续处理当前行.
if (Elements[theRowCount - j - 1, j].IsZero)
{
j++;
}
else
{
RemoveNoZeroElement(Elements, theRowCount - j - 1, j);
}
}
#endregion
else
{
//以元素[theRowCount-j-1,j]为准，向左找.
#region 同一行上向左找.
var theHasFind3 = false;
var theC1 = -1;
var theC2 = -1;
double theSign2 = 1;
//这里需要考虑到对应成比例的问题.
for (int c1 = j; c1 >=0; c1--)
{
if (Elements[theRowCount-j-1, c1].IsZero)
{
continue;
}
else
{
for (int c2 = c1 -1; c2 >= 0; c2--)
{
//零元素可以不考虑.
if (Elements[theRowCount - j - 1, c2].IsZero)
{
continue;
}
else
{
var theTempRatio1 = Elements[theRowCount - j - 1, c1].GetMaxPowerItemRatio();
var theTempRatio2 = Elements[theRowCount - j - 1, c2].GetMaxPowerItemRatio();
var theTempRate = theTempRatio1 / theTempRatio2;
var theTemp = Elements[theRowCount - j - 1, c1] + (theTempRate * Elements[theRowCount - j - 1, c2]);
if (theTemp.IsNotZeroConst || theTemp.IsZero)
{
theHasFind3 = true;
theC1 = c1;
theC2 = c2;
theSign2 = theTempRate;
break;
}
else
{
theTemp = Elements[theRowCount - j - 1, c1] - (theTempRate * Elements[theRowCount - j - 1, c2]);
if (theTemp.IsNotZeroConst || theTemp.IsZero)
{
theHasFind3 = true;
theC1 = c1;
theC2 = c2;
theSign2 = 0.0-theTempRate;
break;
}
}
}
}
if (theHasFind3)
{
break;
}
}
}
if (theHasFind3)
{
//1、先将[theRowCount-j-1,theC1]变为非零常量项
//2、当如果变换使得当前元素为0，则不进行交换
var theContinue3 = true;
if (theC1 == j)
{
var theTempExp = Elements[theRowCount - j - 1, j] + Elements[theRowCount - j - 1, theC2] * theSign2;
if (theTempExp.IsZero)
{
theContinue3 = false;
}
}
if (theContinue3)
{
for (int r = theRowCount - j - 1; r < theRowCount; r++)
{
Elements[r, theC1] = Elements[r, theC1] + Elements[r, theC2] * theSign2;
}
//2、交换j列和c1(theC1)
Elements.SwapCol(j, theC1);
}
//3、以[theRowCount-j-1,j]为准，变其所在行列的非零元素为零元素.
//处理[theRowCount-j-1,j]所在行列的元素
//如果当前元素为0，则继续处理当前行.
if (Elements[theRowCount - j - 1, j].IsZero)
{
j++;
}
else
{
RemoveNoZeroElement(Elements, theRowCount - j - 1, j);
}
}
#endregion
}
}
}
}
/// <summary>
/// 消除Row行,Col列的非零元素
/// </summary>
/// <param name="Elements"></param>
/// <param name="Row"></param>
/// <param name="Col"></param>
private static void RemoveNoZeroElement(TExp[,] Elements,int Row, int Col)
{
var theRowCount = Elements.GetLength(0);
var theRate1 = Elements[Row, Col];
//var theConstRate = 1;// / ((TSymbolExp)theRate1).Ratio;
////处理列[theRowCount-j-1,0..j]
for (int c = 0; c < Col; c++)
{
//零元素不用消元
if (Elements[Row, c].IsZero)
{
continue;
}
var theRate2 = Elements[Row, c];
var theConstRate = theRate2.GetMaxPowerItemRatio() / theRate1.GetMaxPowerItemRatio();
//获取两个表达式的相关系数.
var theTempRatio = TExp.GetRelativRation(theRate1, theRate2);
//c列所在元素都需要变为0
double theAdjustRatio = 0;
uint theAdjustPower = 0;
for (int r = Row; r < theRowCount; r++)
{
if (theTempRatio.IsEqualTo(0))
{
Elements[r, c] = theConstRate * ((Elements[r, c] * theRate1) - (Elements[r, Col] * theRate2));
}
else
{
Elements[r, c] = (Elements[r, c] * theTempRatio) - (Elements[r, Col]);
}
var theTmpPower = Elements[r, c].GetMaxPowerExp().Power;
var theTmpRatio = Elements[r, c].GetMaxPowerItemRatio();
if (theTmpPower > theAdjustPower)
{
theAdjustPower = theTmpPower;
theAdjustRatio = theTmpRatio;
}
else if (theTmpPower == theAdjustPower)
{
if (Math.Abs(theAdjustRatio) < Math.Abs(theTmpRatio))
{
theAdjustRatio = theTmpRatio;
}
}
}
if (theAdjustPower > 0)
{
for (int r = Row; r < theRowCount; r++)
{
Elements[r, c] = Elements[r, c] * (1 / theAdjustRatio);
}
}
}
//处理行[theRowCount-j..theRowCount,j]
for (int r = Row + 1; r < theRowCount; r++)
{
//如果是零元素则继续.
if (Elements[r, Col].IsZero)
{
continue;
}
var theRate2 = Elements[r, Col];
var theConstRate = theRate2.GetMaxPowerItemRatio() / theRate1.GetMaxPowerItemRatio();
double theAdjustRatio = 0;
uint theAdjustPower = 0;
//获取两个表达式的相关系数.
var theTempRatio = TExp.GetRelativRation(theRate1, theRate2);
for (int c = 0; c <= Col; c++)
{
if (theTempRatio.IsEqualTo(0))
{
Elements[r, c] = theConstRate * ((Elements[r, c] * theRate1) - (Elements[Row, c] * theRate2));
}
else
{
Elements[r, c] = (Elements[r, c] * theTempRatio) - (Elements[Row, c]);
}
var theTmpPower = Elements[r, c].GetMaxPowerExp().Power;
var theTmpRatio = Elements[r, c].GetMaxPowerItemRatio();
if (theTmpPower > theAdjustPower)
{
theAdjustPower = theTmpPower;
theAdjustRatio = theTmpRatio;
}
else if (theTmpPower == theAdjustPower)
{
if (Math.Abs(theAdjustRatio) < Math.Abs(theTmpRatio))
{
theAdjustRatio = theTmpRatio;
}
}
}
if (theAdjustPower > 0)
{
for (int c = 0; c <= Col; c++)
{
Elements[r, c] = Elements[r, c] * (1 / theAdjustRatio);
}
}
}
}
}
}
用法：
var theBasis = new List<double[]>();
theBasis.Add(new double[] { 1, 0, 0 });
theBasis.Add(new double[] { 0, 1, 0 });
theBasis.Add(new double[] { 0, 0, 1 });
TExp x = new TSymbolExp() { Power = 1, Ratio = 1, Symbol = "x" };
var theP = new TExp[3, 3]{
{ 2-x, 3.0, 0.0},
{ 3, 2-x, 0.0},
{ 3, -3, 5-x}
};
MyMathLib.PolynomialOfOneBasic.TransToDiagMatrix(theP);
FormatViewArray(theP);这里的x为Lambda,就是λ-矩阵，theP等于(A-λE).计算结果会有小数上的差异，主要是浮点运算引起的，没有完全处理回去，大家注意就行，不影响结果。
线性代数部分就到此尾声了，后面会总结一下概念。
标签：
原文地址：http://blog.csdn.net/hawksoft/article/details/42501661
