In [1]:
# ANS-1
import numpy as np
def solve_linear_system_custom(A_in, b_in):
"""
使用不带行交换的高斯消元法求解 Ax=b，并返回 LU 分解。
这种方法严格对应手动计算的步骤。
"""
# 转换为 float 类型以避免整数除法截断
A = A_in.copy().astype(float)
b = b_in.copy().astype(float)
n = len(b)
# 初始化 L 和 U
# L 初始化为单位矩阵
L = np.eye(n)
# --- 1. 高斯消元过程 (Gaussian Elimination) ---
# 目标：将 A 转换为上三角矩阵 U，并记录乘数到 L
# 对应书中 Algorithm 12
for k in range(n-1): # 遍历每一列 k (作为主元列)
pivot = A[k, k]
if abs(pivot) < 1e-12:
raise ValueError("主元为 0，需要行交换（本代码仅演示不交换的情况）")
for i in range(k+1, n): # 遍历主元下方的每一行 i
# 计算乘数 m
multiplier = A[i, k] / pivot
# 将乘数存入 L 矩阵的对应位置
L[i, k] = multiplier
# 更新第 i 行：Row_i = Row_i - m * Row_k
A[i, k:] = A[i, k:] - multiplier * A[k, k:]
# 同步更新右端项 b (这就相当于把 b 放在增广矩阵里一起变换)
# b[i] = b[i] - m * b[k]
# 注意：实际上我们在 LU 分解思路中通常分步做，
# 但为了直观展示手动过程，我们可以把 b 视为增广列的一部分进行变换
b[i] = b[i] - multiplier * b[k]
# 此时，A 已经被变换成了上三角矩阵 U
U = A
# --- 2. 回代求解 (Back Substitution) ---
# 对应书中 Example 9 [cite: 78-83]
x = np.zeros(n)
for i in range(n-1, -1, -1):
# x[i] = (b[i] - sum(U[i,j]*x[j])) / U[i,i]
sum_ax = sum(U[i, j] * x[j] for j in range(i+1, n))
x[i] = (b[i] - sum_ax) / U[i, i]
return x, L, U
# --- 主程序 ---
print("--- 问题 1 解答 (求解) ---")
# 定义题目中的矩阵 A 和向量 b
A_input = np.array([
[1, 1, 2],
[4, 3, 0],
[0, 2, 2]
])
b_input = np.array([2, 5, 6])
print("--- 原始输入 ---")
print(f"A =\n{A_input}")
print(f"b = {b_input}\n")
# 调用求解函数
x_sol, L_mat, U_mat = solve_linear_system_custom(A_input, b_input)
print("--- 计算结果 ---")
print(f"解向量 x = {x_sol}")
print("(这与我们手工计算的 [-1, 3, 0] 完全一致)\n")
print("--- LU 分解结果 ---")
print(f"L (下三角矩阵) =\n{L_mat}")
print(f"U (上三角矩阵) =\n{U_mat}")
print("\n--- 验证 ---")
print(f"验证 L * U 是否等于 A:\n{np.dot(L_mat, U_mat)}")
print(f"验证 A * x 是否等于 b:\n{np.dot(A_input, x_sol)}")
print("\n" + "="*50 + "\n")
--- 问题 1 解答 (求解) --- --- 原始输入 --- A = [[1 1 2] [4 3 0] [0 2 2]] b = [2 5 6] --- 计算结果 --- 解向量 x = [-1. 3. -0.] (这与我们手工计算的 [-1, 3, 0] 完全一致) --- LU 分解结果 --- L (下三角矩阵) = [[ 1. 0. 0.] [ 4. 1. 0.] [ 0. -2. 1.]] U (上三角矩阵) = [[ 1. 1. 2.] [ 0. -1. -8.] [ 0. 0. -14.]] --- 验证 --- 验证 L * U 是否等于 A: [[1. 1. 2.] [4. 3. 0.] [0. 2. 2.]] 验证 A * x 是否等于 b: [2. 5. 6.] ==================================================
In [2]:
# ANS-2
import numpy as np
def gaussian_elimination_pivot(A_in, b_in, tol=1e-10):
"""
带列主元选择的高斯消元法。
如果主元绝对值小于 tol，则在下方寻找最大元素进行行交换。
"""
# 转换为 float 类型，防止整数除法
A = A_in.copy().astype(float)
b = b_in.copy().astype(float)
n = len(b)
# --- 1. 消元过程 (Forward Elimination) ---
for k in range(n-1):
# === 列主元选择 (Pivoting) ===
# 检查当前主元 A[k,k] 是否过小
if abs(A[k, k]) < tol:
# 在当前列 k 的下方 (k+1 到 n) 寻找绝对值最大的元素
# np.argmax 返回的是相对于切片的索引，所以要加上 k
max_index_rel = np.argmax(np.abs(A[k:, k]))
max_index = max_index_rel + k
# 如果找到的最大值还是很小，说明矩阵奇异
if abs(A[max_index, k]) < tol:
raise ValueError("矩阵是奇异的 (Singular Matrix)")
# 交换行 (Swap rows)
# 注意：要在 A 和 b 上同时进行交换
A[[k, max_index]] = A[[max_index, k]]
b[[k, max_index]] = b[[max_index, k]]
# print(f" [Log] Step {k}: Swapped row {k} with row {max_index}")
# === 标准消元步骤 ===
for i in range(k+1, n):
multiplier = A[i, k] / A[k, k]
A[i, k:] = A[i, k:] - multiplier * A[k, k:]
b[i] = b[i] - multiplier * b[k]
# --- 2. 回代过程 (Back Substitution) ---
x = np.zeros(n)
for i in range(n-1, -1, -1):
sum_ax = sum(A[i, j] * x[j] for j in range(i+1, n))
x[i] = (b[i] - sum_ax) / A[i, i]
return x
# === 题目数据 ===
A = np.array([
[0, 1, 4],
[2, 4, 6],
[5, 6, 0]
])
b = np.array([8, 15, 5])
print("--- 问题 2 解答 ---")
# 1. 使用我们手写的带 Pivoting 的高斯消元
try:
x_custom = gaussian_elimination_pivot(A, b)
print(f"自定义代码结果 (x): {x_custom}")
except Exception as e:
print(f"出错: {e}")
# 2. 使用 Python 内置的 'backslash' (numpy.linalg.solve)
# numpy 的 solve 内部使用 LAPACK 库，默认带有 pivoting
x_numpy = np.linalg.solve(A, b)
print(f"Numpy solve 结果 : {x_numpy}")
# 验证两者是否一致
if np.allclose(x_custom, x_numpy):
print("验证成功：自定义代码与内置函数结果一致。")
else:
print("验证失败：结果不一致。")
print("\n" + "="*50 + "\n")
--- 问题 2 解答 --- 自定义代码结果 (x): [-11. 10. -0.5] Numpy solve 结果 : [-11. 10. -0.5] 验证成功：自定义代码与内置函数结果一致。 ==================================================
In [3]:
# ANS-3
import time
import scipy.sparse as sparse
import scipy.sparse.linalg as splinalg
print("--- 问题 3 解答 (性能测试) ---")
# 定义测试的矩阵大小 n
ns = [10, 100, 1000, 2000, 4000, 8000]
print(f"{'N':<10} | {'Dense Time (s)':<15} | {'Sparse Time (s)':<15} | {'Speedup (x)':<10}")
print("-" * 60)
for n in ns:
# 构造右端项 b (全为 1)
b = np.ones(n)
# === (a) 稠密矩阵 (Dense Matrix) ===
# 使用 np.diag 构建矩阵
# 主对角线 -2，上方偏移1为 1，下方偏移-1为 1
main_diag = -2 * np.ones(n)
upper_diag = np.ones(n-1)
lower_diag = np.ones(n-1)
A_dense = np.diag(main_diag, 0) + np.diag(upper_diag,
1) + np.diag(lower_diag, -1)
# 计时：稠密求解
start_dense = time.time()
x_dense = np.linalg.solve(A_dense, b)
end_dense = time.time()
time_dense = end_dense - start_dense
# === (b) 稀疏矩阵 (Sparse Matrix) ===
# 使用 scipy.sparse.diags 构建 (类似 MATLAB 的 spdiags)
# data 存放对角线数据，offsets 存放对应的偏移量
data = [main_diag, upper_diag, lower_diag]
offsets = [0, 1, -1]
# format='csr' (Compressed Sparse Row) 是求解线性方程组常用的格式
A_sparse = sparse.diags(data, offsets, shape=(n, n), format='csr')
# 计时：稀疏求解
start_sparse = time.time()
# 使用 scipy 的稀疏求解器
x_sparse = splinalg.spsolve(A_sparse, b)
end_sparse = time.time()
time_sparse = end_sparse - start_sparse
# 计算加速比
speedup = time_dense / time_sparse if time_sparse > 0 else 0
print(f"{n:<10} | {time_dense:<15.6f} | {time_sparse:<15.6f} | {speedup:<10.1f}")
print("\n结论：随着 n 的增加，稠密矩阵求解时间呈指数级增长 (约 O(n^3))，")
print("而稀疏矩阵求解非常快，因为通过忽略零元素，它避免了大量的无效计算。")
print("对于 n=8000，加速比通常会达到几百甚至上千倍。")
print("\n" + "="*50 + "\n")
--- 问题 3 解答 (性能测试) --- N | Dense Time (s) | Sparse Time (s) | Speedup (x) ------------------------------------------------------------ 10 | 0.000080 | 0.000186 | 0.4 100 | 0.000153 | 0.000150 | 1.0 1000 | 0.021276 | 0.000856 | 24.8 2000 | 0.116385 | 0.000460 | 253.1 4000 | 0.886252 | 0.000883 | 1003.8 8000 | 6.578838 | 0.001818 | 3618.8 结论：随着 n 的增加，稠密矩阵求解时间呈指数级增长 (约 O(n^3))， 而稀疏矩阵求解非常快，因为通过忽略零元素，它避免了大量的无效计算。 对于 n=8000，加速比通常会达到几百甚至上千倍。 ==================================================
In [4]:
# ANS-4
print("--- 问题 4 解答 (化学配平) ---")
# 已知 y
y = 100.0
# 构建线性方程组 Ax = b
# 未知数向量为 [x, z, w]
# 行 1: C原子平衡 -> 7x - 7z = 0
# 行 2: N原子平衡 -> 3z = y
# 行 3: O原子平衡 -> 6z + w = 3y
A_chem = np.array([
[7, -7, 0], # x, z, w 的系数
[0, 3, 0],
[0, 6, 1]
])
b_chem = np.array([
0, # RHS for C
y, # RHS for N (1*y)
3 * y # RHS for O (3*y)
])
# 求解
try:
solution = np.linalg.solve(A_chem, b_chem)
x_sol, z_sol, w_sol = solution
print(f"当输入 Nitric Acid (y) = {y} 时，计算结果如下：")
print(f"x (Toluene) = {x_sol:.4f}")
print(f"z (TNT) = {z_sol:.4f}")
print(f"w (Water) = {w_sol:.4f}")
# --- 验证步骤 ---
# 我们使用了 C, N, O 的方程求解。现在用 H 的方程来验证结果是否一致。
# H方程: 8x + y = 5z + 2w
# 即: LHS = 8x + y, RHS = 5z + 2w
lhs_h = 8 * x_sol + y
rhs_h = 5 * z_sol + 2 * w_sol
print(f"\n验证 H (氢) 原子守恒:")
print(f"反应前 H 数量: {lhs_h:.4f}")
print(f"反应后 H 数量: {rhs_h:.4f}")
if np.isclose(lhs_h, rhs_h):
print("验证通过！方程组是相容的。")
else:
print("验证失败！")
except Exception as e:
print(f"求解失败: {e}")
--- 问题 4 解答 (化学配平) --- 当输入 Nitric Acid (y) = 100.0 时，计算结果如下： x (Toluene) = 33.3333 z (TNT) = 33.3333 w (Water) = 100.0000 验证 H (氢) 原子守恒: 反应前 H 数量: 366.6667 反应后 H 数量: 366.6667 验证通过！方程组是相容的。