Python案例如何用Scipy做积分计算

wen python案例 2026-07-16 2

本文目录导读：

基础安装与导入
一维定积分 - quad
广义积分
多重积分
实际应用案例
数值积分方法对比
带奇异点的积分处理
常用技巧

我来给你介绍几个使用SciPy进行积分计算的实用案例。

基础安装与导入

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

一维定积分 - `quad`

案例1：计算简单函数积分

# 计算 ∫(0 to 2) x² dx
def f(x):
    return x**2
result, error = integrate.quad(f, 0, 2)
print(f"∫(0 to 2) x² dx = {result:.6f}")
print(f"估计误差: {error:.2e}")
# 理论值: 2.666667

案例2：三角函数积分

# 计算 ∫(0 to π) sin(x) dx
result, error = integrate.quad(np.sin, 0, np.pi)
print(f"∫(0 to π) sin(x) dx = {result:.6f}")
# 理论值: 2.0

案例3：含参数函数积分

# 计算 ∫(0 to 1) 1/(1+x^a) dx
def f(x, a, b):
    return 1 / (1 + x**a) * b
result, error = integrate.quad(f, 0, 1, args=(2, 3))
print(f"∫(0 to 1) 3/(1+x²) dx = {result:.6f}")

广义积分

案例4：无穷积分

# 计算 ∫(0 to ∞) exp(-x²) dx
result, error = integrate.quad(lambda x: np.exp(-x**2), 0, np.inf)
print(f"∫(0 to ∞) exp(-x²) dx = {result:.6f}")
# 理论值: √π/2 ≈ 0.886227

案例5：奇点积分

# 计算 ∫(0 to 1) 1/√x dx (在x=0有奇点)
result, error = integrate.quad(lambda x: 1/np.sqrt(x), 0, 1)
print(f"∫(0 to 1) 1/√x dx = {result:.6f}")
# 理论值: 2.0

多重积分

案例6：二重积分

# 计算 ∫∫(x² + y²) dxdy, 区域: x∈[0,1], y∈[0,1]
def f2(x, y):
    return x**2 + y**2
# 方法1：使用nquad
result, error = integrate.nquad(f2, [[0, 1], [0, 1]])
print(f"二重积分结果: {result:.6f}")
# 方法2：使用dblquad
result2, error2 = integrate.dblquad(lambda y, x: x**2 + y**2, 
                                    0, 1, 
                                    lambda x: 0, lambda x: 1)
print(f"dblquad结果: {result2:.6f}")

案例7：三重积分

# 计算 ∫∫∫(x² + y² + z²) dxdydz
def f3(z, y, x):
    return x**2 + y**2 + z**2
result, error = integrate.nquad(f3, [[0, 1], [0, 1], [0, 1]])
print(f"三重积分结果: {result:.6f}")

实际应用案例

案例8：计算概率分布

# 计算正态分布的概率
from scipy.stats import norm
# P(0 ≤ X ≤ 1.96)，其中X ~ N(0,1)
def normal_pdf(x):
    return norm.pdf(x)
result, error = integrate.quad(normal_pdf, 0, 1.96)
print(f"P(0 ≤ X ≤ 1.96) ≈ {result:.4f}")
print(f"验证: norm.cdf(1.96) - norm.cdf(0) = {norm.cdf(1.96) - norm.cdf(0):.4f}")

案例9：曲线长度计算

# 计算抛物线 y = x² 从 x=0 到 x=2 的曲线长度
# 公式: L = ∫√(1 + (dy/dx)²) dx
def curve_length(x):
    dy_dx = 2*x  # y = x² 的导数
    return np.sqrt(1 + dy_dx**2)
length, error = integrate.quad(curve_length, 0, 2)
print(f"曲线 y = x² 从 x=0 到 x=2 的长度: {length:.4f}")

案例10：物理应用-功的计算

# 计算变力做功: F(x) = 10x², 从 x=0 到 x=3
def force(x):
    return 10 * x**2
work, error = integrate.quad(force, 0, 3)
print(f"变力做功: {work:.2f} 焦耳")

数值积分方法对比

# 对比不同积分方法
import time
def test_function(x):
    return np.sin(np.exp(x))
# 方法1：quad（自适应）
start = time.time()
result1, error1 = integrate.quad(test_function, 0, 2)
time1 = time.time() - start
# 方法2：Romberg积分
start = time.time()
result2 = integrate.romberg(test_function, 0, 2)
time2 = time.time() - start
print(f"quad: 结果={result1:.6f}, 误差={error1:.2e}, 时间={time1:.4f}s")
print(f"romberg: 结果={result2:.6f}, 时间={time2:.4f}s")

带奇异点的积分处理

# 处理带奇异点的积分
def integrand_with_singularity(x):
    # 在x=0有奇点
    return 1 / (np.sqrt(x) * (1 + x))
# 使用weight参数处理
result, error = integrate.quad(integrand_with_singularity, 
                                0, 1, 
                                points=[0])  # 指定奇点位置
print(f"带奇点积分结果: {result:.6f}")

常用技巧

# 1. 提高精度
result, error = integrate.quad(np.sin, 0, np.pi, 
                               epsabs=1e-12,  # 绝对误差
                               epsrel=1e-12)  # 相对误差
# 2. 设置限制
result, error = integrate.quad(np.sin, 0, np.pi,
                               limit=1000)  # 最大子区间数
# 3. 处理振荡函数
def oscillatory_func(x):
    return np.sin(100*x) * np.exp(-x)
result, error = integrate.quad(oscillatory_func, 0, np.inf,
                               limit=1000)

这些案例涵盖了SciPy积分计算的主要应用场景,根据你的具体需求选择合适的方法：