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;

        }