課程筆記

C10 LU Decomposition and Matrix Inversion (LU 分解與矩陣求逆)

10.1 LU Decomposition (LU 分解)

原理：將矩陣分解為兩個矩陣與的乘積。
- [A]{X}-{B}=0[U]{X}-{D}=0
- - [L]{[U]{X}={D}}=[A]{X}-{B}
  - [L][U]=[A]；[L]{D}={B}
- ：下三角矩陣 (Lower Triangular Matrix)
- ：上三角矩陣 (Upper Triangular Matrix)
步驟：
1. 透過高斯消去法，將矩陣前向消去轉化為上三角矩陣。
  - $前向消去$
2. 透過前向消去的過程，將消去的係數存入。
3. 透過回代法，求解。
舉例：
- 解題步驟：
  1. 前向消去找到與。
    - - $前向消去$
      - $前向消去$
    - 可知：; ;
    - 得到：、
  2. 驗證。
  3. 回代，得到：
  4. 回代：

10.2 Matrix Inversion (反矩陣)

定義：
- 其中，為單位矩陣。
方法：
1. 當時：
  - ，求得的第一行。
  - ，求得的第二行。
  - ，求得的第三行。
2. 此時，。
舉例(延續前面的)：求的反矩陣。
1. 透過 LU 分解，求得與。
2. 透過反矩陣的定義，求得：
  1. ：
    - 計算得：
    1. - 計算得：
    2. 此時，的第一行：
  2. ：
    2. 計算得：
    4. 此時，的第二行：
  3. 1. 此時，的第三行：
    - 得：

C11 Special Matrices and Gauss-Seidel Method (特殊矩陣與高斯-賽德法)

11.2 Gauss-Seidel Method (高斯-賽德法)

原理：
- 對於一三次方程組：
- 可將其移項為 [A]{X}={B} 的形式：
- 透過迭代法，將、、逐步求解 (如同固定點迭代法)。
- 重複迭代直至收斂相對百分比誤差小於設定值。
範例：
1. 將方程式轉化為 [A]{X}={B} 的形式：
2. 重複迭代直至收斂相對百分比誤差小於設定值 (Gauss-Seidel Method)：
  1. 迭代第 1 次：
    1. 假設、初值為 0，得到新的：
    2. 帶入初值為 0，以及剛才獲得的，得到新的：
    3. 帶入剛才獲得的、，得到新的：
  2. 迭代第 2 次：
    1. 新的：
    2. 新的：
    3. 新的：
3. 重複迭代直至收斂相對百分比誤差小於設定值 (Jacobi iteration)：
  1. 迭代第 1 次：
    1. 假設、初值為 0，得到新的：
    2. 假設、初值為 0，得到新的：
    3. 假設、初值為 0，得到新的：
  2. 迭代第 2 次：
    1. 新的：
    2. 新的：
    3. 新的：
收斂條件：
- 方程式對角線係數值必須大於其他元素絕對值之總和（充分條件，但非必要條件）：

當收斂慢或不收斂時，可透過鬆弛 (Relaxation) 來加速收斂：
- 鬆弛：
  - 其中，為鬆弛因子，介於 [0, 2]。當介於 [0, 1]，則為 Underrelaxation，減緩收斂速度 (可能讓一些情況可以收斂)；而若介於 [1, 2]，則為 Overrelaxation，加速收斂速度；若恰為 1，則為 Gauss-Seidel Method。

Homework

題目

使用 LU 分解、Gauss-Seidel、Jacobi、Gauss-Siedel with Relaxation 四種方法求解。

解題流程

LU 分解：
1. 透過高斯消去法，求得與。
2. 驗證。
3. 回代，求得。
4. 回代，求得。
Gauss-Seidel：
1. 類似於固定點迭代法，先將方程式轉換為的形式，在迭代過程中求得。
Jacobi：
1. 類似於 Gauss-Seidel，但是在迭代過程中，會同時更新。
Gauss-Siedel with Relaxation：
1. 透過鬆弛因子，調整收斂速度。
- x[i] = lambda*x_new + (1-lambda)*x_old[i]

Code

針對本作業撰寫 Golang 代碼如下所示。其中，定義 calcLUDecompositionMethod()、calcGaussSeidelMethod()、calcJacobiMethod()、calcGaussSeidelMethodRelaxation() 四種方法分別使用 LU Decomposition、Gauss-Seidel、Jacobi、Gauss-Siedel with Relaxation 四種數值方法求解聯立方程組。
本代碼在三個數值方法中先透過 copyMatrix()、copyVector() 方法複製輸入的聯立方程組避免原始矩陣被干擾，並且本代碼建立 printStep()、printStepMatrix() 方法並建立區域變數 Step 紀錄並列印每個執行步驟。
本代碼在 calcGaussSeidelMethod()、calcJacobiMethod()、calcGaussSeidelMethodRelaxation() 方法中添加了 rearrangeMatrix() 方法。主要原因為題目式子存在對角線為零的問題，若直接將其轉換為形式將導致於第二行存在除以零問題使算法無解。因此，本作業透過實踐 rearrangeMatrix() 方法盡可能將對角線調整為絕對值最大的項目。不僅可避免除以零問題，也可以進一步提升方程組被數值方法解決的機會。
而 Homework() 方法定義了於本作業中欲求解的方程組係數矩陣 A、常數矩陣 b，並分別呼叫上述四種方法求解。其中，calcGaussSeidelMethod()、calcJacobiMethod()、calcGaussSeidelMethodRelaxation() 方法於本作業中定義使用相對百分誤差小於 0.5% 作為停止標準。且 calcGaussSeidelMethodRelaxation() 計算中定義本作業 λ 倍數為 1.05。最後，透過 main() 函數呼叫 Homework() 函數執行。

package main
import (
    "fmt"
    "math"
    "time"
)
 
func copyMatrix(A [][]float64) [][]float64 {
    A_Copy := make([][]float64, len(A))
    for i := range A {
        A_Copy[i] = make([]float64, len(A[i]))
        for j := range A[i] {
            A_Copy[i][j] = A[i][j]
        }
    }
    return A_Copy
}
 
func copyVector(b []float64) []float64 {
    b_Copy := make([]float64, len(b))
    copy(b_Copy, b)
    return b_Copy
}
 
func printStep(A [][]float64, b []float64, stepDescription string) {
    fmt.Println(stepDescription)
    n := len(A)
    for i := 0; i < n; i++ {
        for j := 0; j < n; j++ {
            fmt.Printf("%8.4f ", A[i][j])
        }
        fmt.Printf("| %8.4f\n", b[i])
    }
    fmt.Println()
}
 
func printStepMatrix(matrix [][]float64, description string) {
    fmt.Println(description)
    n := len(matrix)
    for i := 0; i < n; i++ {
        for j := 0; j < len(matrix[i]); j++ {
            fmt.Printf("%8.4f ", matrix[i][j])
        }
        fmt.Println()
    }
    fmt.Println()
}
 
func rearrangeMatrix(A [][]float64, b []float64) {
    n := len(A)
    for i := 0; i < n; i++ {
        maxRow := i
        maxValue := math.Abs(A[i][i])
        for k := i + 1; k < n; k++ {
            if math.Abs(A[k][i]) > maxValue {
                maxValue = math.Abs(A[k][i])
                maxRow = k
            }
        }
        if maxRow != i {
            A[i], A[maxRow] = A[maxRow], A[i]
            b[i], b[maxRow] = b[maxRow], b[i]
            fmt.Printf("交換第 %d 行和第 %d 行，以使對角線元素最大。\n", i+1, maxRow+1)
        }
    }
}
 
func calcApproxError(Current, Previous float64) float64 {
    if Current == 0 {
        return math.Inf(1)
    }
    return math.Abs((Current - Previous) / Current)
}
 
func calcLUDecompositionMethod(A [][]float64, b []float64) []float64 {
    start := time.Now() // 記錄開始時間
    defer func() {
        fmt.Printf("運行時間: %v\n", time.Since(start))
    }()
 
    A_Copy := copyMatrix(A)
    b_Copy := copyVector(b)
    n := len(A_Copy)
    x := make([]float64, n)
 
    // 初始化 L 和 U 矩陣為零矩陣
    L := make([][]float64, n)
    U := make([][]float64, n)
    for i := 0; i < n; i++ {
        L[i] = make([]float64, n)
        U[i] = make([]float64, n)
    }
 
    // LU 分解過程
    fmt.Println("LU 分解過程開始...")
    for k := 0; k < n; k++ {
        // 計算 U 矩陣的第 k 行元素
        for j := k; j < n; j++ {
            sum := 0.0
            for s := 0; s < k; s++ {
                sum += L[k][s] * U[s][j]
            }
            U[k][j] = A_Copy[k][j] - sum
            fmt.Printf("U_%d%d = A_%d%d", k+1, j+1, k+1, j+1)
            for s := 0; s < k; s++ {
                fmt.Printf(" - L_%d%d * U_%d%d", k+1, s+1, s+1, j+1)
            }
            fmt.Printf(" = %.4f\n", U[k][j])
        }
 
        // 計算 L 矩陣的第 k 列元素
        for i := k; i < n; i++ {
            if i == k {
                L[i][k] = 1.0 // L 矩陣的對角線元素為 1
                fmt.Printf("L_%d%d = 1.0\n", i+1, k+1)
            } else {
                sum := 0.0
                for s := 0; s < k; s++ {
                    sum += L[i][s] * U[s][k]
                }
                L[i][k] = (A_Copy[i][k] - sum) / U[k][k]
                fmt.Printf("L_%d%d = (A_%d%d", i+1, k+1, i+1, k+1)
                for s := 0; s < k; s++ {
                    fmt.Printf(" - L_%d%d * U_%d%d", i+1, s+1, s+1, k+1)
                }
                fmt.Printf(") / U_%d%d = %.4f\n", k+1, k+1, L[i][k])
            }
        }
 
        // 列印當前的 L 和 U 矩陣
        fmt.Printf("\n步驟 %d: 對第 %d 列進行 LU 分解後的 L 和 U 矩陣\n", k+1, k+1)
        printStepMatrix(L, "L 矩陣：")
        printStepMatrix(U, "U 矩陣：")
    }
 
    // 前向替代，求解 Ld = b
    d := make([]float64, n)
    fmt.Println("前向替代過程開始")
    for i := 0; i < n; i++ {
        sum := b_Copy[i]
        fmt.Printf("d_%d = b_%d", i+1, i+1)
        for j := 0; j < i; j++ {
            sum -= L[i][j] * d[j]
            fmt.Printf(" - L_%d%d * d_%d", i+1, j+1, j+1)
        }
        d[i] = sum
        fmt.Printf(" = %.4f\n", d[i])
    }
    fmt.Printf("d = %v\n", d)
 
    // 後向替代，求解 Ux = d
    fmt.Println("\n後向替代過程開始...")
    for i := n - 1; i >= 0; i-- {
        sum := d[i]
        fmt.Printf("x_%d = (d_%d", i+1, i+1)
        for j := i + 1; j < n; j++ {
            sum -= U[i][j] * x[j]
            fmt.Printf(" - U_%d%d * x_%d", i+1, j+1, j+1)
        }
        x[i] = sum / U[i][i]
        fmt.Printf(") / U_%d%d = %.4f\n", i+1, i+1, x[i])
    }
    fmt.Printf("x (Solution) = %v\n", x)
 
    return x
}
 
func calcGaussSeidelMethod(A [][]float64, b []float64, tolerance float64) []float64 {
    start := time.Now() // 記錄開始時間
    defer func() {
        fmt.Printf("運行時間: %v\n", time.Since(start))
    }()
 
    A_Copy := copyMatrix(A)
    b_Copy := copyVector(b)
    n := len(A_Copy)
    x := make([]float64, n) // 初始解向量，全零
 
    // 重新排列矩陣和向量，使對角線元素最大化
    fmt.Println("對矩陣進行重新排列，使對角線元素最大化...")
    rearrangeMatrix(A_Copy, b_Copy)
    printStep(A_Copy, b_Copy, "重新排列後的矩陣 A 和向量 b：")
 
    fmt.Println("Gauss-Seidel 方法迭代過程開始...")
 
    iteration := 0
    for {
        iteration++
        x_old := copyVector(x)
        fmt.Printf("第 %d 次迭代:\n", iteration)
        for i := 0; i < n; i++ {
            sum := b_Copy[i]
            fmt.Printf("x_%d = (b_%d", i+1, i+1)
            for j := 0; j < n; j++ {
                if j != i {
                    sum -= A_Copy[i][j] * x[j]
                    fmt.Printf(" - A_%d%d * x_%d", i+1, j+1, j+1)
                }
            }
            x[i] = sum / A_Copy[i][i]
            fmt.Printf(") / A_%d%d = %.4f\n", i+1, i+1, x[i])
        }
 
        // 計算誤差
        errors := make([]float64, n)
        for i := 0; i < n; i++ {
            errors[i] = calcApproxError(x[i], x_old[i])
        }
 
        fmt.Printf("Approximation errors: %.5v\n", errors)
 
        // 判斷是否收斂
        maxError := 0.0
        for _, err := range errors {
            if err > maxError {
                maxError = err
            }
        }
        if maxError < tolerance {
            fmt.Println("收斂達到指定的容許誤差。")
            break
        }
    }
 
    fmt.Printf("x (Solution) = %.10v\n", x)
    return x
}
 
func calcJacobiMethod(A [][]float64, b []float64, tolerance float64) []float64 {
    start := time.Now() // 記錄開始時間
    defer func() {
        fmt.Printf("運行時間: %v\n", time.Since(start))
    }()
 
    A_Copy := copyMatrix(A)
    b_Copy := copyVector(b)
    n := len(A_Copy)
    x := make([]float64, n)     // 當前解向量
    x_old := make([]float64, n) // 前一次的解向量
 
    // 重新排列矩陣和向量，使對角線元素最大化
    fmt.Println("對矩陣進行重新排列，使對角線元素最大化...")
    rearrangeMatrix(A_Copy, b_Copy)
    printStep(A_Copy, b_Copy, "重新排列後的矩陣 A 和向量 b：")
 
    fmt.Println("Jacobi 方法迭代過程開始...")
 
    iteration := 0
    for {
        iteration++
        copy(x_old, x)
        fmt.Printf("第 %d 次迭代:\n", iteration)
        for i := 0; i < n; i++ {
            sum := b_Copy[i]
            fmt.Printf("x_%d = (b_%d", i+1, i+1)
            for j := 0; j < n; j++ {
                if j != i {
                    sum -= A_Copy[i][j] * x_old[j]
                    fmt.Printf(" - A_%d%d * x_%d", i+1, j+1, j+1)
                }
            }
            x[i] = sum / A_Copy[i][i]
            fmt.Printf(") / A_%d%d = %.4f\n", i+1, i+1, x[i])
        }
 
        // 計算誤差
        errors := make([]float64, n)
        for i := 0; i < n; i++ {
            errors[i] = calcApproxError(x[i], x_old[i])
        }
 
        fmt.Printf("Approximation errors: %.5v\n", errors)
 
        // 判斷是否收斂
        maxError := 0.0
        for _, err := range errors {
            if err > maxError {
                maxError = err
            }
        }
        if maxError < tolerance {
            fmt.Println("收斂達到指定的容許誤差。")
            break
        }
    }
 
    fmt.Printf("x (Solution) = %v\n", x)
    return x
}
 
func calcGaussSeidelMethodRelaxation(A [][]float64, b []float64, lambda float64, tolerance float64) []float64 {
    start := time.Now() // 記錄開始時間
    defer func() {
        fmt.Printf("運行時間: %v\n", time.Since(start))
    }()
 
    A_Copy := copyMatrix(A)
    b_Copy := copyVector(b)
    n := len(A_Copy)
    x := make([]float64, n) // 初始解向量，全零
 
    // 重新排列矩陣和向量，使對角線元素最大化
    fmt.Println("對矩陣進行重新排列，使對角線元素最大化...")
    rearrangeMatrix(A_Copy, b_Copy)
    printStep(A_Copy, b_Copy, "重新排列後的矩陣 A 和向量 b：")
 
    fmt.Println("帶鬆弛因子的 Gauss-Seidel 方法迭代過程開始...")
 
    iteration := 0
    for {
        iteration++
        x_old := copyVector(x) // 用於計算誤差
        fmt.Printf("第 %d 次迭代:\n", iteration)
        for i := 0; i < n; i++ {
            sum := b_Copy[i]
            fmt.Printf("x_%d = (b_%d", i+1, i+1)
            for j := 0; j < n; j++ {
                if j != i {
                    sum -= A_Copy[i][j] * x[j]
                    fmt.Printf(" - A_%d%d * x_%d", i+1, j+1, j+1)
                }
            }
            x_new := sum / A_Copy[i][i]
            // 應用鬆弛因子進行更新
            x[i] = lambda*x_new + (1-lambda)*x_old[i]
            fmt.Printf(") / A_%d%d = %.4f, 經過鬆弛更新後的 x_%d = %.4f\n", i+1, i+1, x_new, i+1, x[i])
        }
 
        // 計算誤差
        errors := make([]float64, n)
        for i := 0; i < n; i++ {
            errors[i] = calcApproxError(x[i], x_old[i])
        }
 
        fmt.Printf("Approximation errors: %.5v\n", errors)
 
        // 判斷是否收斂
        maxError := 0.0
        for _, err := range errors {
            if err > maxError {
                maxError = err
            }
        }
        if maxError < tolerance {
            fmt.Println("收斂達到指定的容許誤差。")
            break
        }
    }
 
    fmt.Printf("x (Solution) = %v\n", x)
    return x
}
 
func Homework() {
    // A= [2, 4, -1, -2; 4, 0, 2, 1; 1, 3, -2, 0; 3, 2, 0, 5]
    // b= [10; 7; 3; 2]
 
    A := [][]float64{{2, 4, -1, -2}, {1, 3, -2, 0}, {4, 0, 2, 1}, {3, 2, 0, 5}}
    b := []float64{10, 3, 7, 2}
    printStep(A, b, "Original Matrix A and b:")
 
    fmt.Println("----------------------------------------")
    fmt.Println("LU Decomposition Method:")
    calcLUDecompositionMethod(A, b)
    fmt.Println("----------------------------------------")
    fmt.Println("Gauss-Seidel Method:")
    calcGaussSeidelMethod(A, b, 0.005)
    fmt.Println("----------------------------------------")
    fmt.Println("Jacobi Method:")
    calcJacobiMethod(A, b, 0.005)
    fmt.Println("----------------------------------------")
    fmt.Println("Gauss-Seidel Method with Relaxation:")
    lambda := 1.05
    fmt.Printf("Relaxation factor λ = %.2f\n", lambda)
    calcGaussSeidelMethodRelaxation(A, b, lambda, 0.005)
}
 
func main() {
    Homework()
    fmt.Scanln()
}

< Last

Next >

c55Math4833

儲存庫

4. LU分解與矩陣求逆