矩阵计算

为了找到产品运动轨迹的极值，面对一堆离散坐标，便很自然想到先模拟一条曲线图，再根据得到的曲线图，进而求取该曲线的极值。

想法虽然很清晰，但由于多年没有用到高等数学，竟然一时茫然：）。

片刻之后，开始整理思路：

如何模拟一条曲线？从网络中搜索模拟过程，在大量信息中，出现最多的词就是最小二乘法；于是上网找一下最小二乘法的原理，凭着自己的记忆，终于明白了它的意思，于是便开始了考试一样的演算。但得到的结果是，我的演算纸竟然不够长，单项长度太长，在失望中，我又一次茫然。

我回到了起点，看到最初的函数组，在网上寻找函数组的解决方案，依然没有太大收获；无聊之际，我拿起线性代数书浏览，正当我漫无目的的浏览时，一个似曾熟悉的词——矩阵突然出现，在潜意识里，我知道我有新的发现，于是稍微认真的看完了矩阵一章。最后我明白了一个新的解方程组的方法——通过矩阵计算，求解方程组（此时知道，几年的大学是没有白上的）。

根据网上资料，和自己的调整，写出了一个矩阵计算方法：

1）求矩阵的逆矩阵

/// <summary>
/// 求逆矩阵

/// 算法的大致思想是通过行列式初等变换来求。

/// </summary>
/// <param name="dMatrix"></param>
/// <returns></returns>
public static double[,] ReverseMatrix(double[,] dMatrix)
{
//获取矩阵的行数
int Level = dMatrix.GetLength(1);

////求矩阵的行列式的值
//double dMatrixValue = MatrixValue(dMatrix);

//if (dMatrixValue == 0)
//{
// return null;
//}

double[,] dReverseMatrix = new double[Level,2*Level];

//初始化矩阵Level×（2*Level）
for (int i = 0; i < Level; i++)
{
for (int j = 0; j < 2*Level; j++)
{
if (j < Level)
{
dReverseMatrix[i, j] = dMatrix[i, j];
}
else
{
if (j - Level == i)
{
dReverseMatrix[i, j] = 1;
}
else
{
dReverseMatrix[i, j] = 0;
}
}
}
}

for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
{
if (dReverseMatrix[i, j] == 0)
{
if (i == Level - 1)
{
return null;
}
int m = i + 1;

for (; dMatrix[m, j] == 0; m++)
{
if (m == Level - 1)
{
return null;
}
}

if (m == Level)
{
return null;
}
else
{
//把i行和m行相加
for (int n = j; n < 2*Level; n++)
{
dReverseMatrix[i, n] += dReverseMatrix[m, n];
}
}
}

double temp = dReverseMatrix[i, j];

if (temp != 1)
{
//把i行数据，变成以1开始的一行数据
for (int n = j; n < 2*Level; n++)
{
if (dReverseMatrix[i, n] != 0)
{
dReverseMatrix[i, n] /= temp;
}
}
}

//把i行后的所有行的j列变成0
for (int s = Level - 1; s > i; s--)
{
temp = dReverseMatrix[s, j];

for (int t = j; t < 2*Level; t++)
{
dReverseMatrix[s, t] -= (dReverseMatrix[i, t]*temp);
}
}
}

//把矩阵Level×（2*Level）前Level×Level转变为单位矩阵
for (int i = Level - 2; i >= 0; i--)
{
for (int j = i + 1; j < Level; j++)
{
if (dReverseMatrix[i, j] != 0)
{
double tmp = dReverseMatrix[i, j];

for (int n = j; n < 2*Level; n++)
{
dReverseMatrix[i, n] -= (tmp*dReverseMatrix[j, n]);
}
}
}
}

//返回逆矩阵
double[,] dReturn = new double[Level,Level];

for (int i = 0; i < Level; i++)
{
for (int j = 0; j < Level; j++)
{
dReturn[i, j] = dReverseMatrix[i, j + Level];
}
}
return dReturn;
}

2）求矩阵对应行列式的值

/// <summary>
/// 求矩阵对应行列式的值
/// </summary>
/// <param name="MatrixList"></param>
/// <returns></returns>
public static double MatrixValue(double[,] MatrixList)
{
int Level = MatrixList.GetLength(1);

double[,] dMatrix = new double[Level,Level];

for (int i = 0; i < Level; i++)
{
for (int j = 0; j < Level; j++)
{
dMatrix[i, j] = MatrixList[i, j];
}
}

int sign = 1;

for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
{
//判断改行dMatrix[i, j]是否为0，若是，则寻找i后的行（m,m>i，切dMatrix[m, j]!=0)进行交换
if (dMatrix[i, j] == 0)
{
if (i == Level - 1)
{
return 0;
}

int m = i + 1;

//获取一个dMatrix[m, j]不为为0的行
for (; dMatrix[m, j] == 0; m++)
{
if (m == Level - 1)
{
return 0;
}
}

//判断是否达到矩阵的最大行，若是，则返回0

//把i行和m行调换
double temp;
for (int n = j; n < Level; n++)
{
temp = dMatrix[i, n];
dMatrix[i, n] = dMatrix[m, n];
dMatrix[m, n] = temp;
}
sign *= (-1);
}

//把当前行以后的行所对应的列变成0
double tmp;
for (int s = Level - 1; s > i; s--)
{
tmp = dMatrix[s, j];
//j行后面的所有行
for (int t = j; t < Level; t++)
{
dMatrix[s, t] -= dMatrix[i, t]*(tmp/dMatrix[i, j]);
}
}
}

double result = 1;
for (int i = 0; i < Level; i++)
{
if (dMatrix[i, i] != 0)
{
result *= dMatrix[i, i];
}
else
{
return 0;
}
}

return sign*result;
}

3）求两个矩阵相乘的结果

/// <summary>
/// 矩阵乘法实现，即获取两个矩阵相乘的结果
/// </summary>
/// <param name="firstMatrix"></param>
/// <param name="secondMatrix"></param>
/// <returns></returns>
public static double[,] MatrixMultiplication(double[,] firstMatrix, double[,] secondMatrix)
{
double[,] resultMatrix = new double[firstMatrix.GetLength(0),secondMatrix.GetLength(1)];

//判断相乘矩阵是否合法，即第一个矩阵的列要等于第二个矩阵的行
if (firstMatrix.GetLength(1) != secondMatrix.GetLength(0))
{
return null;
}

//求结果矩阵
for (int rowIndex = 0; rowIndex < firstMatrix.GetLength(0); rowIndex++)
{
for (int colIndex = 0; colIndex < secondMatrix.GetLength(1); colIndex++)
{
//初始化结果矩阵的元素
resultMatrix[rowIndex, colIndex] = 0;

for (int i = 0; i < firstMatrix.GetLength(1); i++)
{
//求结果矩阵的元素值
resultMatrix[rowIndex, colIndex] += firstMatrix[rowIndex, i] * secondMatrix[i, colIndex];
}
}
}

return resultMatrix;
}