2.1 数值积分与微分方程求解

2.1 数值积分与微分方程求解

本节摘要：数值积分是用有限个采样点估计定积分的技术。scipy.integrate 的 quad 基于自适应 Gauss-Kronrod 求积，自动在函数变化剧烈处加密采样并返回误差估计，是一维定积分的首选工具；dblquad、tplquad 把同一套思想扩展到二重与三重积分；solve_ivp 求解常微分方程初值问题，内置 RK45 等变步长方法并支持事件检测。本节从能跑的最小例子开始，逐个讲清用法与取舍。

核心问题

阅读完本节，你应当能够：

1. 用 quad 计算一维定积分，并正确解读返回值中的误差估计；
2. 面对奇异点、振荡被积函数，说出该直接算、拆区间还是换工具；
3. 用 dblquad、tplquad 把二重、三重积分区域写成嵌套边界；
4. 用 solve_ivp 求解常微分方程初值问题，设置精度参数与事件函数；
5. 说出数值积分与数值求解微分方程在思想上的共同点。

一、问题与直觉

先看一个真实场景。你要估算一条河流某断面的年输沙量，手里只有每隔一段时间的含沙量观测值，曲线下的面积就是总量。手算的办法是画网格数格子，或者把曲线切成一条条细长条求和——这就是数值积分最初的直觉：用有限个点去近似无限细的分割。

问题在于"切成多细才够"。切粗了，函数在条内剧烈起伏时误差大得离谱；切细了，平坦区域纯属浪费计算量。更麻烦的是，我们事先并不知道函数哪里起伏剧烈。第 1 章提过，手写 Simpson 公式只要十几行代码，但真碰上带奇异点或高频振荡的被积函数，朴素实现会安静地给出错误答案，连个警告都没有。

quad 的思路是让计算机自己决定在哪里加密。它先在一个子区间上同时用低阶和高阶公式各算一遍，两个结果的差就是该子区间误差的估计；误差超标就二分这个子区间，递归处理，直到每个子区间都满足容差。这种"哪里误差大就加密哪里"的策略叫自适应求积，是数值积分从"会算"走向"算得准"的关键一步。它和后面要讲的变步长微分方程求解是同一棵树上结的两个果子。

二、核心原理

2.1 quad：自适应求积的最小例子

import numpy as np
from scipy.integrate import quad

def signal(x):
    return np.exp(-0.5 * x) * np.sin(6 * x)

area, err = quad(signal, 0, 4)
print(f"积分值: {area:.6f}，误差估计: {err:.2e}")
​

两行调用，返回两个值：积分值和绝对误差估计。这个例子的被积函数是指数衰减乘以高频振荡，固定步长方法要踩很多点才能跟住振荡；quad 内部用的是 Gauss-Kronrod 嵌套公式——同一个区间上用两套不同阶数的节点同时求积，高阶结果作参考值，两套结果的差当作误差，差太大就切半区间重来，直到每个子区间都达标。

一个常被忽略的细节：第二个返回值不是"保证上界"，而是算法对误差的估计。对光滑函数它通常相当保守，对振荡剧烈或接近奇异的情况可能偏乐观。工程上的习惯是把误差值打印出来看数量级，而不是默认它可信。quad 还能接受 epsabs 与 epsrel 两个容差参数，默认都在 1e-8 附近，需要更快或更准时先动它们。

2.2 奇异点、振荡与无穷区间

数值积分最怕两件事：区间内有不可积奇点，以及函数剧烈振荡。先看端点奇点：

area, err = quad(lambda x: 1 / np.sqrt(x), 0, 1)
print(f"奇异积分值: {area:.6f}，误差: {err:.2e}")
​

被积函数在零点趋于无穷，但积分本身收敛，quad 直接就能处理。如果奇点出现在区间内部，比如被积函数含 x 减 0.5 的绝对值的倒数，那就需要把区间在奇点处拆成两段，让每个子区间只含端点奇点，再分别求积相加。

振荡函数是另一类麻烦。quad 对有限次振荡还能应付，振荡频率高到一定程度会触发"积分可能发散或收敛缓慢"的警告。这时候先别急着加大 limit 参数放宽递归层数，更有效的做法是把区间拆成若干段，每段振荡次数有限，再逐段求积求和。无穷区间也不是问题：quad 支持把积分上限写成正无穷，它内部通过变量代换把无穷区间映射到有限区间，精度照样有保障。

area, err = quad(lambda x: np.exp(-x), 0, np.inf)
print(f"无穷区间积分值: {area:.6f}，误差: {err:.2e}")
​

无穷区间的实用场景很多：概率密度函数的期望与方差、衰减过程的累积量，都涉及从某点到无穷的积分。指数衰减这类函数收敛很快，几个积分点就够；但如果被积函数在无穷远处只是缓慢衰减，比如 1 除以 x 的平方，收敛会慢得多，这时先把积分拆成"有限段加尾巴段"分别求，反而更快更稳。

2.3 dblquad 与 tplquad：嵌套求积

二重积分不过是把一维求积套两层。dblquad 的签名里，被积函数先收 y 再收 x，y 的上下界各是一个函数（可以是 x 的函数）。比如计算单位圆第一象限上函数 x 平方加 y 平方的积分：

from scipy.integrate import dblquad

def integrand(y, x):
    return x**2 + y**2

val, err = dblquad(integrand, 0, 1, 0, lambda x: np.sqrt(1 - x**2))
print(f"二重积分值: {val:.6f}，误差: {err:.2e}")
​

三重积分 tplquad 同理，多一层嵌套，内层边界依次是外层变量的函数。写法上最容易翻车的是参数顺序：内层变量写在最前面，边界函数按从外到内的顺序给出。顺序写反程序不报错，只给你一个看起来合理的错数，这类 bug 最难抓。核对签名永远第一步，写完后用对称性、量纲或已知特例验算一遍是第二步。

嵌套积分在真实问题里出现得比想象中频繁：概率论里二维分布的归一化常数、物理里刚体的转动惯量、流体力学里截面的流量，全是二重或三重积分。它们的共同套路是先把积分区域画出来，确定"外层变量扫什么范围、内层边界怎么随外层变化"，再照抄进代码。区域画不清楚就写代码，等于闭着眼睛开车。

2.4 solve_ivp：常微分方程初值问题

积分解决"总量"，微分方程解决"演化"。很多物理过程写成方程就是：状态对时间的导数等于某个函数，再加上初始时刻的状态，构成初值问题。solve_ivp 的默认方法是 RK45——四阶 Runge-Kutta 配五阶误差估计，同样是"算两步比较误差、自动变步长"的套路，和自适应求积一脉相承。

from scipy.integrate import solve_ivp

def pendulum(t, s):
    theta, omega = s
    return [omega, -np.sin(theta)]

sol = solve_ivp(pendulum, [0, 10], [1.2, 0.0], rtol=1e-8, atol=1e-10)
print(sol.success, sol.nfev)
print(sol.t[-1], sol.y[0, -1])
​

注意返回对象的几个字段：success 表示求解是否成功，nfev 是函数求值次数，t 和 y 是输出网格与对应状态。精度控制靠 rtol 与 atol 这对相对、绝对容差，而不是步长——步长由算法自己决定，你只表达"我要多准"。想在任意时刻取值，打开 dense_output 获得连续解，供后续插值使用。

事件检测是 solve_ivp 的杀手锏。想记录摆锤第一次摆回竖直位置的时刻，把"角度等于零"写成事件函数，再把 direction 设为 1 表示只记向上穿越：

def cross_zero(t, s):
    return s[0]

cross_zero.direction = 1

sol = solve_ivp(pendulum, [0, 10], [1.2, 0.0], events=(cross_zero,), rtol=1e-8)
print("第一次过零时刻:", sol.t_events[0][0])
​

事件函数返回零即触发，算法在步进过程中用根查找精确定位穿越点，不会因为步长过大而错过事件。这是手写欧拉法完全做不到的——自己写循环得时刻盯着符号变化再二分，远不如把事件函数交给求解器省心。轨道交汇、碰撞时刻、阈值报警，这类问题用事件检测都能干净利落地解决。

三、工程实践要点

3.1 什么时候别用 quad

quad 不是万能的。被积函数只有离散采样点、没有连续表达式时，应该用 numpy 的 trapezoid 或 scipy.integrate 的 simpson 做固定网格求积——它们直接吃数组，不需要函数对象。维度很高时（比如超过四重），嵌套求积的代价指数爆炸，应当考虑蒙特卡洛方法。下表给出快速选型。

3.2 精度与性能的平衡

quad 的 epsabs 与 epsrel 分别控制绝对与相对容差，默认约 1e-8。对精度要求不高的场景，把容差放宽到 1e-6 往往能让计算快好几倍；反过来，被积函数量级很小或很大时，要同时调绝对容差，否则相对误差会被量级吃光。limit 参数限制递归深度，振荡函数频繁报警告时，先想到拆区间而不是无脑加大 limit。微分方程求解也类似：先默认容差跑通，再按物理需求收紧，一上来就追求 1e-12 只会让求解器在数值噪声里空转。

⚠️ 常见坑：dblquad 与 tplquad 的参数顺序是"内层变量在前，边界函数从外向内"。写反了程序不报错，只输出一个看着合理其实错误的数。写完务必用已知特例验算，比如把区域换成矩形、函数换成常数，心算结果一对照就露馅。
💡 关键直觉：自适应求积和变步长求解微分方程是同一套思想——用低阶与高阶结果的差估计误差，误差大就加密，误差小就大步走。理解这一个机制，quad、solve_ivp 乃至第 3 节优化的收敛判定都能串起来，学一个顶三个。

3.3 初值问题的工程习惯

给 solve_ivp 设置容差时，先默认 rtol 1e-8 跑通，再按需收紧；需要固定输出时刻时用 t_eval 参数，别依赖算法自己的网格。事件函数尽量写成连续、光滑的表达式，direction 方向信息能显著减少无关触发。对刚性方程——化学反应、电路这类时间尺度相差悬殊的系统——默认 RK45 会因稳定性限制把步长压到无法接受，此时把 method 换成 Radau 或 LSODA 往往立竿见影。判断刚性的实用信号：RK45 步长小得离谱、nfev 爆炸，但换隐式方法后瞬间变快，那就是典型的刚性症状。

重点提炼

• 要点一：quad 基于自适应 Gauss-Kronrod 嵌套求积，自动加密采样并返回误差估计。
• 要点二：误差估计是参考值不是保证上界，打印出来看数量级是基本习惯。
• 要点三：端点奇点交给 quad 处理，区间内部奇点先拆区间再分别求积。
• 要点四：dblquad 与 tplquad 的边界函数按从外到内的顺序嵌套，参数顺序写反不报错。
• 要点五：solve_ivp 默认 RK45，精度靠 rtol 与 atol 控制，步长由算法自适应决定。
• 要点六：事件函数返回零触发，配 direction 精确捕获轨迹穿越时刻。
• 要点七：刚性方程要换 Radau 或 LSODA，硬用 RK45 会慢到怀疑人生。
• 要点八：只有离散采样点或维度很高的积分，别用 quad，换固定网格求积或蒙特卡洛。

下一节我们从"算面积"转向"找极值"：优化与根查找。拟合参数、解方程、求最小值，本质是同一类搜索问题，只是搜索策略不同——选对策略，正是下一节的主线。