数值积分上机课写的 Python 代码,实现线性方程组的直接解法。
import numpy as np
import time
def gaussian_elimination(A, b):
if A.shape[0] != A.shape[1]:
raise ValueError('A must be a square matirx', A)
if b.shape[1] != 1:
raise ValueError('B must be a column vector', b)
if A.shape[0] != b.shape[0]:
raise ValueError('A and b must have the same number of rows', A, b)
n = A.shape[0]
if np.linalg.matrix_rank(A) != n:
raise ValueError('A is a singular matrix', A)
X = np.matrix(np.empty(n)).T
if A[0, 0] == 0:
for i in range(1, len(A)):
if A[i, 0] != 0:
A[[0, i], :] = A[[i, 0], :]
break
for j in range(n):
for i in range(j + 1, n):
if A[j, j] == 0:
for i in range(j + 1, len(A)):
if A[i, j] != 0:
A[[j, i], :] = A[[i, j], :]
break
m = A[i, j] / A[j, j]
A[i] = A[i] - A[j] * m
b[i] = b[i] - b[j] * m
X[-1] = b[-1] / A[-1, -1]
for i in reversed(range(n - 1)):
sum = 0
for j in range(i + 1, n):
sum += A[i, j] * X[j]
X[i] = (b[i] - sum) / A[i, i]
return X
def lu_decomposition(A, b):
def LU(A):
L = np.eye(len(A))
U = np.zeros(np.shape(A))
for r in range(1, len(A)):
U[0, r - 1] = A[0, r - 1]
L[r, 0] = A[r, 0] / A[0, 0]
U[0, -1] = A[0, -1]
for r in range(1, len(A)):
for i in range(r, len(A)):
delta = 0
for k in range(0, r):
delta += L[r, k] * U[k, i]
U[r, i] = A[r, i] - delta
for i in range(r + 1, len(A)):
theta = 0
for k in range(0, r):
theta += L[i, k] * U[k, r]
L[i, r] = (A[i, r] - theta) / U[r, r]
return L,U
if A.shape[0] != A.shape[1]:
raise ValueError('A must be a square matirx', A)
if b.shape[1] != 1:
raise ValueError('B must be a column vector', b)
if A.shape[0] != b.shape[0]:
raise ValueError('A and b must have the same number of rows', A, b)
n = A.shape[0]
if np.linalg.matrix_rank(A) != n:
raise ValueError('A is a singular matrix', A)
L, U = LU(A)
n = len(A)
y = np.zeros((n, 1))
b = np.array(b).reshape(n, 1)
for i in range(len(A)):
t = 0
for j in range(i):
t += L[i][j] * y[j][0]
y[i][0] = b[i][0] - t
X = np.zeros((n, 1))
for i in range(len(A) - 1, -1, -1):
t = 0
for j in range(i + 1, len(A)):
t += U[i][j] * X[j][0]
t = y[i][0] - t
if t != 0 and U[i][i] == 0:
return 0
X[i] = t / U[i][i]
return X