Gauss消元法的原理及Java實現

補充知識：

• 正定矩陣
• 奇異矩陣
• 嚴格對角占優

 

要理解Gauss消去法，首先來看一個例子：

從上例子可以看出，高斯消去法實際上就是我們初中學的階二元一次方程組，只不過那里的未知數個數$n=2$

$n>2$時，Gauss消去法的思路實際上和解二元一次方程組是一樣的，方法如下：

• 將$n$方程組中的$n-1$個方程通過消元，形成一個與原方程組等價的一個新方程組，新方程組中的$n-1$個方程僅包含$n-1$個未知數。
• 故問題就轉化為了求解$n-1$元的方程組，這樣我們可以繼續消元，以次類推，直到最后一個方程組為一元一次方程組
• 從最后一個一元一次方程組求解出最后一個未知量，然后逐步回代入之前的方程組，從而得到所有的未知數。
• 我們可以看到Gauss實際上就分為兩步：消去和回代

下面通過一般化得到Gauss消元法的求解過程

 

 

以上就是Gauss消去法的基本步驟，我們再回過頭看看有沒有什么問題？

我們在求比例$l_{ik}= \frac{a_{ik}^{\left (k-1  \right )}}{a_{kk}^{\left (k-1  \right )}}$時，如果分母很小，即：

 

$l_{ik}\rightarrow \infty$，那么

總結一下，能否使用Gauss消元法的情況

為了解決這個問題，我們可以使用列主元Gauss消元法。

參考了一些網上的代碼，這里給出Gauss的Java實現

 

package peterxiazhe;

import java.util.Scanner;

public class Gauss {
    /**
     * 列主元高斯消去法
     */
    static double A[][];
    static double b[];
    static double x[];
    
    static int n;    //n表示未知數的個數
    static int n_2;    //記錄換行的次數
    
    public static void main(String[] args) {
        System.out.println("--------------輸入方程組未知數的個數---------------");
        Scanner sc = new Scanner(System.in);
        n = sc.nextInt();
        
        A = new double[n][n];
        b = new double[n];
        x = new double[n];
        
        System.out.println("--------------輸入方程組的系數矩陣A:---------------");
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                A[i][j] = sc.nextDouble();
            }
        }
        
        System.out.println("--------------輸入方程組的常量向量b:---------------");
        for(int i = 0; i < n; i++) {
                b[i] = sc.nextDouble();
            }
        
        Elimination();
        BackSubstitution();
        PrintRoot();
    }
    
    
    //消元法
    public static void Elimination() {
        PrintA();
        for(int k = 0; k < n; k++) {
            WrapRow(k);
            for(int i = k+1; i < n; i++) {
                double l = A[i][k] / A[k][k];
                A[i][k] = 0;
                
                for(int j = k+1; j < n; j++) {
                    A[i][j] = A[i][j] - l * A[k][j];
                }
                b[i] = b[i] - l * b[k];
            }
            //System.out.println("第" + k + "次消元后:");
            //PrintA();
        }
    }
    
    //回代法
    public static void    BackSubstitution() {
        x[n-1] = b[n-1] / A[n-1][n-1];
        for(int i = n - 2; i >= 0; i--) {
            x[i] = (b[i] - solve(i)) / A[i][i];
        }
    }
    
    public static double solve(int i) {
        double result = 0.0;
        for(int j = i; j < n; j++)
            result += A[i][j] * x[j];
        return result;
    }
    
    
    //輸出方程組的根
    public static void PrintRoot() {
        System.out.println("--------------方程組的根為---------------");
        for(int  i = 0; i < n; i++) {
            System.out.println("x" + (i+1) + " = " + x[i]);
        }
    }
    
    //交換Swap函數???
    public static void Swap(double[] ar, int x, int y) {
        Double tmp = ar[x];
        ar[x] = ar[y];
        ar[y] = tmp;
    }
    
    public static void PrintA() {    //輸出A的增廣矩陣
        //System.out.println("--------------增廣矩陣---------------");
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                System.out.print(A[i][j] + " ");
            }
            System.out.println(b[i]);
        }
    }
    
    //交換矩陣的行
    public static void WrapRow(int k) {    //k表示第k+1輪消元
        double maxElement = Math.abs(A[k][k]);
        
        int WrapRowIndex = k;    //    記住要交換的行
        for(int i = k + 1; i < n; i++) {
            if (Math.abs(A[i][k]) > maxElement) {
                WrapRowIndex = i;
                maxElement = A[i][k];
            }
        }
        if (WrapRowIndex != k) {    //交換求得最大主元
            n_2 += 1;
            System.out.println("k = " + k + "時，" + "要交換的行為" + k + "和"+ WrapRowIndex);
            
            //先交換A
            for(int j = k; j < n; j++) {
                double[] arr = {A[k][j], A[WrapRowIndex][j]};
                Swap(arr, 0, 1);
                A[k][j] = arr[0]; A[WrapRowIndex][j] = arr[1];
//                double tmp = A[k][j];
//                A[k][j] = A[WrapRowIndex][j];
//                A[WrapRowIndex][j] = tmp;
            }
            
            //再交換b
            double[] arr = {b[k], b[WrapRowIndex]};
            Swap(arr, 0, 1);
            b[k] = arr[0]; b[WrapRowIndex] = arr[1];
//            double tmp = b[k];
//            b[k] = b[WrapRowIndex];
//            b[WrapRowIndex] = tmp;
            System.out.println("--------------交換后---------------");
            PrintA();
        }        
    }
}

 

 

注意：由於Java不支持對基本數據類型的引用傳遞，這里使用了一個小技巧。