Agent 22 - Python科学计算实战（第一性原理版）
费曼原则：如果你不能用代码实现它，你就没有真正理解它。
核心哲学：代码即证明
数学概念往往是抽象的，但Python让抽象变得可触摸、可验证、可可视化。本教程遵循"第一性原理"思维：
- 用代码验证数学：不写数学伪代码，写真实可运行的Python
- 让抽象变具体：每个概念都有可视化呈现
- 从数值到符号：从数值实验过渡到严格数学证明
第一层：Python科研环境
1.1 环境搭建
Anaconda安装与配置
# 安装Anaconda后，创建隔离的数学计算环境
conda create -n math python=3.11
conda activate math
# 安装核心科学计算库
conda install numpy scipy sympy matplotlib pandas jupyter
pip install plotly ipywidgets
Jupyter Notebook/Lab配置
# 在~/.jupyter/jupyter_notebook_config.py中添加
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from scipy import *
# 自动导入设置
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
VS Code科研配置
// .vscode/settings.json
{
"python.analysis.typeCheckingMode": "basic",
"jupyter.askForKernelRestart": false,
"editor.formatOnSave": true,
"python.analysis.autoImportCompletions": true
}
1.2 核心库速览
| 库 | 功能 | 典型应用 |
|---|---|---|
| NumPy | 数组与矩阵运算 | 向量化计算、线性代数 |
| SciPy | 科学计算 | 积分、优化、微分方程 |
| SymPy | 符号计算 | 解析解、公式推导 |
| Matplotlib | 可视化 | 函数图像、数据分析 |
| Plotly | 交互式可视化 | 3D图形、动态演示 |
NumPy基础（向量化编程入门）
import numpy as np
# 创建数组
x = np.linspace(0, 2*np.pi, 1000) # 1000个点均匀分布在[0, 2π]
y = np.sin(x) # 向量化计算，无需循环
# 矩阵运算
A = np.array([[1, 2], [3, 4]])
B = np.linalg.inv(A) # 矩阵求逆
eigenvalues = np.linalg.eigvals(A) # 特征值
# 性能对比：向量化 vs 循环
def compute_loop(n):
"""Python循环版本（慢）"""
result = []
for i in range(n):
result.append(np.sin(i/n))
return result
def compute_vectorized(n):
"""NumPy向量化版本（快）"""
x = np.arange(n) / n
return np.sin(x)
# 向量化通常快100倍以上
第二层：高数概念验证
2.1 极限与连续
2.1.1 数列收敛可视化
费曼视角：极限不是一个遥远的目标，而是一个越来越近的过程。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
def visualize_sequence_convergence():
"""
可视化数列 a_n = 1/n 收敛到 0
费曼解释：想象一只蚂蚁每次走一半，它最终会无限接近但永不触碰
"""
fig, ax = plt.subplots(figsize=(10, 6))
# 生成数列 a_n = 1/n
n = np.arange(1, 101)
a_n = 1 / n
# 绘制收敛过程
ax.scatter(n, a_n, c=n, cmap='viridis', s=50, alpha=0.7)
ax.axhline(y=0, color='red', linestyle='--', label='极限 L=0')
ax.set_xlabel('n', fontsize=12)
ax.set_ylabel('$a_n = \\frac{1}{n}$', fontsize=12)
ax.set_title('数列收敛可视化：当 n→∞ 时，a_n→0', fontsize=14)
ax.legend()
ax.grid(True, alpha=0.3)
# 添加收敛带
epsilon = 0.05
ax.axhspan(-epsilon, epsilon, alpha=0.2, color='green', label=f'ε={epsilon}')
N_epsilon = int(np.ceil(1/epsilon))
ax.axvline(x=N_epsilon, color='green', linestyle=':', alpha=0.7)
ax.text(N_epsilon+2, 0.5, f'N={N_epsilon}', fontsize=10)
return fig
# 更多收敛数列示例
def demonstrate_convergence_types():
"""
展示不同类型的收敛
"""
n = np.arange(1, 51)
sequences = {
'1/n (单调收敛)': 1/n,
'(-1)^n/n (振荡收敛)': (-1)**n / n,
'1 - 1/n (单调递增)': 1 - 1/n,
'n/(n+1) (有理函数)': n/(n+1),
'(1+1/n)^n → e': (1 + 1/n)**n
}
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
for idx, (name, seq) in enumerate(sequences.items()):
ax = axes[idx]
ax.plot(n, seq, 'bo-', markersize=4, alpha=0.7)
# 标注极限
limit = seq[-1]
ax.axhline(y=limit, color='red', linestyle='--', alpha=0.7, label=f'极限≈{limit:.4f}')
ax.set_title(name)
ax.set_xlabel('n')
ax.set_ylabel('$a_n$')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
return fig
2.1.2 ε-δ定义验证
def verify_epsilon_delta():
"""
验证 lim(x→2) x² = 4 的 ε-δ 定义
对于任意 ε > 0，存在 δ > 0，使得当 0 < |x-2| < δ 时，|x²-4| < ε
费曼视角：无论多窄的"误差带"（ε），我总能找到足够小的"邻域"（δ）
"""
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左图：直观理解
ax1 = axes[0]
x = np.linspace(0, 4, 1000)
y = x**2
epsilon = 0.5
delta = epsilon / 5 # 因为导数在x=2处为4，所以 δ ≈ ε/4
ax1.plot(x, y, 'b-', linewidth=2, label='$f(x) = x^2$')
ax1.axhline(y=4, color='gray', linestyle='-', alpha=0.5)
ax1.axvline(x=2, color='gray', linestyle='-', alpha=0.5)
# ε-带
ax1.axhspan(4-epsilon, 4+epsilon, alpha=0.2, color='red', label=f'ε={epsilon}')
# δ-带
ax1.axvspan(2-delta, 2+delta, alpha=0.2, color='green', label=f'δ={delta:.3f}')
ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.set_title('ε-δ 定义的直观理解')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, 4])
ax1.set_ylim([0, 16])
# 右图：验证函数
ax2 = axes[1]
epsilons = np.logspace(-3, 0, 50) # 从0.001到1
# 精确计算：|x²-4| < ε ⟹ |x-2||x+2| < ε
# 在 x≈2 附近，|x+2|≈4，所以 δ = ε/4
deltas = epsilons / 4
ax2.loglog(epsilons, deltas, 'b-', linewidth=2, label='δ = ε/4')
ax2.loglog(epsilons, epsilons, 'r--', alpha=0.5, label='δ = ε (参考线)')
ax2.set_xlabel('ε')
ax2.set_ylabel('δ')
ax2.set_title('ε-δ 对应关系')
ax2.legend()
ax2.grid(True, alpha=0.3)
return fig
def find_delta_for_epsilon(f, c, L, epsilon, x_range=(0, 5)):
"""
数值方法：给定ε，寻找合适的δ
Args:
f: 目标函数
c: 趋近点
L: 极限值
epsilon: 误差限
x_range: 搜索范围
Returns:
delta: 找到的δ值
"""
x = np.linspace(x_range[0], x_range[1], 10000)
y = f(x)
# 找到满足 |f(x)-L| < ε 的 x 区间
mask = np.abs(y - L) < epsilon
valid_x = x[mask]
if len(valid_x) == 0:
return None
# 找到包含 c 的最大连通区间
left = np.min(valid_x[valid_x < c]) if np.any(valid_x < c) else c
right = np.max(valid_x[valid_x > c]) if np.any(valid_x > c) else c
delta = min(c - left, right - c)
return delta
# 使用示例
# epsilon = 0.01
# delta = find_delta_for_epsilon(lambda x: x**2, 2, 4, epsilon)
# print(f"当 ε={epsilon} 时，取 δ={delta:.6f} 可保证 |x²-4| < {epsilon}")
2.1.3 函数连续性检测
def check_continuity(f, x, tol=1e-6):
"""
数值检测函数在点x处的连续性
原理：lim(h→0) f(x+h) = f(x)
"""
h_values = np.logspace(-10, -1, 100)
left_limits = [f(x - h) for h in h_values]
right_limits = [f(x + h) for h in h_values]
f_x = f(x)
left_continuous = np.allclose(left_limits[::-1], f_x, atol=tol)
right_continuous = np.allclose(right_limits, f_x, atol=tol)
return {
'left_continuous': left_continuous,
'right_continuous': right_continuous,
'continuous': left_continuous and right_continuous,
'left_limit': left_limits[-1] if left_limits else None,
'right_limit': right_limits[-1] if right_limits else None,
'function_value': f_x
}
# 示例：检测各类不连续点
def visualize_discontinuities():
"""
可视化不同类型的不连续性
"""
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
x = np.linspace(-2, 2, 1000)
# 1. 可去间断点
ax1 = axes[0, 0]
x1 = x[x != 1]
y1 = (x1**2 - 1) / (x1 - 1) # = x+1 当 x≠1
ax1.plot(x1, y1, 'b-', linewidth=2)
ax1.plot(1, 2, 'ro', markersize=10, fillstyle='none', markeredgewidth=2, label='可去间断点')
ax1.set_title('可去间断点: $f(x) = \\frac{x^2-1}{x-1}$')
ax1.grid(True, alpha=0.3)
ax1.legend()
# 2. 跳跃间断点
ax2 = axes[0, 1]
y2 = np.where(x < 0, -1, 1)
ax2.plot(x[x < 0], y2[x < 0], 'b-', linewidth=2)
ax2.plot(x[x >= 0], y2[x >= 0], 'b-', linewidth=2)
ax2.plot(0, 0, 'ro', markersize=10, fillstyle='none', markeredgewidth=2)
ax2.set_title('跳跃间断点: 符号函数')
ax2.set_ylim([-2, 2])
ax2.grid(True, alpha=0.3)
# 3. 无穷间断点
ax3 = axes[1, 0]
x3 = x[x != 0]
y3 = 1 / x3
ax3.plot(x3[x3 < 0], y3[x3 < 0], 'b-', linewidth=2)
ax3.plot(x3[x3 > 0], y3[x3 > 0], 'b-', linewidth=2)
ax3.axvline(x=0, color='red', linestyle='--', label='无穷间断点')
ax3.set_title('无穷间断点: $f(x) = \\frac{1}{x}$')
ax3.set_ylim([-10, 10])
ax3.grid(True, alpha=0.3)
ax3.legend()
# 4. 振荡间断点
ax4 = axes[1, 1]
x4 = x[x != 0]
y4 = np.sin(1/x4)
ax4.plot(x4, y4, 'b-', linewidth=0.5)
ax4.axvline(x=0, color='red', linestyle='--', label='振荡间断点')
ax4.set_title('振荡间断点: $f(x) = \\sin(\\frac{1}{x})$')
ax4.set_ylim([-1.5, 1.5])
ax4.grid(True, alpha=0.3)
ax4.legend()
plt.tight_layout()
return fig
2.2 导数与微分
2.2.1 数值微分方法
import numpy as np
class NumericalDifferentiation:
"""
数值微分：当解析解难求时，用有限差分近似导数
费曼理解：导数是"变化率的放大镜"
"""
@staticmethod
def forward_diff(f, x, h=1e-5):
"""前向差分: f'(x) ≈ [f(x+h) - f(x)] / h
误差: O(h)
"""
return (f(x + h) - f(x)) / h
@staticmethod
def backward_diff(f, x, h=1e-5):
"""后向差分: f'(x) ≈ [f(x) - f(x-h)] / h
误差: O(h)
"""
return (f(x) - f(x - h)) / h
@staticmethod
def central_diff(f, x, h=1e-5):
"""中心差分: f'(x) ≈ [f(x+h) - f(x-h)] / (2h)
误差: O(h²) —— 更精确！
"""
return (f(x + h) - f(x - h)) / (2 * h)
@staticmethod
def five_point_stencil(f, x, h=1e-5):
"""五点公式: f'(x) ≈ [-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)] / (12h)
误差: O(h⁴) —— 高精度
"""
return (-f(x + 2*h) + 8*f(x + h) - 8*f(x - h) + f(x - 2*h)) / (12 * h)
@staticmethod
def second_derivative(f, x, h=1e-5):
"""二阶导数: f''(x) ≈ [f(x+h) - 2f(x) + f(x-h)] / h²
物理意义：加速度、曲率
"""
return (f(x + h) - 2*f(x) + f(x - h)) / (h**2)
# 精度验证示例
def verify_numerical_differentiation():
"""
用已知导数的函数验证数值微分的精度
"""
# 测试函数：f(x) = e^x * sin(x)
# 精确导数：f'(x) = e^x * (sin(x) + cos(x))
def f(x):
return np.exp(x) * np.sin(x)
def df_exact(x):
return np.exp(x) * (np.sin(x) + np.cos(x))
x_test = 1.0
exact = df_exact(x_test)
# 测试不同步长
h_values = np.logspace(-10, -1, 50)
methods = {
'Forward': NumericalDifferentiation.forward_diff,
'Backward': NumericalDifferentiation.backward_diff,
'Central': NumericalDifferentiation.central_diff,
'Five-point': NumericalDifferentiation.five_point_stencil
}
errors = {name: [] for name in methods}
for h in h_values:
for name, method in methods.items():
approx = method(f, x_test, h)
error = abs(approx - exact)
errors[name].append(error)
# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
for name, errs in errors.items():
ax.loglog(h_values, errs, 'o-', label=name, markersize=3)
ax.axhline(y=np.finfo(float).eps, color='k', linestyle='--', alpha=0.5, label='Machine ε')
ax.set_xlabel('Step size h')
ax.set_ylabel('Absolute Error')
ax.set_title('数值微分方法的误差比较')
ax.legend()
ax.grid(True, alpha=0.3)
return fig
2.2.2 梯度计算与自动微分
import numpy as np
def gradient_2d(f, x, y, h=1e-5):
"""
计算二元函数的梯度 ∇f = [∂f/∂x, ∂f/∂y]
费曼几何解释：梯度指向函数增长最快的方向
"""
df_dx = (f(x + h, y) - f(x - h, y)) / (2 * h)
df_dy = (f(x, y + h) - f(x, y - h)) / (2 * h)
return np.array([df_dx, df_dy])
def visualize_gradient_field():
"""
可视化梯度场：
- 等高线表示函数值
- 箭头表示梯度方向（垂直于等高线，指向高处）
"""
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 创建网格
x = np.linspace(-3, 3, 50)
y = np.linspace(-3, 3, 50)
X, Y = np.meshgrid(x, y)
# 测试函数1：f(x,y) = x² + y² (碗形)
Z1 = X**2 + Y**2
ax1 = axes[0]
contour1 = ax1.contour(X, Y, Z1, levels=15, cmap='viridis')
ax1.clabel(contour1, inline=True, fontsize=8)
# 计算并绘制梯度场
skip = 3
dZ1_dx = 2 * X
dZ1_dy = 2 * Y
ax1.quiver(X[::skip, ::skip], Y[::skip, ::skip],
dZ1_dx[::skip, ::skip], dZ1_dy[::skip, ::skip],
scale=50, color='red', alpha=0.7)
ax1.set_title('$f(x,y) = x^2 + y^2$ 的梯度场')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_aspect('equal')
# 测试函数2：f(x,y) = sin(x) * cos(y) (波纹形)
Z2 = np.sin(X) * np.cos(Y)
ax2 = axes[1]
contour2 = ax2.contour(X, Y, Z2, levels=15, cmap='viridis')
ax2.clabel(contour2, inline=True, fontsize=8)
dZ2_dx = np.cos(X) * np.cos(Y)
dZ2_dy = -np.sin(X) * np.sin(Y)
ax2.quiver(X[::skip, ::skip], Y[::skip, ::skip],
dZ2_dx[::skip, ::skip], dZ2_dy[::skip, ::skip],
scale=20, color='red', alpha=0.7)
ax2.set_title('$f(x,y) = \\sin(x)\\cos(y)$ 的梯度场')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_aspect('equal')
plt.tight_layout()
return fig
2.2.3 泰勒展开验证
import numpy as np
import sympy as sp
def taylor_series_visualization():
"""
可视化泰勒展开：用多项式逼近任意函数
费曼洞察：泰勒展开就像用"本地信息"（函数值和各阶导数）
重建"全局行为"
"""
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
x = np.linspace(-np.pi, np.pi, 500)
# 1. sin(x) 在 x=0 处的泰勒展开
ax1 = axes[0, 0]
ax1.plot(x, np.sin(x), 'k-', linewidth=2, label='$\\sin(x)$')
colors = plt.cm.viridis(np.linspace(0, 1, 5))
for n, color in zip([1, 3, 5, 7, 9], colors):
# 泰勒展开: sin(x) ≈ Σ (-1)^k * x^(2k+1) / (2k+1)!
taylor = sum([(-1)**k * x**(2*k+1) / np.math.factorial(2*k+1)
for k in range((n+1)//2)])
ax1.plot(x, taylor, '--', color=color, alpha=0.8,
label=f'n={n}')
ax1.set_title('$\\sin(x)$ 的泰勒展开（在 x=0 处）')
ax1.set_ylim([-2, 2])
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)
# 2. e^x 的泰勒展开
ax2 = axes[0, 1]
ax2.plot(x, np.exp(x), 'k-', linewidth=2, label='$e^x$')
for n, color in zip([1, 2, 3, 4, 5], colors):
taylor = sum([x**k / np.math.factorial(k) for k in range(n+1)])
ax2.plot(x, taylor, '--', color=color, alpha=0.8,
label=f'n={n}')
ax2.set_title('$e^x$ 的泰勒展开')
ax2.set_ylim([-2, 10])
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)
# 3. 收敛半径可视化（以 1/(1-x) 为例）
ax3 = axes[1, 0]
x3 = np.linspace(-2, 2, 500)
x3_valid = x3[np.abs(x3) < 1]
ax3.plot(x3_valid, 1/(1-x3_valid), 'k-', linewidth=2, label='$\\frac{1}{1-x}$')
for n, color in zip([3, 5, 10, 20], colors):
taylor = sum([x3**k for k in range(n+1)]) # 几何级数
ax3.plot(x3, taylor, '--', color=color, alpha=0.8,
label=f'n={n}')
ax3.axvline(x=1, color='red', linestyle=':', alpha=0.7, label='收敛边界 |x|=1')
ax3.axvline(x=-1, color='red', linestyle=':', alpha=0.7)
ax3.set_title('收敛半径示例：$\\frac{1}{1-x}$')
ax3.set_ylim([-5, 10])
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.axhline(y=0, color='k', linewidth=0.5)
ax3.axvline(x=0, color='k', linewidth=0.5)
# 4. 误差分析
ax4 = axes[1, 1]
x0 = 1.0 # 在 x=1 处评估 sin(x) 的泰勒展开误差
n_values = range(1, 20)
exact = np.sin(x0)
errors = []
for n in n_values:
approx = sum([(-1)**k * x0**(2*k+1) / np.math.factorial(2*k+1)
for k in range((n+1)//2)])
errors.append(abs(approx - exact))
ax4.semilogy(n_values, errors, 'bo-', markersize=6)
ax4.set_xlabel('泰勒展开的阶数 n')
ax4.set_ylabel('绝对误差')
ax4.set_title('泰勒展开的收敛速度')
ax4.grid(True, alpha=0.3)
ax4.axhline(y=np.finfo(float).eps, color='red', linestyle='--',
label='机器精度')
ax4.legend()
plt.tight_layout()
return fig
# 使用 SymPy 进行符号泰勒展开
def symbolic_taylor_expansion():
"""
使用 SymPy 进行精确的符号泰勒展开
"""
x = sp.Symbol('x')
functions = {
'sin(x)': sp.sin(x),
'cos(x)': sp.cos(x),
'e^x': sp.exp(x),
'ln(1+x)': sp.log(1+x)
}
results = {}
for name, func in functions.items():
# 在 x=0 处展开到 6 阶
taylor = sp.series(func, x, 0, 7).removeO()
results[name] = taylor
print(f"{name} ≈ {taylor}")
return results
2.3 积分计算
2.3.1 数值积分方法
import numpy as np
from scipy import integrate
class NumericalIntegration:
"""
数值积分：当原函数难求时，用数值方法近似定积分
费曼理解：积分是"无穷小累加的艺术"
"""
@staticmethod
def rectangle_rule(f, a, b, n=1000):
"""
矩形法（左端点）：∫[a,b]f(x)dx ≈ h * Σ f(x_i)
误差: O(1/n)
"""
h = (b - a) / n
x = np.linspace(a, b - h, n) # 左端点
return h * np.sum(f(x))
@staticmethod
def trapezoidal_rule(f, a, b, n=1000):
"""
梯形法：∫[a,b]f(x)dx ≈ h * [f(a)/2 + Σ f(x_i) + f(b)/2]
误差: O(1/n²)
"""
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = f(x)
return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])
@staticmethod
def simpson_rule(f, a, b, n=1000):
"""
Simpson法则（抛物线近似）：
∫[a,b]f(x)dx ≈ h/3 * [f(a) + 4Σf(x_odd) + 2Σf(x_even) + f(b)]
误差: O(1/n⁴) —— 高精度！
注意：n 必须是偶数
"""
if n % 2 == 1:
n += 1
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = f(x)
result = h/3 * (y[0] +
4 * np.sum(y[1:-1:2]) +
2 * np.sum(y[2:-2:2]) +
y[-1])
return result
@staticmethod
def gaussian_quadrature(f, a, b, n=5):
"""
高斯-勒让德求积：使用最优节点和权重
对于多项式，2n-1 阶以下精确
"""
# 使用 scipy 的实现
result, error = integrate.quad(f, a, b)
return result
# 误差比较实验
def compare_integration_methods():
"""
比较不同数值积分方法的精度和收敛速度
"""
# 测试函数：∫[0,π] sin(x) dx = 2
def f(x):
return np.sin(x)
exact = 2.0
a, b = 0, np.pi
n_values = [10, 20, 50, 100, 200, 500, 1000, 2000]
methods = {
'Rectangle': NumericalIntegration.rectangle_rule,
'Trapezoidal': NumericalIntegration.trapezoidal_rule,
'Simpson': NumericalIntegration.simpson_rule
}
errors = {name: [] for name in methods}
for n in n_values:
for name, method in methods.items():
approx = method(f, a, b, n)
error = abs(approx - exact)
errors[name].append(error)
# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
for name, errs in errors.items():
ax.loglog(n_values, errs, 'o-', label=name, markersize=6)
# 参考线
n_array = np.array(n_values)
ax.loglog(n_array, 1/n_array, 'k--', alpha=0.5, label='O(1/n)')
ax.loglog(n_array, 1/n_array**2, 'k:', alpha=0.5, label='O(1/n²)')
ax.loglog(n_array, 1/n_array**4, 'k-.', alpha=0.5, label='O(1/n⁴)')
ax.set_xlabel('Number of subintervals n')
ax.set_ylabel('Absolute Error')
ax.set_title('数值积分方法的收敛速度比较')
ax.legend()
ax.grid(True, alpha=0.3)
return fig
2.3.2 符号积分验证
import sympy as sp
def symbolic_integration_examples():
"""
使用 SymPy 进行精确符号积分
验证数值结果的准确性
"""
x = sp.Symbol('x')
# 基本积分
integrals = [
('多项式', x**3 + 2*x**2 - x + 1),
('三角函数', sp.sin(x)**2),
('指数函数', sp.exp(x) * sp.sin(x)),
('有理函数', 1/(1 + x**2)),
('根式函数', sp.sqrt(1 - x**2)),
]
results = []
for name, func in integrals:
# 不定积分
indefinite = sp.integrate(func, x)
# 验证：求导后应恢复原函数
verified = sp.simplify(sp.diff(indefinite, x) - func) == 0
results.append({
'name': name,
'function': func,
'integral': indefinite,
'verified': verified
})
print(f"\n{name}:")
print(f" ∫ {func} dx = {indefinite}")
print(f" 验证通过: {verified}")
return results
def definite_integral_examples():
"""
定积分示例：计算面积、体积、质心等
"""
x = sp.Symbol('x')
examples = [
('圆面积', sp.sqrt(1 - x**2), -1, 1, sp.pi/2),
('高斯积分', sp.exp(-x**2), -sp.oo, sp.oo, sp.sqrt(sp.pi)),
('反常积分', 1/x**2, 1, sp.oo, 1),
]
for name, func, a, b, expected in examples:
result = sp.integrate(func, (x, a, b))
print(f"\n{name}: ∫[{a},{b}] {func} dx = {result}")
if expected is not None:
print(f" 期望值: {expected}")
print(f" 相等: {sp.simplify(result - expected) == 0}")
2.3.3 反常积分处理
import numpy as np
from scipy import integrate
def improper_integration():
"""
处理反常积分（无穷区间或无界函数）
费曼洞察：反常积分是"可驯服的无穷"
"""
# 1. 无穷区间积分
def infinite_interval_examples():
examples = [
('高斯积分', lambda x: np.exp(-x**2), 0, np.inf, np.sqrt(np.pi)/2),
('p-积分 (p>1)', lambda x: 1/x**2, 1, np.inf, 1),
('指数衰减', lambda x: np.exp(-x), 0, np.inf, 1),
]
results = []
for name, f, a, b, expected in examples:
result, error = integrate.quad(f, a, b)
results.append({
'name': name,
'numerical': result,
'expected': expected,
'error': error,
'diff': abs(result - expected)
})
return results
# 2. 瑕积分（无界函数）
def singular_integration():
examples = [
('瑕积分', lambda x: 1/np.sqrt(x), 0, 1, 2),
('对数积分', lambda x: -np.log(x), 0, 1, 1),
]
results = []
for name, f, a, b, expected in examples:
# scipy 可以自动处理某些奇点
result, error = integrate.quad(f, a, b)
results.append({
'name': name,
'numerical': result,
'expected': expected
})
return results
# 3. 柯西主值
def cauchy_principal_value():
"""
计算 1/x 在 -1 到 1 的柯西主值
P.V. ∫[-1,1] 1/x dx = 0（奇函数对称）
"""
# 使用 scipy 的特殊参数
result, error, info = integrate.quad(lambda x: 1/x, -1, 1,
limit=100,
weight='cauchy', wvar=0)
return result
return {
'infinite': infinite_interval_examples(),
'singular': singular_integration(),
'cauchy_pv': cauchy_principal_value()
}
第三层：实战项目
项目1：函数可视化系统
# project1_function_visualizer.py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
import sympy as sp
class FunctionVisualizer:
"""
交互式函数可视化系统
功能：
- 绘制函数图像
- 显示导数与原函数关系
- 泰勒展开动画逼近
"""
def __init__(self):
self.x_range = (-10, 10)
self.fig = None
self.axes = None
def interactive_plot(self, func_str='sin(x)'):
"""
创建交互式函数绘图界面
"""
x = np.linspace(self.x_range[0], self.x_range[1], 1000)
# 将字符串转换为函数
x_sym = sp.Symbol('x')
func_sym = sp.sympify(func_str)
f = sp.lambdify(x_sym, func_sym, 'numpy')
fig, ax = plt.subplots(figsize=(12, 6))
plt.subplots_adjust(bottom=0.25)
y = f(x)
line, = ax.plot(x, y, 'b-', linewidth=2, label=f'y = {func_str}')
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)
ax.grid(True, alpha=0.3)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('函数可视化')
ax.legend()
# 添加导数显示
df_sym = sp.diff(func_sym, x_sym)
df = sp.lambdify(x_sym, df_sym, 'numpy')
# 创建滑块控制显示范围
axcolor = 'lightgoldenrodyellow'
ax_xmin = plt.axes([0.15, 0.1, 0.3, 0.03], facecolor=axcolor)
ax_xmax = plt.axes([0.6, 0.1, 0.3, 0.03], facecolor=axcolor)
sxmin = Slider(ax_xmin, 'X min', -20, 0, valinit=self.x_range[0])
sxmax = Slider(ax_xmax, 'X max', 0, 20, valinit=self.x_range[1])
def update(val):
xmin = sxmin.val
xmax = sxmax.val
x_new = np.linspace(xmin, xmax, 1000)
y_new = f(x_new)
line.set_xdata(x_new)
line.set_ydata(y_new)
ax.set_xlim([xmin, xmax])
ax.relim()
ax.autoscale_view()
fig.canvas.draw_idle()
sxmin.on_changed(update)
sxmax.on_changed(update)
return fig
def derivative_relationship(self, func_str='x**3 - 3*x'):
"""
可视化导数与原函数的关系
- 函数上升时，导数为正
- 函数下降时，导数为负
- 极值点，导数为零
"""
x_sym = sp.Symbol('x')
func_sym = sp.sympify(func_str)
df_sym = sp.diff(func_sym, x_sym)
f = sp.lambdify(x_sym, func_sym, 'numpy')
df = sp.lambdify(x_sym, df_sym, 'numpy')
x = np.linspace(-3, 3, 500)
y = f(x)
dy = df(x)
fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True)
# 原函数
ax1 = axes[0]
ax1.plot(x, y, 'b-', linewidth=2, label=f'$f(x) = {sp.latex(func_sym)}$')
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_ylabel('f(x)')
ax1.set_title('原函数')
# 导数
ax2 = axes[1]
ax2.plot(x, dy, 'r-', linewidth=2, label=f"$f'(x) = {sp.latex(df_sym)}$")
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.fill_between(x, 0, dy, where=(dy > 0), alpha=0.3, color='green', label='f 上升')
ax2.fill_between(x, 0, dy, where=(dy < 0), alpha=0.3, color='red', label='f 下降')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_ylabel("f'(x)")
ax2.set_title('导数')
# 对应关系
ax3 = axes[2]
ax3.plot(x, y, 'b-', linewidth=2, alpha=0.7, label='f(x)')
# 标记导数为零的点（极值点）
zero_crossings = x[:-1][np.diff(np.sign(dy)) != 0]
for zc in zero_crossings:
ax3.axvline(x=zc, color='purple', linestyle='--', alpha=0.5)
ax3.plot(zc, f(zc), 'ro', markersize=10)
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_ylabel('f(x)')
ax3.set_xlabel('x')
ax3.set_title('极值点（导数为零）')
plt.tight_layout()
return fig
def taylor_animation(self, func_str='sin(x)', center=0, max_degree=20):
"""
泰勒展开逼近动画
"""
x_sym = sp.Symbol('x')
func_sym = sp.sympify(func_str)
x = np.linspace(-np.pi, np.pi, 500)
y_exact = sp.lambdify(x_sym, func_sym, 'numpy')(x)
fig, ax = plt.subplots(figsize=(12, 6))
line_exact, = ax.plot(x, y_exact, 'k-', linewidth=2, label='Exact')
line_taylor, = ax.plot([], [], 'r--', linewidth=2, label='Taylor')
ax.set_xlim([-np.pi, np.pi])
ax.set_ylim([-2, 2])
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)
ax.grid(True, alpha=0.3)
ax.legend()
ax.set_title(f'Taylor Series Approximation of {func_str} at x={center}')
# 预计算各阶泰勒展开
taylor_funcs = []
for n in range(1, max_degree + 1):
taylor_expr = sp.series(func_sym, x_sym, center, n+1).removeO()
taylor_func = sp.lambdify(x_sym, taylor_expr, 'numpy')
taylor_funcs.append(taylor_func)
def animate(frame):
y_taylor = taylor_funcs[frame](x)
line_taylor.set_data(x, y_taylor)
ax.set_title(f'Taylor Series (n={frame+1})')
return line_taylor,
from matplotlib.animation import FuncAnimation
anim = FuncAnimation(fig, animate, frames=max_degree, interval=500, blit=True)
return fig, anim
# 使用示例
if __name__ == '__main__':
viz = FunctionVisualizer()
viz.interactive_plot('sin(x)')
plt.show()
项目2：数值积分器
# project2_numerical_integrator.py
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import sympy as sp
class NumericalIntegrator:
"""
通用数值积分器
功能：
- 多种积分方法实现
- 误差分析和自适应步长
- 可视化积分过程
"""
METHODS = {
'rectangle': '矩形法',
'trapezoidal': '梯形法',
'simpson': 'Simpson法',
'romberg': 'Romberg加速',
'gaussian': '高斯求积'
}
def integrate(self, f, a, b, method='simpson', n=1000, **kwargs):
"""
通用积分接口
Args:
f: 被积函数
a, b: 积分区间
method: 积分方法
n: 子区间数
"""
if method == 'rectangle':
return self._rectangle(f, a, b, n)
elif method == 'trapezoidal':
return self._trapezoidal(f, a, b, n)
elif method == 'simpson':
return self._simpson(f, a, b, n)
elif method == 'romberg':
return self._romberg(f, a, b, **kwargs)
elif method == 'gaussian':
return self._gaussian(f, a, b, **kwargs)
else:
raise ValueError(f"Unknown method: {method}")
def _rectangle(self, f, a, b, n):
"""左矩形法"""
h = (b - a) / n
x = np.linspace(a, b - h, n)
return h * np.sum(f(x))
def _trapezoidal(self, f, a, b, n):
"""梯形法"""
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = f(x)
return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])
def _simpson(self, f, a, b, n):
"""Simpson 1/3 法则"""
if n % 2 == 1:
n += 1
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = f(x)
return h/3 * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-2:2]) + y[-1])
def _romberg(self, f, a, b, max_iter=10, tol=1e-10):
"""
Romberg 积分：Richardson 外推加速
原理：通过不同步长的结果外推，消除低阶误差项
"""
R = np.zeros((max_iter, max_iter))
# 初始梯形法
n = 1
h = b - a
R[0, 0] = h * (f(a) + f(b)) / 2
for i in range(1, max_iter):
# 细分步长
n *= 2
h /= 2
# 计算新的梯形估计
x_new = np.linspace(a + h, b - h, n // 2)
R[i, 0] = R[i-1, 0] / 2 + h * np.sum(f(x_new))
# Richardson 外推
for j in range(1, i + 1):
R[i, j] = (4**j * R[i, j-1] - R[i-1, j-1]) / (4**j - 1)
# 检查收敛
if abs(R[i, i] - R[i-1, i-1]) < tol:
return R[i, i]
return R[max_iter-1, max_iter-1]
def _gaussian(self, f, a, b, n=5):
"""高斯-勒让德求积"""
# 使用 scipy
result, _ = integrate.quad(f, a, b, limit=n)
return result
def adaptive_simpson(self, f, a, b, tol=1e-6, max_depth=20):
"""
自适应 Simpson 积分
只在误差大的区域细分，提高效率
"""
def adaptive_recursion(f, a, b, eps, whole, fa, fb, fm, depth):
m = (a + b) / 2
lm = (a + m) / 2
rm = (m + b) / 2
flm = f(lm)
frm = f(rm)
left = (fa + 4*flm + fm) * (m - a) / 6
right = (fm + 4*frm + fb) * (b - m) / 6
delta = left + right - whole
if depth >= max_depth or abs(delta) <= 15 * eps:
return left + right + delta / 15
return (adaptive_recursion(f, a, m, eps/2, left, fa, fm, flm, depth+1) +
adaptive_recursion(f, m, b, eps/2, right, fm, fb, frm, depth+1))
fa, fb = f(a), f(b)
m = (a + b) / 2
fm = f(m)
whole = (fa + 4*fm + fb) * (b - a) / 6
return adaptive_recursion(f, a, b, tol, whole, fa, fb, fm, 0)
def visualize_integration(self, f, a, b, n=10, method='trapezoidal'):
"""
可视化积分过程
"""
fig, ax = plt.subplots(figsize=(12, 6))
# 绘制函数
x_fine = np.linspace(a, b, 1000)
ax.plot(x_fine, f(x_fine), 'b-', linewidth=2, label='f(x)')
# 绘制积分近似
x_coarse = np.linspace(a, b, n + 1)
y_coarse = f(x_coarse)
if method == 'rectangle':
for i in range(n):
rect = plt.Rectangle((x_coarse[i], 0),
x_coarse[i+1] - x_coarse[i],
y_coarse[i],
fill=True, alpha=0.3, color='red',
edgecolor='red')
ax.add_patch(rect)
elif method == 'trapezoidal':
for i in range(n):
x_trap = [x_coarse[i], x_coarse[i], x_coarse[i+1], x_coarse[i+1]]
y_trap = [0, y_coarse[i], y_coarse[i+1], 0]
ax.fill(x_trap, y_trap, alpha=0.3, color='red', edgecolor='red')
elif method == 'simpson':
# Simpson 使用抛物线近似
for i in range(0, n, 2):
if i + 2 <= n:
x_seg = np.linspace(x_coarse[i], x_coarse[i+2], 50)
# 抛物线插值
coeffs = np.polyfit([x_coarse[i], x_coarse[i+1], x_coarse[i+2]],
[y_coarse[i], y_coarse[i+1], y_coarse[i+2]], 2)
y_parabola = np.polyval(coeffs, x_seg)
ax.fill_between(x_seg, 0, y_parabola, alpha=0.3, color='red')
ax.axhline(y=0, color='k', linewidth=0.5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'{self.METHODS[method]} 可视化 (n={n})')
ax.legend()
ax.grid(True, alpha=0.3)
return fig
def error_analysis(self, f, a, b, exact_value):
"""
不同方法的误差分析
"""
n_values = np.logspace(1, 4, 20, dtype=int)
results = {}
for method in ['rectangle', 'trapezoidal', 'simpson']:
errors = []
for n in n_values:
approx = self.integrate(f, a, b, method, n)
error = abs(approx - exact_value)
errors.append(error)
results[method] = errors
# 绘制误差图
fig, ax = plt.subplots(figsize=(10, 6))
for method, errors in results.items():
ax.loglog(n_values, errors, 'o-', label=self.METHODS[method], markersize=4)
# 参考线
ax.loglog(n_values, 1/n_values, 'k--', alpha=0.5, label='O(1/n)')
ax.loglog(n_values, 1/n_values**2, 'k:', alpha=0.5, label='O(1/n²)')
ax.loglog(n_values, 1/n_values**4, 'k-.', alpha=0.5, label='O(1/n⁴)')
ax.set_xlabel('Number of subintervals')
ax.set_ylabel('Absolute Error')
ax.set_title('Integration Methods Error Comparison')
ax.legend()
ax.grid(True, alpha=0.3)
return fig
# 使用示例
if __name__ == '__main__':
integrator = NumericalIntegrator()
# 测试函数：高斯函数积分
def gaussian(x):
return np.exp(-x**2)
exact = np.sqrt(np.pi) # ∫[-∞,∞] e^(-x²) dx = √π
# 比较不同方法
for method in ['rectangle', 'trapezoidal', 'simpson', 'romberg']:
result = integrator.integrate(gaussian, -5, 5, method, n=1000)
print(f"{method}: {result:.10f}, error: {abs(result - exact):.2e}")
项目3：微分方程求解器
# project3_differential_equations.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
class ODESolver:
"""
常微分方程数值求解器
功能：
- 实现经典数值方法（欧拉、龙格-库塔）
- 绘制相平面图
- 比较不同方法的精度和稳定性
"""
def __init__(self):
self.methods = {
'euler': self.euler_method,
'improved_euler': self.improved_euler,
'rk4': self.rk4_method,
'rk45': self.rk45_scipy
}
def euler_method(self, f, y0, t):
"""
显式欧拉方法（一阶）
y_{n+1} = y_n + h * f(t_n, y_n)
简单但不稳定，仅用于教学
"""
y = np.zeros((len(t), len(y0) if hasattr(y0, '__len__') else 1))
y[0] = y0
for i in range(len(t) - 1):
h = t[i+1] - t[i]
y[i+1] = y[i] + h * np.array(f(t[i], y[i]))
return y.squeeze()
def improved_euler(self, f, y0, t):
"""
改进欧拉方法（Heun方法，二阶）
k1 = f(t_n, y_n)
k2 = f(t_n + h, y_n + h*k1)
y_{n+1} = y_n + h/2 * (k1 + k2)
"""
y = np.zeros((len(t), len(y0) if hasattr(y0, '__len__') else 1))
y[0] = y0
for i in range(len(t) - 1):
h = t[i+1] - t[i]
k1 = np.array(f(t[i], y[i]))
k2 = np.array(f(t[i+1], y[i] + h*k1))
y[i+1] = y[i] + h/2 * (k1 + k2)
return y.squeeze()
def rk4_method(self, f, y0, t):
"""
经典龙格-库塔方法（四阶，RK4）
k1 = f(t_n, y_n)
k2 = f(t_n + h/2, y_n + h/2*k1)
k3 = f(t_n + h/2, y_n + h/2*k2)
k4 = f(t_n + h, y_n + h*k3)
y_{n+1} = y_n + h/6 * (k1 + 2k2 + 2k3 + k4)
精度高，是实际应用中最常用的方法
"""
y = np.zeros((len(t), len(y0) if hasattr(y0, '__len__') else 1))
y[0] = y0
for i in range(len(t) - 1):
h = t[i+1] - t[i]
k1 = np.array(f(t[i], y[i]))
k2 = np.array(f(t[i] + h/2, y[i] + h/2*k1))
k3 = np.array(f(t[i] + h/2, y[i] + h/2*k2))
k4 = np.array(f(t[i+1], y[i] + h*k3))
y[i+1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
return y.squeeze()
def rk45_scipy(self, f, y0, t_span, **kwargs):
"""
使用 scipy 的自适应 RK45 方法
"""
sol = solve_ivp(f, t_span, [y0] if np.isscalar(y0) else y0,
method='RK45', dense_output=True, **kwargs)
t_dense = np.linspace(t_span[0], t_span[1], 1000)
y_dense = sol.sol(t_dense)
return t_dense, y_dense.T
def solve(self, f, y0, t, method='rk4'):
"""
通用求解接口
"""
if method in self.methods:
return self.methods[method](f, y0, t)
else:
raise ValueError(f"Unknown method: {method}")
def harmonic_oscillator(self, omega=1.0, method='rk4'):
"""
简谐振子：y'' + ω²y = 0
转化为方程组：
y1' = y2
y2' = -ω²*y1
解析解：y(t) = A*cos(ωt) + B*sin(ωt)
"""
def f(t, y):
return [y[1], -omega**2 * y[0]]
y0 = [1.0, 0.0] # 初始位移1，初始速度0
t = np.linspace(0, 4*np.pi/omega, 1000)
solution = self.solve(f, y0, t, method)
# 解析解
exact = np.cos(omega * t)
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
# 时域图
ax1 = axes[0]
ax1.plot(t, solution[:, 0], 'b-', label='Numerical')
ax1.plot(t, exact, 'r--', alpha=0.7, label='Exact: cos(t)')
ax1.set_xlabel('t')
ax1.set_ylabel('y')
ax1.set_title('Harmonic Oscillator Solution')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 相平面图（速度-位移）
ax2 = axes[1]
ax2.plot(solution[:, 0], solution[:, 1], 'b-', linewidth=1)
ax2.set_xlabel('Position (y)')
ax2.set_ylabel('Velocity (y\')')
ax2.set_title('Phase Portrait')
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')
plt.tight_layout()
return fig
def damped_pendulum(self, gamma=0.1, omega0=1.0):
"""
阻尼单摆：θ'' + γθ' + ω₀²sin(θ) = 0
"""
def f(t, y):
theta, theta_dot = y
return [theta_dot, -gamma*theta_dot - omega0**2*np.sin(theta)]
# 不同初始条件的相平面图
fig, ax = plt.subplots(figsize=(10, 8))
initial_conditions = np.linspace(-2*np.pi, 2*np.pi, 20)
colors = plt.cm.viridis(np.linspace(0, 1, len(initial_conditions)))
t = np.linspace(0, 20, 2000)
for theta0, color in zip(initial_conditions, colors):
y0 = [theta0, 0.0]
sol = self.solve(f, y0, t, 'rk4')
ax.plot(sol[:, 0], sol[:, 1], color=color, alpha=0.7, linewidth=0.8)
ax.set_xlabel('θ (Angle)')
ax.set_ylabel('θ\' (Angular Velocity)')
ax.set_title('Damped Pendulum Phase Portrait')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)
return fig
def lorenz_attractor(self, sigma=10, rho=28, beta=8/3):
"""
Lorenz吸引子（混沌系统示例）
dx/dt = σ(y - x)
dy/dt = x(ρ - z) - y
dz/dt = xy - βz
"""
def f(t, state):
x, y, z = state
return [
sigma * (y - x),
x * (rho - z) - y,
x * y - beta * z
]
# 微小差异的初始条件
t_span = (0, 50)
dt = 0.01
t = np.arange(t_span[0], t_span[1], dt)
fig = plt.figure(figsize=(14, 5))
# 3D轨迹
ax1 = fig.add_subplot(131, projection='3d')
colors = ['red', 'blue']
for i, delta in enumerate([0, 1e-5]):
y0 = [1.0 + delta, 1.0, 1.0]
sol = self.solve(f, y0, t, 'rk4')
ax1.plot(sol[:, 0], sol[:, 1], sol[:, 2],
color=colors[i], alpha=0.7, linewidth=0.5,
label=f'δ={delta}')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')
ax1.set_title('Lorenz Attractor (Chaos)')
ax1.legend()
# X随时间变化
ax2 = fig.add_subplot(132)
sol1 = self.solve(f, [1.0, 1.0, 1.0], t, 'rk4')
sol2 = self.solve(f, [1.00001, 1.0, 1.0], t, 'rk4')
ax2.plot(t, sol1[:, 0], 'r-', alpha=0.7, label='x(0)=1.0')
ax2.plot(t, sol2[:, 0], 'b-', alpha=0.7, label='x(0)=1.00001')
ax2.set_xlabel('t')
ax2.set_ylabel('x')
ax2.set_title('Sensitive Dependence on Initial Conditions')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 误差增长
ax3 = fig.add_subplot(133)
error = np.abs(sol1[:, 0] - sol2[:, 0])
ax3.semilogy(t, error, 'g-', linewidth=1)
ax3.set_xlabel('t')
ax3.set_ylabel('|x₁ - x₂|')
ax3.set_title('Error Growth (Exponential)')
ax3.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def compare_methods_accuracy(self):
"""
比较不同数值方法的精度
"""
# 测试方程：y' = -y, y(0) = 1
# 解析解：y(t) = e^(-t)
def f(t, y):
return -y
y0 = 1.0
t_end = 5
exact = np.exp(-t_end)
h_values = np.logspace(-3, -1, 20)
errors = {method: [] for method in ['euler', 'improved_euler', 'rk4']}
for h in h_values:
t = np.arange(0, t_end + h, h)
for method in errors.keys():
sol = self.solve(f, y0, t, method)
error = abs(sol[-1] - exact)
errors[method].append(error)
fig, ax = plt.subplots(figsize=(10, 6))
for method, errs in errors.items():
ax.loglog(h_values, errs, 'o-', label=method, markersize=4)
ax.loglog(h_values, h_values, 'k--', alpha=0.5, label='O(h)')
ax.loglog(h_values, h_values**2, 'k:', alpha=0.5, label='O(h²)')
ax.loglog(h_values, h_values**4, 'k-.', alpha=0.5, label='O(h⁴)')
ax.set_xlabel('Step size h')
ax.set_ylabel('Error at t=5')
ax.set_title('ODE Solver Accuracy Comparison')
ax.legend()
ax.grid(True, alpha=0.3)
return fig
# 使用示例
if __name__ == '__main__':
solver = ODESolver()
solver.harmonic_oscillator()
plt.show()
项目4：优化问题求解器
# project4_optimization.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, minimize_scalar
from mpl_toolkits.mplot3d import Axes3D
class OptimizationSolver:
"""
数值优化求解器
实现：
- 梯度下降法
- 牛顿法
- 约束优化
- 可视化收敛过程
"""
def gradient_descent(self, f, grad_f, x0, lr=0.1, tol=1e-6, max_iter=1000):
"""
梯度下降法
x_{k+1} = x_k - lr * ∇f(x_k)
最直观的优化方法，沿着最陡下降方向移动
"""
x = np.array(x0, dtype=float)
trajectory = [x.copy()]
for i in range(max_iter):
gradient = grad_f(x)
x_new = x - lr * gradient
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
trajectory.append(x.copy())
return x, np.array(trajectory)
def newton_method(self, f, grad_f, hess_f, x0, tol=1e-6, max_iter=100):
"""
牛顿法
x_{k+1} = x_k - H⁻¹(x_k) * ∇f(x_k)
利用二阶信息（Hessian），收敛更快但需要计算逆矩阵
"""
x = np.array(x0, dtype=float)
trajectory = [x.copy()]
for i in range(max_iter):
gradient = grad_f(x)
hessian = hess_f(x)
# 求解线性方程组 H * d = -∇f
direction = np.linalg.solve(hessian, -gradient)
x_new = x + direction
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
trajectory.append(x.copy())
return x, np.array(trajectory)
def newton_line_search(self, f, grad_f, hess_f, x0, tol=1e-6, max_iter=100):
"""
带线搜索的牛顿法
"""
x = np.array(x0, dtype=float)
trajectory = [x.copy()]
for i in range(max_iter):
gradient = grad_f(x)
hessian = hess_f(x)
direction = np.linalg.solve(hessian, -gradient)
# Armijo 线搜索
alpha = 1.0
c = 0.5
while f(x + alpha * direction) > f(x) + c * alpha * np.dot(gradient, direction):
alpha *= 0.5
x_new = x + alpha * direction
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
trajectory.append(x.copy())
return x, np.array(trajectory)
def visualize_2d_optimization(self, f, grad_f, x_range, y_range, method='gd'):
"""
可视化2D优化过程
"""
fig = plt.figure(figsize=(14, 6))
# 创建网格
x = np.linspace(x_range[0], x_range[1], 100)
y = np.linspace(y_range[0], y_range[1], 100)
X, Y = np.meshgrid(x, y)
Z = np.array([[f(np.array([xi, yi])) for xi in x] for yi in y])
# 3D表面图
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('f(x,y)')
ax1.set_title('Objective Function')
# 等高线图 + 优化轨迹
ax2 = fig.add_subplot(122)
contour = ax2.contour(X, Y, Z, levels=20, cmap='viridis')
ax2.clabel(contour, inline=True, fontsize=8)
# 运行优化
x0 = np.array([-1.5, 1.5])
if method == 'gd':
x_opt, traj = self.gradient_descent(f, grad_f, x0, lr=0.1)
title = 'Gradient Descent'
elif method == 'newton':
# 需要Hessian
def hess_f(x):
return np.array([
[2 + 4*x[0]**2 - 2*(1-x[0]**2), -2*x[0]*x[1]],
[-2*x[0]*x[1], 2]
])
x_opt, traj = self.newton_method(f, grad_f, hess_f, x0)
title = 'Newton Method'
# 绘制轨迹
ax2.plot(traj[:, 0], traj[:, 1], 'r.-', markersize=8, linewidth=1.5)
ax2.plot(traj[0, 0], traj[0, 1], 'go', markersize=12, label='Start')
ax2.plot(traj[-1, 0], traj[-1, 1], 'r*', markersize=15, label='End')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title(f'{title} Trajectory ({len(traj)} iterations)')
ax2.legend()
ax2.set_aspect('equal')
plt.tight_layout()
return fig
def rosenbrock_function(self):
"""
Rosenbrock函数（优化测试经典函数）
f(x,y) = (1-x)² + 100(y-x²)²
全局最小值在 (1, 1)，位于狭长山谷中
"""
def f(x):
return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2
def grad_f(x):
return np.array([
-2*(1-x[0]) - 400*x[0]*(x[1]-x[0]**2),
200*(x[1]-x[0]**2)
])
def hess_f(x):
return np.array([
[2 + 800*x[0]**2 - 400*(x[1]-x[0]**2), -400*x[0]],
[-400*x[0], 200]
])
return f, grad_f, hess_f
def constrained_optimization_demo(self):
"""
约束优化示例
最小化 f(x,y) = x² + y²
约束: x + y >= 1
"""
def objective(x):
return x[0]**2 + x[1]**2
# 不等式约束: g(x) >= 0
def constraint(x):
return x[0] + x[1] - 1
# 使用 scipy 求解
x0 = [0.5, 0.5]
cons = {'type': 'ineq', 'fun': constraint}
result = minimize(objective, x0, method='SLSQP', constraints=cons)
# 可视化
fig, ax = plt.subplots(figsize=(10, 8))
# 可行域
x = np.linspace(-1, 2, 100)
y = np.linspace(-1, 2, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + Y**2
# 等高线
contour = ax.contour(X, Y, Z, levels=15, cmap='viridis')
ax.clabel(contour, inline=True, fontsize=8)
# 约束边界
x_boundary = np.linspace(-1, 2, 100)
y_boundary = 1 - x_boundary
ax.plot(x_boundary, y_boundary, 'r-', linewidth=2, label='Constraint: x+y=1')
ax.fill_between(x_boundary, -1, y_boundary, where=(y_boundary > -1),
alpha=0.2, color='red', label='Infeasible region')
# 最优解
ax.plot(result.x[0], result.x[1], 'r*', markersize=20, label=f'Optimum {result.x.round(4)}')
ax.plot(x0[0], x0[1], 'go', markersize=12, label='Start')
ax.set_xlim([-1, 2])
ax.set_ylim([-1, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Constrained Optimization')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
return fig
def compare_convergence(self):
"""
比较梯度下降和牛顿法的收敛速度
"""
f, grad_f, hess_f = self.rosenbrock_function()
x0 = np.array([-1.2, 1.0])
# 梯度下降
_, traj_gd = self.gradient_descent(f, grad_f, x0, lr=0.001, max_iter=5000)
errors_gd = [f(x) for x in traj_gd]
# 牛顿法
_, traj_newton = self.newton_method(f, grad_f, hess_f, x0)
errors_newton = [f(x) for x in traj_newton]
fig, ax = plt.subplots(figsize=(10, 6))
ax.semilogy(errors_gd, 'b-', label=f'Gradient Descent ({len(traj_gd)} iters)')
ax.semilogy(errors_newton, 'r-', label=f'Newton Method ({len(traj_newton)} iters)')
ax.set_xlabel('Iteration')
ax.set_ylabel('Objective Value (log scale)')
ax.set_title('Convergence Speed Comparison on Rosenbrock Function')
ax.legend()
ax.grid(True, alpha=0.3)
return fig
# 使用示例
if __name__ == '__main__':
opt = OptimizationSolver()
# Rosenbrock函数优化
f, grad_f, hess_f = opt.rosenbrock_function()
opt.visualize_2d_optimization(f, grad_f, [-2, 2], [-1, 3], method='gd')
plt.show()
项目5：傅里叶分析工具
# project5_fourier_analysis.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq, fftshift
from scipy.signal import butter, filtfilt
class FourierAnalyzer:
"""
傅里叶分析工具箱
功能：
- DFT/FFT 实现
- 频谱分析
- 滤波器设计
- 信号重构
"""
def __init__(self, sample_rate=1000):
self.fs = sample_rate # 采样率 (Hz)
self.dt = 1 / sample_rate
def dft_slow(self, x):
"""
慢速 DFT（用于教学理解）
X[k] = Σ x[n] * e^(-j*2π*k*n/N)
复杂度 O(N²)，实际使用 FFT
"""
N = len(x)
X = np.zeros(N, dtype=complex)
for k in range(N):
for n in range(N):
X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
return X
def analyze_signal(self, signal, t=None):
"""
完整信号分析
Returns:
时域、频域、功率谱密度
"""
N = len(signal)
if t is None:
t = np.arange(N) * self.dt
# FFT
X = fft(signal)
# 频率轴
freqs = fftfreq(N, self.dt)
# 功率谱密度
psd = np.abs(X)**2 / (self.fs * N)
psd[1:-1] *= 2 # 双边谱转单边谱（除直流和奈奎斯特频率）
# 相位谱
phase = np.angle(X)
return {
'time': t,
'signal': signal,
'frequencies': freqs[:N//2],
'amplitude': np.abs(X)[:N//2] * 2 / N,
'power': psd[:N//2],
'phase': phase[:N//2]
}
def visualize_analysis(self, signal, t=None):
"""
可视化信号分析结果
"""
result = self.analyze_signal(signal, t)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 时域波形
ax1 = axes[0, 0]
ax1.plot(result['time'], result['signal'], 'b-', linewidth=1)
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Amplitude')
ax1.set_title('Time Domain Signal')
ax1.grid(True, alpha=0.3)
# 幅频谱
ax2 = axes[0, 1]
ax2.stem(result['frequencies'], result['amplitude'], basefmt=' ')
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Amplitude')
ax2.set_title('Amplitude Spectrum')
ax2.grid(True, alpha=0.3)
# 功率谱密度
ax3 = axes[1, 0]
ax3.semilogy(result['frequencies'], result['power'], 'g-', linewidth=1)
ax3.set_xlabel('Frequency (Hz)')
ax3.set_ylabel('Power/Frequency (dB/Hz)')
ax3.set_title('Power Spectral Density')
ax3.grid(True, alpha=0.3)
# 相位谱
ax4 = axes[1, 1]
ax4.plot(result['frequencies'], np.unwrap(result['phase']), 'r-', linewidth=1)
ax4.set_xlabel('Frequency (Hz)')
ax4.set_ylabel('Phase (rad)')
ax4.set_title('Phase Spectrum')
ax4.grid(True, alpha=0.3)
plt.tight_layout()
return fig, result
def filter_design(self, filter_type='lowpass', cutoff=100, order=4):
"""
Butterworth 滤波器设计
Args:
filter_type: 'lowpass', 'highpass', 'bandpass', 'bandstop'
cutoff: 截止频率 (Hz)，带通/带阻时为 [low, high]
order: 滤波器阶数
"""
nyquist = self.fs / 2
if isinstance(cutoff, (list, tuple)):
normalized = [c / nyquist for c in cutoff]
else:
normalized = cutoff / nyquist
b, a = butter(order, normalized, btype=filter_type)
return b, a
def apply_filter(self, signal, b, a):
"""
应用滤波器（零相位滤波，避免相位失真）
"""
return filtfilt(b, a, signal)
def denoising_demo(self):
"""
信号去噪演示
"""
# 生成含噪声信号
t = np.linspace(0, 1, self.fs, endpoint=False)
# 原始信号：多个正弦波
signal_clean = (np.sin(2 * np.pi * 50 * t) + # 50 Hz
0.5 * np.sin(2 * np.pi * 120 * t) + # 120 Hz
0.3 * np.sin(2 * np.pi * 200 * t)) # 200 Hz
# 添加噪声
noise = 0.5 * np.random.randn(len(t))
signal_noisy = signal_clean + noise
# 分析
fig, axes = plt.subplots(3, 2, figsize=(14, 12))
# 时域
ax1 = axes[0, 0]
ax1.plot(t[:200], signal_clean[:200], 'b-', label='Clean')
ax1.plot(t[:200], signal_noisy[:200], 'r-', alpha=0.5, label='Noisy')
ax1.legend()
ax1.set_title('Time Domain (First 200 samples)')
ax1.grid(True, alpha=0.3)
# 频域
freqs = fftfreq(len(t), self.dt)[:len(t)//2]
psd_clean = np.abs(fft(signal_clean))[:len(t)//2]**2
psd_noisy = np.abs(fft(signal_noisy))[:len(t)//2]**2
ax2 = axes[0, 1]
ax2.semilogy(freqs, psd_clean, 'b-', label='Clean')
ax2.semilogy(freqs, psd_noisy, 'r-', alpha=0.5, label='Noisy')
ax2.set_xlim([0, 300])
ax2.legend()
ax2.set_title('Power Spectrum')
ax2.grid(True, alpha=0.3)
# 低通滤波
b, a = self.filter_design('lowpass', cutoff=150, order=6)
signal_filtered = self.apply_filter(signal_noisy, b, a)
ax3 = axes[1, 0]
ax3.plot(t[:200], signal_clean[:200], 'b-', label='Clean', linewidth=2)
ax3.plot(t[:200], signal_filtered[:200], 'g-', label='Filtered', alpha=0.8)
ax3.legend()
ax3.set_title('Lowpass Filtered')
ax3.grid(True, alpha=0.3)
psd_filtered = np.abs(fft(signal_filtered))[:len(t)//2]**2
ax4 = axes[1, 1]
ax4.semilogy(freqs, psd_clean, 'b-', label='Clean')
ax4.semilogy(freqs, psd_filtered, 'g-', label='Filtered')
ax4.set_xlim([0, 300])
ax4.legend()
ax4.set_title('Filtered Spectrum')
ax4.grid(True, alpha=0.3)
# 带通滤波（只保留 100-150 Hz）
b_bp, a_bp = self.filter_design('bandpass', cutoff=[80, 160], order=6)
signal_bandpass = self.apply_filter(signal_noisy, b_bp, a_bp)
ax5 = axes[2, 0]
ax5.plot(t[:200], signal_bandpass[:200], 'purple', label='Bandpass')
ax5.legend()
ax5.set_title('Bandpass Filter (80-160 Hz)')
ax5.grid(True, alpha=0.3)
psd_bp = np.abs(fft(signal_bandpass))[:len(t)//2]**2
ax6 = axes[2, 1]
ax6.semilogy(freqs, psd_bp, 'purple', label='Bandpass')
ax6.set_xlim([0, 300])
ax6.legend()
ax6.set_title('Bandpass Spectrum')
ax6.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def reconstruction_demo(self):
"""
信号重构演示：从频域回到时域
"""
# 方波的傅里叶级数重构
t = np.linspace(-np.pi, np.pi, 1000)
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
axes = axes.flatten()
n_harmonics = [1, 3, 5, 10, 20, 50]
for idx, n in enumerate(n_harmonics):
ax = axes[idx]
# 方波傅里叶级数：Σ (4/π) * sin((2k+1)x)/(2k+1)
reconstruction = np.zeros_like(t)
for k in range(n):
harmonic = 2*k + 1
reconstruction += (4/np.pi) * np.sin(harmonic * t) / harmonic
ax.plot(t, np.sign(np.sin(t)), 'k--', linewidth=1, label='Square wave', alpha=0.5)
ax.plot(t, reconstruction, 'b-', linewidth=1.5, label=f'{n} harmonics')
ax.set_ylim([-1.5, 1.5])
ax.set_title(f'Fourier Series Reconstruction (n={n})')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
return fig
def spectrogram(self, signal, window_size=256, overlap=128):
"""
短时傅里叶变换（时频分析）
"""
hop = window_size - overlap
n_windows = (len(signal) - window_size) // hop + 1
freqs = fftfreq(window_size, self.dt)[:window_size//2]
times = np.arange(n_windows) * hop * self.dt
spectrogram = np.zeros((len(freqs), n_windows))
window = np.hanning(window_size)
for i in range(n_windows):
start = i * hop
segment = signal[start:start+window_size] * window
spectrum = np.abs(fft(segment))[:window_size//2]
spectrogram[:, i] = spectrum
return times, freqs, spectrogram
def visualize_spectrogram(self, signal, t=None):
"""
可视化时频谱
"""
times, freqs, spec = self.spectrogram(signal)
fig, ax = plt.subplots(figsize=(12, 6))
im = ax.pcolormesh(times, freqs, 20*np.log10(spec + 1e-10),
shading='gouraud', cmap='viridis')
ax.set_ylabel('Frequency (Hz)')
ax.set_xlabel('Time (s)')
ax.set_title('Spectrogram (STFT)')
ax.set_ylim([0, self.fs/2])
plt.colorbar(im, ax=ax, label='Magnitude (dB)')
return fig
# 使用示例
if __name__ == '__main__':
analyzer = FourierAnalyzer(sample_rate=1000)
# 生成测试信号
t = np.linspace(0, 1, 1000)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)
analyzer.visualize_analysis(signal, t)
plt.show()
代码规范与最佳实践
4.1 向量化编程
import numpy as np
# ❌ 慢：Python循环
def slow_sum(arr):
result = 0
for x in arr:
result += x
return result
# ✅ 快：NumPy向量化
def fast_sum(arr):
return np.sum(arr)
# ❌ 慢：逐元素操作
result = []
for i in range(len(x)):
result.append(np.sin(x[i]) * np.cos(x[i]))
# ✅ 快：向量化操作
result = np.sin(x) * np.cos(x)
# ✅ 广播机制
# (100, 1) + (1, 50) → (100, 50)
A = np.random.randn(100, 1)
B = np.random.randn(1, 50)
C = A + B # 自动广播
4.2 函数封装
class MathVisualizer:
"""
数学可视化工具类
提供常见数学概念的可视化方法
"""
def __init__(self, figsize=(10, 6), style='seaborn-v0_8'):
"""
初始化可视化器
Args:
figsize: 默认图像尺寸
style: Matplotlib样式
"""
self.figsize = figsize
plt.style.use(style)
def plot_function(self, f, x_range, **kwargs):
"""
绘制函数图像
Args:
f: 可调用函数
x_range: (xmin, xmax) 元组
**kwargs: 传递给plt.plot的参数
"""
x = np.linspace(x_range[0], x_range[1], 1000)
y = f(x)
fig, ax = plt.subplots(figsize=self.figsize)
ax.plot(x, y, **kwargs)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid(True, alpha=0.3)
return fig, ax
4.3 文档字符串规范
def numerical_derivative(f, x, h=1e-5, method='central'):
"""
计算函数的数值导数
使用有限差分方法近似计算函数在某点的导数。
中心差分法的误差为 O(h²)，优于前向差分的 O(h)。
Parameters
----------
f : callable
要求导的函数，接受标量输入返回标量输出
x : float
求导点
h : float, optional
差分步长，默认为 1e-5
过小会导致浮点误差，过大会增加截断误差
method : str, optional
差分方法：'forward', 'backward', 'central'
默认为 'central'（精度最高）
Returns
-------
float
导数的数值近似值
Raises
------
ValueError
当 method 不是支持的值时
Examples
--------
>>> f = lambda x: x**2
>>> numerical_derivative(f, 2.0)
4.000000000000001 # 接近理论值 4
>>> numerical_derivative(np.sin, 0.0)
1.0 # cos(0) = 1
References
----------
.. [1] 数值分析，李庆扬等，第5章
"""
if method == 'forward':
return (f(x + h) - f(x)) / h
elif method == 'backward':
return (f(x) - f(x - h)) / h
elif method == 'central':
return (f(x + h) - f(x - h)) / (2 * h)
else:
raise ValueError(f"Unknown method: {method}")
4.4 单元测试
import unittest
import numpy as np
from numpy.testing import assert_almost_equal, assert_allclose
class TestNumericalMethods(unittest.TestCase):
"""
数值方法单元测试
"""
def setUp(self):
"""测试前的准备工作"""
self.tolerance = 1e-6
def test_numerical_derivative_polynomial(self):
"""测试多项式函数的数值微分"""
f = lambda x: 3*x**3 + 2*x**2 + x + 1
df_exact = lambda x: 9*x**2 + 4*x + 1
x_test = 2.0
from numerical_methods import numerical_derivative
df_numerical = numerical_derivative(f, x_test)
assert_almost_equal(df_numerical, df_exact(x_test), decimal=5)
def test_numerical_derivative_trig(self):
"""测试三角函数的数值微分"""
f = np.sin
df_exact = np.cos
x_test = np.pi / 4
from numerical_methods import numerical_derivative
df_numerical = numerical_derivative(f, x_test)
assert_almost_equal(df_numerical, df_exact(x_test), decimal=8)
def test_numerical_integration(self):
"""测试数值积分"""
from numerical_methods import NumericalIntegration
integrator = NumericalIntegration()
# ∫[0,1] x² dx = 1/3
f = lambda x: x**2
result = integrator.simpson(f, 0, 1, n=100)
assert_almost_equal(result, 1/3, decimal=8)
if __name__ == '__main__':
unittest.main()
附录：快速参考表
常用NumPy操作
| 操作 | 代码 |
|---|---|
| 创建等差数列 | np.linspace(a, b, n) |
| 创建等比数列 | np.logspace(a, b, n) |
| 矩阵乘法 | A @ B 或 np.dot(A, B) |
| 元素乘法 | A * B |
| 矩阵求逆 | np.linalg.inv(A) |
| 特征值 | np.linalg.eigvals(A) |
| 条件筛选 | A[A > 0] |
| 广播 | A[:, None] + B[None, :] |
常用SciPy模块
from scipy import integrate # 积分
from scipy import optimize # 优化
from scipy import linalg # 线性代数
from scipy import interpolate # 插值
from scipy import fft # 傅里叶变换
from scipy import signal # 信号处理
from scipy import stats # 统计
常用SymPy操作
import sympy as sp
x, y, z = sp.symbols('x y z')
f = sp.sin(x)**2 + sp.cos(x)**2 # 符号表达式
sp.simplify(f) # 化简 → 1
sp.diff(f, x) # 求导
sp.integrate(f, x) # 积分
sp.series(f, x, 0, 5) # 泰勒展开
sp.solve(x**2 - 1, x) # 解方程 → [-1, 1]
sp.lambdify(x, f, 'numpy') # 转为Python函数
总结
核心原则
- 代码验证数学：用数值实验验证理论结果
- 可视化优先：一张图胜过千言万语
- 模块化设计：封装可复用的数学工具
- 测试驱动：确保数值方法的正确性
学习路径
- 第1周：熟悉Python科学计算环境
- 第2-3周：实现基本数值方法（微分、积分）
- 第4-5周：完成实战项目
- 第6周：优化代码，添加文档和测试
扩展阅读
- 《Numerical Recipes》(William H. Press)
- 《Python科学计算》(张若愚)
- 《Scipy Lecture Notes》(在线免费)
- MIT 18.06 (线性代数)
- MIT 18.085 (计算科学与工程)
"用代码验证数学，让抽象变具体。这是第一性原理的力量。" —— 费曼精神