高斯消元法求方程组的解

高斯消元法（也叫行化简法），简单地来说就是构造方程组的增广矩阵，将增广矩阵化为阶梯型，将每一行主元位置化为1，主元下面元素都为0。最后从下往上依次代入求解。

至于主元的选择，最简单的就是从第一行开始依次选取主对角线上的元素，但是这样做的有一点需要注意就是如果主元非常小，就可能除了之后产生非常大的数字，可能产生精度的损失。所以每一次选取一列中最大的数字作为主元，通过行交换把改数字放在主对角线上。

public class GaussProgram {

    /// <summary>
    /// selectMainElement 用来选取第k个主元 其所在列下面所有元素中最大的元素作为主元
    /// </summary>
    /// <param name="n">传递过来的矩阵a的行数</param>
    /// <param name="k">现在主元所在的行数</param>
    /// <param name="a">传递过来用来处理的矩阵</param>
    public static void selectMainElement(int n, int k, double[,] a)
    {
        double t, mainElement;
        int l; // 用于保存主元所在的行号
        mainElement = Math.Abs(a[k, k]);

        // 从第k行到第n行寻找第k列的主元素，记下主元素mainElement和所在的行号l
        l = k;
        for (int i = k; i < n; ++i) 
        {
            if (mainElement < Math.Abs (a[i, k])) 
            {
                mainElement = Math.Abs (a[i, k]);
                l = i;
            }
        }

        // 将l行与k行交换，每行前面的k个元素都是0，不必交换
        if (l != k) 
        {
            for (int j = k; j <= n; ++j) 
            {
                t = a [k, j];
                a[k, j] = a[l, j];
                a[l, j] = t;
            }
        }
    }

    /// <summary>
    /// Gauss 列主消元法求n元一次方程组的解
    /// </summary>
    /// <param name="n">传递过来的矩阵a的行数</param>
    /// <param name="a">方程组的增广矩阵</param>
    public static double[] Gauss(int n, double[,] a)
    {
        double[] x = new double[n];
        // 消元过程 最外层循环是找主元
        for (int k = 0; k < n - 1; ++k) 
        {
            // 找主元并交换矩阵相关行
            selectMainElement (n, k, a);

            for (int i = k; i < n - 1; ++i) // i为行号，每次循环将主元素下面的所有元素化为0
            { 
                double m = a[i + 1, k] / a[k, k]; // m为要化为0需要给主元乘上的因子
                for (int j = k; j <= n; ++j) 
                {
                    a [i + 1, j] = a [i + 1, j] - m * a [k, j];
                }
            }
        }

        // 回代求解过程
        for (int k = n - 1; k >= 0; --k) 
        {
            double addResult = 0.0;
            for (int j = k + 1; j < n; ++j) 
            {
                addResult += x[j] * a[k, j]; // k为行数，j为列数，从所求未知数的后一列开始
            }
            x[k] = (a[k, n] - addResult) / a[k, k];
        }
        return x;
    }
}