Crout分解法

前面介绍了Cholesky分解法，但是Cholesky分解法只是用于对称正定矩阵，对于不是对称矩阵、或不是正定矩阵时，要计算线性方程组时，在矩阵阶数不是很大时可以采用Cramer法则，也可以采用高斯消元法来求解！顺便介绍一下高斯消元法，对于一个n阶矩阵，用高斯消元法不计加减运算，只记乘除运算，要运算n的立方/3+n的平方-n/3次，计算量比较大，而且精度较差。所以后来有了改进的高斯消元法——高斯主元素消去法，即先找出主元，然后利用高斯消元法来计算，也是实际计算中常用的一直方法。

但是，如果先将矩阵进行分解成一个下三角矩阵与一个上三角矩阵的程序，然后计算线性方程组的工作量将会大大减少，这就引出了矩阵的LU分解。Crout分解就是其中一种常用的方法，其思路比较适合计算机程序设计。

Crout分解法

Crout分解法可用于非对称、非正定的矩阵。当然，如果矩阵是对称正定矩阵，Crout分解法当然也适用。

2）程序设计

程序设计主要分为计算LU矩阵、计算Y矩阵和计算X矩阵三个部分。

（1）计算LU矩阵

计算LU矩阵的程序主要根据式（5-6）和（5-7）来设计，其代码如下：

//计算LU矩阵

double[,] L = new double[arrayCount, arrayCount];

double[,] U = new double[arrayCount, arrayCount];

for (i = 0; i < arrayCount; i++)

{

U[i, i] = 1;

}

for (k = 0; k < arrayCount; k++)

{

for (i = k; i < arrayCount; i++)

{

temp = 0.0;

for (r = 0; r < k ; r++)

{

temp+=L[i,r]*U[r,k];

}

L[i, k] = A[i, k] - temp;

}

for (j = k+1; j < arrayCount; j++)

{

temp = 0.0;

for (r = 0; r < k ; r++)

{

temp += L[k, r] * U[r, j];

}

U[k, j] = (A[k, j] - temp) / L[k, k];

}

}

（2）计算Y矩阵

计算Y矩阵主要根据式（5-8），其代码如下：

//计算Y矩阵

double[] Y = new double[arrayCount];

for (i = 0; i < arrayCount; i++)

{

temp = 0.0;

for (j = 0; j < i; j++)

{

temp += L[i, j] * Y[j];

}

Y[i]=(B[i]-temp)/L[i,i];

}

（3）计算X矩阵

计算X矩阵主要根据式（5-9），其代码如下：

//计算X矩阵

double[] XX = new double[arrayCount];

for (i = arrayCount-1; i >=0 ; i--)

{

temp = 0.0;

for (j = i; j < arrayCount; j++)

{

temp += U[i, j] * XX[j];

}

XX[i] = Y[i] - temp;

}

至此，整个Crout分解法的函数如下：

private void CalFoundation2(double[,] A, out double[] X, double[] B)

{

int arrayCount = B.Length;//矩阵的行、列数

int i, j, k, r;

double temp;

//计算LU矩阵

double[,] L = new double[arrayCount, arrayCount];

double[,] U = new double[arrayCount, arrayCount];

for (i = 0; i < arrayCount; i++)

{

U[i, i] = 1;

}

for (k = 0; k < arrayCount; k++)

{

for (i = k; i < arrayCount; i++)

{

temp = 0.0;

for (r = 0; r < k ; r++)

{

temp+=L[i,r]*U[r,k];

}

L[i, k] = A[i, k] - temp;

}

for (j = k+1; j < arrayCount; j++)

{

temp = 0.0;

for (r = 0; r < k ; r++)

{

temp += L[k, r] * U[r, j];

}

U[k, j] = (A[k, j] - temp) / L[k, k];

}

}

//计算Y矩阵

double[] Y = new double[arrayCount];

for (i = 0; i < arrayCount; i++)

{

temp = 0.0;

for (j = 0; j < i; j++)

{

temp += L[i, j] * Y[j];

}

Y[i]=(B[i]-temp)/L[i,i];

}

//计算X矩阵

double[] XX = new double[arrayCount];

for (i = arrayCount-1; i >=0 ; i--)

{

temp = 0.0;

for (j = i; j < arrayCount; j++)

{

temp += U[i, j] * XX[j];

}

XX[i] = Y[i] - temp;

}

X = XX;

}