3. 直接三角分解法(LU Decomposition)

下面是使用Python实现直接三角分解法(LU Decomposition)求解线性方程组的代码。LU分解将矩阵分解为一个下三角矩阵 ( L ) 和一个上三角矩阵 ( U )。求解时,可以通过先解 ( Ly = b ) 然后解 ( Ux = y ) 来得到最终解 ( x )。

实现步骤

  1. LU分解:将矩阵 ( A ) 分解为下三角矩阵 ( L ) 和上三角矩阵 ( U )。

  2. 前向替代(Forward Substitution):解决方程 ( Ly = b )。

  3. 回代(Backward Substitution):解决方程 ( Ux = y )。

代码实现

import numpy as np

def lu_decomposition(A):
    """
    Perform LU decomposition of matrix A.

    Parameters:
    A : 2D list or numpy array
        Coefficient matrix

    Returns:
    L : numpy array
        Lower triangular matrix
    U : numpy array
        Upper triangular matrix
    """
    n = len(A)
    L = np.zeros((n, n))
    U = np.zeros((n, n))

    for i in range(n):
        # Upper Triangular
        for k in range(i, n):
            U[i, k] = A[i, k] - sum(L[i, j] * U[j, k] for j in range(i))
        
        # Lower Triangular
        for k in range(i, n):
            if i == k:
                L[i, i] = 1  # Diagonal as 1
            else:
                L[k, i] = (A[k, i] - sum(L[k, j] * U[j, i] for j in range(i))) / U[i, i]

    return L, U

def forward_substitution(L, b):
    """
    Perform forward substitution to solve Ly = b.

    Parameters:
    L : 2D list or numpy array
        Lower triangular matrix
    b : 1D list or numpy array
        Right-hand side vector

    Returns:
    y : numpy array
        Solution vector
    """
    n = len(b)
    y = np.zeros(n)

    for i in range(n):
        y[i] = (b[i] - sum(L[i, j] * y[j] for j in range(i))) / L[i, i]

    return y

def backward_substitution(U, y):
    """
    Perform backward substitution to solve Ux = y.

    Parameters:
    U : 2D list or numpy array
        Upper triangular matrix
    y : 1D list or numpy array
        Solution vector from forward substitution

    Returns:
    x : numpy array
        Solution vector
    """
    n = len(y)
    x = np.zeros(n)

    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - sum(U[i, j] * x[j] for j in range(i + 1, n))) / U[i, i]

    return x

def solve_lu_decomposition(A, b):
    """
    Solve the system of linear equations Ax = b using LU Decomposition.

    Parameters:
    A : 2D list or numpy array
        Coefficient matrix
    b : 1D list or numpy array
        Right-hand side vector

    Returns:
    x : numpy array
        Solution vector
    """
    L, U = lu_decomposition(A)
    y = forward_substitution(L, b)
    x = backward_substitution(U, y)
    return x

# Example usage
A = np.array([[2, -1, 1],
              [3, 3, 9],
              [3, 3, 5]], dtype=float)
b = np.array([2, -1, 2], dtype=float)

x = solve_lu_decomposition(A, b)
print("Solution:", x)

代码说明

  1. LU分解(lu_decomposition函数)

    • 将矩阵 ( A ) 分解为下三角矩阵 ( L ) 和上三角矩阵 ( U )。

  2. 前向替代(forward_substitution函数)

    • 使用下三角矩阵 ( L ) 和右端向量 ( b ) 求解中间变量 ( y )。

  3. 回代(backward_substitution函数)

    • 使用上三角矩阵 ( U ) 和中间变量 ( y ) 求解最终解 ( x )。

  4. 求解线性方程组(solve_lu_decomposition函数)

    • 综合上述步骤,完成LU分解并求解线性方程组。

示例

给定系数矩阵 (A) 和右端向量 (b):

A = [[2, -1, 1],
     [3, 3, 9],
     [3, 3, 5]]
b = [2, -1, 2]

运行上述代码,求解线性方程组 (Ax = b),可以得到解向量 x

Last updated