高斯消去法

高斯消去法

      1、高斯消去法

       以前在学习线性代数时，经常通过初等变换将一个矩阵化为行阶梯形矩阵，然后再进行

各种处理等。最近使用高斯消去法求解线性方程组时，对矩阵也要进行这种处理。将方程组

对应的增广矩阵化为行阶梯型矩阵后，就能很方便地求出方程组的解。

      废话就不多说了，下面先贴上代码(java实现)

public class GaussianElimination {
	/**
	 * Get solutions of the equations using Gaussian elimination.
	 * @param a - Coefficient matrix of the equations
	 * @param b - right-hand side of the equations
	 * @return solution of the equations
	 */
	public static double[] solve(double[][] a, double[] b) {		
		elimination(a, b);		//Eliminate the matrix
		return backSubstitution(a, b);   //solve the equations
	}

	/**
	 * Solve the equations
	 * @param a - Coefficient matrix of the equations
	 * @param b - constant matrix of the equations
	 * @return root of the equations
	 */
	private static double[] backSubstitution(double[][] a, double[] b) {
		double[] result = new double[a[0].length];
		for(int i=a.length-1; i>=0; i--) {
			double equationRight = b[i];
			for(int j=i+1; j<a[i].length; j++) {
				equationRight -= a[i][j] * result[j];
			}
			result[i] = equationRight / a[i][i];
		}
		return result;
	}
	
	/**
	 * Eliminate the matrix to simplify it
	 * @param a - Coefficient matrix of the equations
	 * @param b - right-hand side of the equations
	 */
	private static void elimination(double[][] a, double[] b) {
		for(int j=0; j<a[0].length - 1; j++) {	
			if(a[j][j] == 0) {
				throw new IllegalArgumentException("zero pivot encountered.");
			}
			
			for(int i=j+1; i<a.length; i++) {
				double mult = a[i][j] / a[j][j]; 
				for(int k=j; k<a[i].length; k++) {
					a[i][k] = a[i][k] - a[j][k] * mult;
				}
				b[i] = b[i] - b[j] * mult;
			}
		}
	}
}

     上面的函数   elimination 就做了一件事情，将一个矩阵通过矩阵初等变换转化成行阶梯形矩阵。

但是当方程组对应的系数矩阵对角线上的元素中有0时，高斯消去法就失效了，待会再讨论改

进方法。将矩阵化成行阶梯形矩阵后，矩阵就是一个上三角矩阵了，这样就可以很方便地

求出方程的解。 而函数  backSubstitution 则完成了回代过程，从上三角的最下端的 Xn开始，

依次求出Xn, Xn-1, ....x2, x1(每一次求解都会用到之前求出的解)。

      2、改进

      上面讲到当线性方程组对应的系数矩阵中对角线元素出现0时，高斯消去法就会失效。

这里，我们将通过在原有的经典高斯消去法基础上加上部分选元法来解决这一问题。

      2.1部分选元法

        部分消元法在执行每一步消元步骤之前都要比较这一列上的对角线元素及该列该对角线

元素以下的元素的大小，确定最大元素所在行，并将该行与该列对角线上元素所在行交换。

进行交换之后，再与经典高斯消去法一样进行消元。这样就可以解决经典高斯消去法的

局限了。

      举个例子如下：

      对于方程组：

                          

      其对应的 增广矩阵为

                     

      根据高斯消去法，进行第一次消元时，需要比较 |A11|=1 与 |A21|=1及 |A31|=2,选择 A31

作为新的主元，这通过交换第一行与第三行即可实现。交换之后，第二行减去 (-1/2) 乘以

第一行，第三行减去 (1/2) 乘以第一行，即完成了第一次消元。

        完成第一次消元后，矩阵变成了

                    

        后面的消元与这类似，就不再详述了。

        2.2 程序实现

         在上面的经典高斯消去法代码中进行消去之前，对现在的矩阵进行行变换，使得每一次

进行消元之前新的主元出现在对角线上。实现这一过程的函数代码如下:

/**
	 * Select max abs element as current partial pivot
	 * @param a - Coefficient matrix of the equations
	 * @param b - constant matrix of the equations
	 * @param j - index of current row to eliminate
	 */
	private static void getMaxAbsPivot(double[][] a, double[] b, int j) {
		int maxIndex = j;
		double maxValue = Math.abs(a[j][j]);
		for(int i=j+1; i<a.length; i++) {
			if(Math.abs(a[i][j]) > maxValue) {
				maxIndex = i;
				maxValue = Math.abs(a[i][j]);
			}
		}
		
		if(maxIndex != j) {
			for(int m=0; m<a[0].length; m++) {
				double temp = a[j][m];
				a[j][m] = a[maxIndex][m];
				a[maxIndex][m] = temp; 
			}
			double tempB = b[j];
			b[j] = b[maxIndex];
			b[maxIndex] = tempB;
		}
	}	

       3.后记

              上面的部分选元法解决了方程的系数矩阵的对角线上元素为0时，高斯消去法失效的问题,

       但是每进行一次消元 ，都要进行多次矩阵初等行变换 ，效率有点低。对于一般要求，当遇到

      为0的主元时，我们可以简单的将这个主元的值设置为一个很小的小数，例如 1e-20。反映到

      代码上来，就是将上面的 elimination 函数中 

                                       if(a[j][j] == 0) {.....}

      这一条语句改为   

                                      if(a[j][j] == 0) { a[j][j] = 1e-20; }

              这样不但对解的精确度的影响几乎可以忽略，而且还省去了部分选元法中的初等行变换。