NumPy 与 SciPy 科学计算

SciPy 是建立在 NumPy 基础上的科学计算库，提供了许多高级算法和函数，涵盖优化、积分、插值、信号处理、线性代数、统计等模块。NumPy 提供高效的多维数组，而 SciPy 则利用这些数组实现科学计算算法。两者结合，构成了 Python 科学计算的核心生态。

1. SciPy 概述

SciPy 包含多个子模块，每个模块专注于特定领域。常用的模块包括：

scipy.integrate：数值积分，求解常微分方程
scipy.optimize：优化算法，曲线拟合，求根
scipy.interpolate：插值（一维和多维）
scipy.linalg：线性代数（扩展了 NumPy 的 linalg）
scipy.stats：统计分布和检验
scipy.signal：信号处理（滤波、卷积等）
scipy.sparse：稀疏矩阵处理

所有 SciPy 函数都接受 NumPy 数组作为输入，并返回 NumPy 数组或标量。

2. 数值积分（scipy.integrate）

使用 quad 计算定积分，使用 odeint 或 solve_ivp 求解常微分方程。

2.1 定积分

import numpy as np
from scipy.integrate import quad

# 定义被积函数
def f(x):
    return x ** 2

# 计算 ∫₀¹ x² dx
result, error = quad(f, 0, 1)
print(f"积分结果: {result}, 估计误差: {error}")

输出：

积分结果: 0.33333333333333337, 估计误差: 3.700743415417188e-15

2.2 常微分方程

from scipy.integrate import odeint

# 定义微分方程 dy/dt = -2y
def model(y, t):
    dydt = -2 * y
    return dydt

# 时间点
t = np.linspace(0, 5, 100)
# 初始条件
y0 = 1
# 求解
y = odeint(model, y0, t)

# 绘制结果（可选）
import matplotlib.pyplot as plt
plt.plot(t, y)
plt.xlabel('t')
plt.ylabel('y')
plt.title('dy/dt = -2y')
plt.show()

3. 优化（scipy.optimize）

包含最小化、曲线拟合、求根等功能。

3.1 曲线拟合

from scipy.optimize import curve_fit

# 生成带噪声的数据
xdata = np.linspace(0, 10, 50)
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

np.random.seed(42)
ydata = func(xdata, 2.5, 1.3, 0.5) + 0.2 * np.random.normal(size=len(xdata))

# 拟合
popt, pcov = curve_fit(func, xdata, ydata, p0=[1, 1, 1])
print("拟合参数:", popt)   # 应该接近 [2.5, 1.3, 0.5]

3.2 最小化

from scipy.optimize import minimize

# 定义 Rosenbrock 函数
def rosen(x):
    return sum(100.0 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)

# 初始猜测
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead')
print(res.x)

4. 线性代数（scipy.linalg）

scipy.linalg 包含所有 numpy.linalg 的函数，并增加了更多高级功能（如 LU 分解、Schur 分解等）。通常 SciPy 的线性代数函数更快、更稳定。

from scipy import linalg

A = np.array([[1, 2], [3, 4]])
# 矩阵逆
A_inv = linalg.inv(A)
print("逆矩阵:\n", A_inv)

# 行列式
det = linalg.det(A)
print("行列式:", det)

# LU 分解
P, L, U = linalg.lu(A)
print("L:\n", L)
print("U:\n", U)

5. 插值（scipy.interpolate）

提供多种插值方法，如一维插值 interp1d、样条插值 splrep/splev、网格插值 griddata 等。

from scipy.interpolate import interp1d

x = np.linspace(0, 10, 10)
y = np.sin(x)

# 创建线性插值函数
f_linear = interp1d(x, y, kind='linear')
# 创建三次样条插值
f_cubic = interp1d(x, y, kind='cubic')

# 在新点上评估
x_new = np.linspace(0, 10, 50)
y_linear = f_linear(x_new)
y_cubic = f_cubic(x_new)

6. 统计（scipy.stats）

scipy.stats 提供了大量概率分布、统计检验和描述性统计函数。

from scipy import stats

# 生成正态分布随机数
rvs = stats.norm.rvs(loc=0, scale=1, size=1000, random_state=42)

# 计算描述性统计
mean = np.mean(rvs)
std = np.std(rvs)
print(f"样本均值: {mean:.3f}, 样本标准差: {std:.3f}")

# 进行 t 检验（检验均值是否为 0）
t_stat, p_value = stats.ttest_1samp(rvs, 0)
print(f"t 统计量: {t_stat:.3f}, p 值: {p_value:.3f}")

# 拟合数据到正态分布
params = stats.norm.fit(rvs)
print("拟合参数 (loc, scale):", params)

7. 信号处理（scipy.signal）

包括滤波、卷积、频谱分析等。

from scipy import signal

# 生成一个带噪声的信号
t = np.linspace(0, 1, 500)
clean = np.sin(2 * np.pi * 5 * t)
noise = np.random.normal(0, 0.5, len(t))
noisy = clean + noise

# 设计一个低通滤波器
sos = signal.butter(10, 10, 'lowpass', fs=500, output='sos')
filtered = signal.sosfilt(sos, noisy)

# 绘制对比
import matplotlib.pyplot as plt
plt.plot(t, noisy, label='带噪声')
plt.plot(t, filtered, label='滤波后')
plt.legend()
plt.show()

8. 稀疏矩阵（scipy.sparse）

当矩阵中零元素占绝大多数时，使用稀疏矩阵可以节省内存和计算时间。SciPy 支持多种稀疏格式（CSR、CSC、COO 等），并能与 NumPy 数组相互转换。

from scipy.sparse import csr_matrix

# 创建密集矩阵
dense = np.array([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
# 转换为 CSR 稀疏矩阵
sparse = csr_matrix(dense)
print("稀疏矩阵:\n", sparse)

# 转换回密集
dense_again = sparse.toarray()
print("恢复为密集:\n", dense_again)

提示： 在使用 SciPy 的函数时，大多数情况下直接传入 NumPy 数组即可。SciPy 与 NumPy 配合紧密，返回的结果也通常是 NumPy 数组，便于后续处理。

总结

SciPy 扩展了 NumPy 的功能，提供了大量科学计算所需的算法。通过将 NumPy 数组作为数据容器，SciPy 使得复杂的数值计算变得简单而高效。掌握这些常用子模块，可以让你在处理科学和工程问题时游刃有余。

上一章: 与 Matplotlib 可视化下一章: 并行计算基础

NumPy教程