SciPy

目录

  1. 简介
  2. 安装与配置
  3. 核心模块概览
  4. scipy.optimize - 优化与拟合
  5. scipy.integrate - 积分计算
  6. scipy.interpolate - 插值
  7. scipy.linalg - 线性代数
  8. scipy.stats - 统计
  9. scipy.signal - 信号处理
  10. scipy.sparse - 稀疏矩阵
  11. 综合实例
  12. 总结与资源

简介

SciPy（Scientific Python）是一个开源的Python库，用于科学和技术计算。它构建在NumPy之上，提供了大量用于数学、科学和工程计算的算法和工具。

SciPy 与 NumPy 的关系

  • NumPy: 提供基础的多维数组对象和基本操作
  • SciPy: 在NumPy基础上提供更高级的数学算法和工具

主要特点

  • 丰富的科学计算功能
  • 与NumPy无缝集成
  • 高效的C/Fortran底层实现
  • 广泛的文档和社区支持

安装与配置

安装方法

# 使用 pip 安装
pip install scipy

# 使用 conda 安装（推荐）
conda install scipy

# 安装完整科学计算环境
pip install numpy scipy matplotlib pandas

验证安装

import scipy
print(scipy.__version__)

基本导入方式

# 导入整个库
import scipy

# 导入特定模块
from scipy import optimize, integrate, stats

# 常用组合导入
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

核心模块概览

模块功能描述
scipy.optimize优化和根查找算法
scipy.integrate积分和微分方程求解
scipy.interpolate插值和平滑样条
scipy.linalg线性代数运算
scipy.stats统计分布和检验
scipy.signal信号处理工具
scipy.sparse稀疏矩阵处理
scipy.fft快速傅里叶变换
scipy.special特殊数学函数

scipy.optimize - 优化与拟合

1. 函数最小化

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

# 定义目标函数：Rosenbrock函数（测试优化算法的经典函数）
def rosenbrock(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

# 初始猜测点
x0 = np.array([-1, 2])

# 使用 BFGS 算法进行无约束最小化
result = optimize.minimize(rosenbrock, x0, method='BFGS')

print(f"最小值点: {result.x}")
print(f"最小值: {result.fun}")
print(f"迭代次数: {result.nit}")
print(f"是否成功: {result.success}")

2. 曲线拟合

from scipy.optimize import curve_fit

# 定义拟合函数模型
def func(x, a, b, c):
    """指数衰减函数模型"""
    return a * np.exp(-b * x) + c

# 生成模拟数据
x_data = np.linspace(0, 4, 50)
y_true = func(x_data, 2.5, 1.3, 0.5)
noise = np.random.normal(0, 0.1, len(x_data))
y_data = y_true + noise

# 进行曲线拟合
popt, pcov = curve_fit(func, x_data, y_data)

print(f"拟合参数: a={popt[0]:.3f}, b={popt[1]:.3f}, c={popt[2]:.3f}")

# 计算参数的标准差
perr = np.sqrt(np.diag(pcov))
print(f"参数标准差: {perr}")

3. 根查找

from scipy.optimize import fsolve, brentq

# 定义方程
def equation(x):
    return x**3 - 2*x**2 - 5*x + 6

# 使用 fsolve 求解（需要初始猜测）
root1 = fsolve(equation, x0=0)
print(f"根1: {root1[0]:.6f}")

# 使用 brentq 在区间内查找（更稳健）
root2 = brentq(equation, -2, 0)
root3 = brentq(equation, 0, 3)
print(f"根2: {root2:.6f}")
print(f"根3: {root3:.6f}")

4. 约束优化

# 定义目标函数
def objective(x):
    return x[0]**2 + x[1]**2

# 定义约束条件
def constraint1(x):
    return x[0] + x[1] - 1  # 等式约束: x[0] + x[1] = 1

# 约束定义
constraints = {'type': 'eq', 'fun': constraint1}
bounds = [(0, None), (0, None)]  # x[0] >= 0, x[1] >= 0

# 求解
result = optimize.minimize(
    objective, 
    x0=[0.5, 0.5], 
    method='SLSQP',
    bounds=bounds,
    constraints=constraints
)

print(f"约束优化结果: x = {result.x}")
print(f"目标函数值: {result.fun}")

scipy.integrate - 积分计算

1. 数值积分

from scipy import integrate
import numpy as np

# 定义被积函数
def integrand(x):
    return np.exp(-x**2)

# 一重积分
result, error = integrate.quad(integrand, 0, np.inf)
print(f"∫[0,∞] e^(-x²) dx = {result:.10f}")
print(f"误差估计: {error:.2e}")
print(f"理论值 √π/2 = {np.sqrt(np.pi)/2:.10f}")

# 二重积分
def integrand2(y, x):
    return x * y

result2, error2 = integrate.dblquad(integrand2, 0, 1, 0, 1)
print(f"二重积分结果: {result2}")

2. 常微分方程求解

from scipy.integrate import odeint, solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# 定义ODE: dy/dt = -2y
def model(y, t):
    return -2 * y

# 初始条件和时间点
y0 = 1
t = np.linspace(0, 5, 100)

# 使用 odeint 求解
y = odeint(model, y0, t)

# 使用 solve_ivp（推荐新方法）
sol = solve_ivp(lambda t, y: -2*y, [0, 5], [y0], t_eval=t)

# 绘制结果
plt.plot(t, y, 'b-', label='odeint')
plt.plot(sol.t, sol.y[0], 'r--', label='solve_ivp')
plt.plot(t, np.exp(-2*t), 'g:', label='解析解')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.title('ODE 求解比较')

3. 洛伦兹吸引子（复杂ODE示例）

def lorenz(t, state, sigma=10, rho=28, beta=8/3):
    """洛伦兹系统"""
    x, y, z = state
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

# 初始条件
initial_state = [1.0, 1.0, 1.0]
t_span = [0, 50]
t_eval = np.linspace(0, 50, 10000)

# 求解
sol = solve_ivp(lorenz, t_span, initial_state, t_eval=t_eval, method='RK45')

# 3D 可视化
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(sol.y[0], sol.y[1], sol.y[2], linewidth=0.5)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Lorenz Attractor')

scipy.interpolate - 插值

1. 一维插值

from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

# 原始数据点
x = np.linspace(0, 10, 10)
y = np.sin(x)

# 更精细的插值点
x_new = np.linspace(0, 10, 100)

# 线性插值
f_linear = interpolate.interp1d(x, y, kind='linear')
y_linear = f_linear(x_new)

# 三次样条插值
cs = interpolate.CubicSpline(x, y)
y_cubic = cs(x_new)

# 绘制比较
plt.scatter(x, y, color='red', s=50, label='原始数据')
plt.plot(x_new, np.sin(x_new), 'k-', label='真实函数')
plt.plot(x_new, y_linear, 'b--', label='线性插值')
plt.plot(x_new, y_cubic, 'g:', label='三次样条')
plt.legend()
plt.title('插值方法比较')

2. 二维插值

from scipy.interpolate import griddata

# 生成随机散点数据
np.random.seed(42)
points = np.random.rand(100, 2)
values = np.sin(points[:, 0] * 10) * np.cos(points[:, 1] * 10)

# 创建规则网格
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:100j]

# 不同插值方法
grid_z_nearest = griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z_linear = griddata(points, values, (grid_x, grid_y), method='linear')
grid_z_cubic = griddata(points, values, (grid_x, grid_y), method='cubic')

# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
methods = ['Nearest', 'Linear', 'Cubic']
grids = [grid_z_nearest, grid_z_linear, grid_z_cubic]

for ax, method, grid in zip(axes, methods, grids):
    ax.imshow(grid.T, extent=(0, 1, 0, 1), origin='lower', cmap='viridis')
    ax.scatter(points[:, 0], points[:, 1], c='red', s=10)
    ax.set_title(f'{method} Interpolation')

scipy.linalg - 线性代数

1. 基本矩阵运算

from scipy import linalg
import numpy as np

# 创建矩阵
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 10]])

# 行列式
det = linalg.det(A)
print(f"行列式: {det:.4f}")

# 逆矩阵
A_inv = linalg.inv(A)
print(f"逆矩阵:\n{A_inv}")

# 验证
print(f"A @ A_inv =\n{np.round(A @ A_inv, 10)}")

# 矩阵范数
print(f"Frobenius范数: {linalg.norm(A)}")
print(f"2-范数: {linalg.norm(A, 2)}")

2. 特征值分解

# 特征值和特征向量
eigenvalues, eigenvectors = linalg.eig(A)

print("特征值:")
print(eigenvalues)
print("\n特征向量:")
print(eigenvectors)

# 验证: A @ v = λ * v
for i in range(len(eigenvalues)):
    lhs = A @ eigenvectors[:, i]
    rhs = eigenvalues[i] * eigenvectors[:, i]
    print(f"\n特征向量 {i+1} 验证: {np.allclose(lhs, rhs)}")

3. 线性方程组求解

# 求解 Ax = b
b = np.array([1, 2, 3])

# 方法1: 使用 solve
x = linalg.solve(A, b)
print(f"解 x: {x}")
print(f"验证 Ax = b: {np.allclose(A @ x, b)}")

# 方法2: LU分解
lu, piv = linalg.lu_factor(A)
x_lu = linalg.lu_solve((lu, piv), b)
print(f"LU分解解: {x_lu}")

# 方法3: 对于大规模稀疏系统，使用稀疏矩阵求解器
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve

A_sparse = csr_matrix(A)
x_sparse = spsolve(A_sparse, b)
print(f"稀疏求解器解: {x_sparse}")

4. 奇异值分解 (SVD)

# SVD分解
U, s, Vh = linalg.svd(A)

print(f"U 矩阵:\n{U}")
print(f"奇异值: {s}")
print(f"Vh 矩阵:\n{Vh}")

# 重构原矩阵
S = np.diag(s)
A_reconstructed = U @ S @ Vh
print(f"\n重构矩阵:\n{A_reconstructed}")
print(f"重构误差: {np.max(np.abs(A - A_reconstructed)):.2e}")

scipy.stats - 统计

1. 概率分布

from scipy import stats
import numpy as np
import matplotlib.pyplot as plt

# 正态分布
mu, sigma = 0, 1
normal_dist = stats.norm(mu, sigma)

# 概率密度函数 (PDF)
x = np.linspace(-5, 5, 100)
pdf = normal_dist.pdf(x)

# 累积分布函数 (CDF)
cdf = normal_dist.cdf(x)

# 随机采样
samples = normal_dist.rvs(size=1000)

# 绘制
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(x, pdf, 'b-', label='PDF')
axes[0].hist(samples, bins=30, density=True, alpha=0.7, label='Samples')
axes[0].legend()
axes[0].set_title('Normal Distribution PDF')

axes[1].plot(x, cdf, 'r-', label='CDF')
axes[1].legend()
axes[1].set_title('Normal Distribution CDF')

2. 统计检验

# 生成两组数据
np.random.seed(42)
group1 = np.random.normal(100, 10, 100)  # 均值100，标准差10
group2 = np.random.normal(105, 10, 100)  # 均值105，标准差10

# t检验
t_stat, p_value = stats.ttest_ind(group1, group2)
print(f"t统计量: {t_stat:.4f}")
print(f"p值: {p_value:.4f}")
print(f"显著性水平 α=0.05: {'拒绝原假设' if p_value < 0.05 else '不拒绝原假设'}")

# 单样本t检验
t_stat2, p_value2 = stats.ttest_1samp(group1, 100)
print(f"\n单样本t检验:")
print(f"t统计量: {t_stat2:.4f}, p值: {p_value2:.4f}")

# 卡方检验
observed = np.array([10, 20, 30])
expected = np.array([15, 15, 30])
chi2, p_chi2 = stats.chisquare(observed, expected)
print(f"\n卡方检验: χ²={chi2:.4f}, p={p_chi2:.4f}")

3. 描述性统计

# 生成数据
data = np.random.normal(50, 15, 1000)

# 描述性统计
print(f"均值: {np.mean(data):.2f}")
print(f"中位数: {np.median(data):.2f}")
print(f"标准差: {np.std(data):.2f}")
print(f"方差: {np.var(data):.2f}")
print(f"偏度: {stats.skew(data):.4f}")
print(f"峰度: {stats.kurtosis(data):.4f}")

# 百分位数
percentiles = [25, 50, 75]
print(f"\n百分位数 {percentiles}: {np.percentile(data, percentiles)}")

# 综合描述统计
desc = stats.describe(data)
print(f"\n综合统计:\n{desc}")

scipy.signal - 信号处理

1. 滤波器设计

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# 生成含噪声的信号
t = np.linspace(0, 1, 1000, endpoint=False)
clean_signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 50 * t)
noise = np.random.normal(0, 0.3, len(t))
noisy_signal = clean_signal + noise

# 设计低通滤波器
sos = signal.butter(10, 20, 'low', fs=1000, output='sos')
filtered_signal = signal.sosfilt(sos, noisy_signal)

# 绘制
plt.figure(figsize=(12, 6))
plt.plot(t, noisy_signal, 'b-', alpha=0.5, label='含噪声信号')
plt.plot(t, filtered_signal, 'r-', label='滤波后信号')
plt.plot(t, clean_signal, 'g--', label='原始信号')
plt.xlabel('时间 [s]')
plt.ylabel('幅值')
plt.legend()
plt.title('低通滤波效果')

2. 快速傅里叶变换 (FFT)

from scipy.fft import fft, fftfreq

# 采样参数
fs = 1000  # 采样频率
T = 1/fs   # 采样周期
t = np.linspace(0, 1, fs, endpoint=False)

# 复合信号
signal_data = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

# FFT
yf = fft(signal_data)
xf = fftfreq(fs, T)[:fs//2]

# 绘制频谱
plt.figure(figsize=(12, 4))
plt.plot(xf, 2.0/fs * np.abs(yf[0:fs//2]))
plt.xlabel('频率 [Hz]')
plt.ylabel('幅值')
plt.title('频谱分析')
plt.grid()

3. 卷积和相关

# 信号
x = np.array([1, 2, 3, 4, 5])
h = np.array([1, 0, -1])

# 卷积
conv_result = signal.convolve(x, h, mode='full')
print(f"卷积结果: {conv_result}")

# 互相关
corr_result = signal.correlate(x, h, mode='full')
print(f"互相关结果: {corr_result}")

# 自相关
auto_corr = signal.correlate(x, x, mode='full')
print(f"自相关结果: {auto_corr}")

scipy.sparse - 稀疏矩阵

1. 稀疏矩阵创建

from scipy.sparse import csr_matrix, csc_matrix, coo_matrix
import numpy as np

# 密集矩阵
dense_matrix = np.array([
    [1, 0, 0, 2],
    [0, 3, 0, 0],
    [0, 0, 4, 0],
    [5, 0, 0, 6]
])

# 转换为稀疏矩阵格式
sparse_csr = csr_matrix(dense_matrix)
sparse_csc = csc_matrix(dense_matrix)
sparse_coo = coo_matrix(dense_matrix)

print(f"原始矩阵大小: {dense_matrix.nbytes} bytes")
print(f"CSR格式数据: {sparse_csr.data}")
print(f"CSR格式索引: {sparse_csr.indices}")
print(f"CSR格式指针: {sparse_csr.indptr}")

# 从COO格式创建（适合构建）
row = np.array([0, 1, 2, 3, 3])
col = np.array([0, 1, 2, 0, 3])
data = np.array([1, 3, 4, 5, 6])
coo = coo_matrix((data, (row, col)), shape=(4, 4))
print(f"\nCOO矩阵:\n{coo.toarray()}")

2. 稀疏矩阵运算

from scipy.sparse.linalg import spsolve, eigsh

# 创建大型稀疏矩阵
n = 1000
diagonals = [np.ones(n), -2*np.ones(n), np.ones(n)]
A_sparse = sparse.diags(diagonals, [-1, 0, 1], format='csr')

# 稀疏矩阵向量乘法
x = np.random.rand(n)
b = A_sparse @ x

# 求解稀疏线性系统
x_solved = spsolve(A_sparse, b)
print(f"求解误差: {np.max(np.abs(x - x_solved)):.2e}")

# 计算部分特征值（对于大型矩阵更高效）
eigenvalues, eigenvectors = eigsh(A_sparse, k=5, which='SM')
print(f"\n最小5个特征值: {eigenvalues}")

综合实例

实例1：数据拟合与预测

import numpy as np
from scipy.optimize import curve_fit
from scipy import stats
import matplotlib.pyplot as plt

# 模拟实验数据（带噪声）
np.random.seed(42)
x_data = np.linspace(0, 10, 50)
true_params = [2.5, -1.0, 0.5]  # 真实参数

# 多项式模型 + 噪声
def polynomial(x, a, b, c):
    return a * x**2 + b * x + c

y_true = polynomial(x_data, *true_params)
y_data = y_true + np.random.normal(0, 2, len(x_data))

# 拟合
popt, pcov = curve_fit(polynomial, x_data, y_data)

# 计算R²
y_pred = polynomial(x_data, *popt)
ss_res = np.sum((y_data - y_pred)**2)
ss_tot = np.sum((y_data - np.mean(y_data))**2)
r_squared = 1 - (ss_res / ss_tot)

# 置信区间
perr = np.sqrt(np.diag(pcov))

print("拟合结果:")
print(f"参数: a={popt[0]:.3f}±{perr[0]:.3f}, b={popt[1]:.3f}±{perr[1]:.3f}, c={popt[2]:.3f}±{perr[2]:.3f}")
print(f"R² = {r_squared:.4f}")

# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='实验数据', alpha=0.6)
plt.plot(x_data, y_true, 'g--', label='真实曲线')
plt.plot(x_data, y_pred, 'r-', label='拟合曲线')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('多项式拟合')

实例2：物理系统仿真

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# 单摆系统（非线性）
def pendulum(t, y, g=9.81, L=1.0):
    """
    y[0] = 角度 θ
    y[1] = 角速度 ω
    """
    theta, omega = y
    dtheta_dt = omega
    domega_dt = -(g/L) * np.sin(theta)
    return [dtheta_dt, domega_dt]

# 初始条件：大角度（非线性效应明显）
y0 = [np.pi/2, 0]  # 90度，静止释放
t_span = [0, 10]
t_eval = np.linspace(0, 10, 1000)

# 求解
sol = solve_ivp(pendulum, t_span, y0, t_eval=t_eval, method='RK45')

# 可视化
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# 角度随时间变化
axes[0].plot(sol.t, sol.y[0])
axes[0].set_xlabel('时间 [s]')
axes[0].set_ylabel('角度 [rad]')
axes[0].set_title('单摆运动')
axes[0].grid()

# 相空间图（角度 vs 角速度）
axes[1].plot(sol.y[0], sol.y[1])
axes[1].set_xlabel('角度 [rad]')
axes[1].set_ylabel('角速度 [rad/s]')
axes[1].set_title('相空间图')
axes[1].grid()

plt.tight_layout()

实例3：图像处理（结合scipy.ndimage）

from scipy import ndimage
import numpy as np
import matplotlib.pyplot as plt

# 创建测试图像
np.random.seed(42)
image = np.zeros((100, 100))
image[30:70, 30:70] = 1  # 正方形
image += np.random.normal(0, 0.1, image.shape)

# 高斯滤波
filtered = ndimage.gaussian_filter(image, sigma=2)

# 边缘检测
sobel_x = ndimage.sobel(image, axis=0)
sobel_y = ndimage.sobel(image, axis=1)
sobel_magnitude = np.sqrt(sobel_x**2 + sobel_y**2)

# 形态学操作
from scipy.ndimage import binary_erosion, binary_dilation
binary = image > 0.5
eroded = binary_erosion(binary)
dilated = binary_dilation(binary)

# 可视化
fig, axes = plt.subplots(2, 3, figsize=(12, 8))
axes[0, 0].imshow(image, cmap='gray')
axes[0, 0].set_title('原始图像')
axes[0, 1].imshow(filtered, cmap='gray')
axes[0, 1].set_title('高斯滤波')
axes[0, 2].imshow(sobel_magnitude, cmap='gray')
axes[0, 2].set_title('Sobel边缘检测')
axes[1, 0].imshow(binary, cmap='gray')
axes[1, 0].set_title('二值化')
axes[1, 1].imshow(eroded, cmap='gray')
axes[1, 1].set_title('腐蚀')
axes[1, 2].imshow(dilated, cmap='gray')
axes[1, 2].set_title('膨胀')

for ax in axes.flat:
    ax.axis('off')

总结与资源

SciPy 核心要点

  1. 优化 (optimize): 函数最小化、曲线拟合、根查找
  2. 积分 (integrate): 数值积分、ODE求解
  3. 插值 (interpolate): 一维/二维插值、样条
  4. 线性代数 (linalg): 矩阵运算、分解、求解
  5. 统计 (stats): 分布、检验、描述统计
  6. 信号处理 (signal): 滤波、FFT、卷积
  7. 稀疏矩阵 (sparse): 高效存储和运算

学习资源

  • 官方文档: https://docs.scipy.org/doc/scipy/
  • SciPy Lecture Notes: https://scipy-lectures.org/
  • GitHub: https://github.com/scipy/scipy
  • Stack Overflow: 搜索 “scipy” 标签

最佳实践

  1. 优先使用 NumPy: 基础数组操作使用 NumPy
  2. 选择合适的算法: 根据问题特点选择优化方法
  3. 注意数值稳定性: 处理病态问题时谨慎
  4. 利用稀疏性: 大规模稀疏问题使用稀疏矩阵
  5. 验证结果: 重要计算进行结果验证

与其他库的配合

库用途
NumPy基础数组操作
Matplotlib数据可视化
Pandas数据处理
Scikit-learn机器学习
SymPy符号计算

本教程基于 SciPy 1.10+ 版本编写，建议保持库版本更新以获得最佳体验。

# scipy
Logo
AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

cover

开源鸿蒙跨平台直播｜KMP/CMP 鸿蒙适配，从共享代码到真机运行

avatar AtomGit开源社区
cover

中国电信 Xing4.0-29B-A4B 上线 AtomGit AI

avatar AtomGit开源社区
cover

奖金加码直达70万！中国电信星辰杯高校 AI 算子开发挑战赛火热进行中

avatar AtomGit开源社区