下面是使用Python实现直接三角分解法(LU Decomposition)求解线性方程组的代码。LU分解将矩阵分解为一个下三角矩阵 ( L ) 和一个上三角矩阵 ( U )。求解时,可以通过先解 ( Ly = b ) 然后解 ( Ux = y ) 来得到最终解 ( x )。
实现步骤
LU分解:将矩阵 ( A ) 分解为下三角矩阵 ( L ) 和上三角矩阵 ( U )。
前向替代(Forward Substitution):解决方程 ( Ly = b )。
回代(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)
代码说明
LU分解(lu_decomposition函数):
将矩阵 ( A ) 分解为下三角矩阵 ( L ) 和上三角矩阵 ( U )。
前向替代(forward_substitution函数):
使用下三角矩阵 ( L ) 和右端向量 ( b ) 求解中间变量 ( y )。
回代(backward_substitution函数):
使用上三角矩阵 ( U ) 和中间变量 ( y ) 求解最终解 ( x )。
求解线性方程组(solve_lu_decomposition函数):
综合上述步骤,完成LU分解并求解线性方程组。
示例
给定系数矩阵 (A) 和右端向量 (b):
A = [[2, -1, 1],
[3, 3, 9],
[3, 3, 5]]
b = [2, -1, 2]