跳到正文
格致开物MATHWIKI

数值积分

数值积分(numerical integration)用有限个函数值或观测值近似计算定积分。最常见的思路是把曲线下的区域分成小段,用矩形、梯形或容易积分的多项式代替,再将各段结果相加。

只用几个数,怎样估计曲线下的面积

考虑 I=01x2dx。它的精确值为三分之一,适合用来检查近似方法。先把区间分成 [0,1/2][1/2,1] 两段,每段宽为 h=1/2

如果每段用中点的函数值作矩形高度,就在 1/43/4 处采样: M2=12[(14)2+(34)2]=12(116+916)=516. 这里字母 M 表示中点公式,下标二表示分成两段。结果约为 0.3125,比三分之一略小。

另一种方法连接各段两个端点,形成梯形。两个梯形面积分别为 12(0+1/4)/212(1/4+1)/2,相加得 T2=12[0+12+14]=38. 结果为 0.375,比真值略大。函数 x2 向上弯,端点连线在曲线上方,这解释了梯形高估的方向。

同一平方函数的中点矩形、两段梯形与辛普森抛物线比较,分别得到十六分之五、八分之三、三分之一
对同一积分使用不同局部近似。三幅图横纵坐标相同,面积可以直接比较。

左图每段矩形的顶边经过曲线中点;中图每条直线连接曲线端点;右图使用通过零、二分之一、一这三个位置的抛物线,它正好就是原函数。右图对应的辛普森公式将在下面推出。

从两小段写成一般公式

[a,b] 等分为 n 段,令 h=(ba)/n,节点为 xi=a+ih。复合中点公式为 Mn=hi=0n1f(a+(i+12)h). 每项是一个中点矩形的面积。复合梯形公式为 Tn=h[f(a)+f(b)2+i=1n1f(xi)]. 内部节点是左右两个梯形共用的端点,各贡献一半,合计权重一;最外两端只有一个梯形,因此权重为一半。

若被积函数表示流量,单位为升每秒,那么 h 以秒计,乘积单位是升,得到总流量。函数为负时也保留符号,例如速度积分表示净位移,而非全部路程。求积规则计算的是这种带符号的累积量。

误差为什么与步长平方有关

在一小段 [u,v] 内,以直线 L(x) 连接两个端点的函数值。如果 f 二次连续可微,线性插值的余项为 f(x)L(x)=f(ξx)2(xu)(xv), 其中 ξx[u,v]。设整个积分区间上 |f|M2,对两者之差积分并取上界: |uv(fL)dx|M22uv(xu)(vx)dx. 令段长为 h=vu,再令 t=xu,右侧的小积分成为 0ht(ht)dt=h32h33=h36. 所以每段误差不超过 M2h3/12。一共有 n=(ba)/h 段,将这些上界相加,得到 |ITn|(ba)M2h212. 局部误差是三次方,累加的段数与 1/h 成正比,因此总误差成为二次方。

中点公式也可用泰勒展开估计。以中点为中心,线性项在左右两侧积分相消;二阶余项的绝对值不超过 M2t2/2。积分 t[h/2,h/2] 后,每段上界为 M2h3/24,于是 |IMn|(ba)M2h224. 开头例子中 f=2h=1/2,两个界分别为 1/241/48,恰好等于实际误差。OpenStax:数值积分与误差界

辛普森的一、四、一权重从何而来

在对称区间 [h,h] 上,用二次多项式 p(t)=A+Bt+Ct2 近似函数。它的精确积分是 hhp(t)dt=2hA+23h3C, 因为奇函数项 Bt 的积分为零。另一方面,三个节点的加权和满足 h3[p(h)+4p(0)+p(h)]=h3[6A+2Ch2]=2hA+23h3C. 因此权重一、四、一准确积分任意二次多项式。对三次项 Dt3,积分与对称节点的加权和也都为零,所以公式还能准确处理三次多项式。

将每两小段作为一组,段数 n 取正偶数,就得到复合辛普森公式: Sn=h3[f(x0)+f(xn)+4j=0n/21f(x2j+1)+2j=1n/21f(x2j)]. 每组中心的权重为四,两组共用端点的权重合为二,最外端点权重为一。对于 x2 和两小段,结果为 S2=16[0+4×14+1]=13. 这就是图中右侧精确贴合的原因。

在四次连续可微且 |f(4)|M4 时,辛普森的误差保证为 |ISn|(ba)M4h4/180。这个标准余项界需要四阶光滑性,证明与完整条件见上述教材。它比中点、梯形使用了更强的函数信息。

计算对数积分,并提前确定需要的段数

考虑 I=12x1dx=log2。取四段,h=1/4,节点 1,5/4,3/2,7/4,2 的函数值分别为 1,4/5,2/3,4/7,1/2。按辛普森权重相加: S4=112[1+445+223+447+12]=174725200.6932539683. 四阶导数为 24/x5,区间上不超过二十四,所以 |IS4|2418044=119200.00052083.log20.6931471806 比较,实际误差约为 0.00010679,确实在保证范围内。

若事先要求误差不超过百万分之一,则应选偶数 n 使 24/(180n4)106。取 n=20 即可,误差界约为 8.33×107。这种精度规划使用已知导数界,不需要先知道准确积分值。

加密哪些地方,以及何时停下

均匀划分在每处投入相同采样量。若函数只在局部变化剧烈,可以只细分那些区间,这称为自适应求积。常见做法比较一个辛普森面板的粗结果 Sh 与细分后的 Sh/2。如果已进入误差主要为 Ch4 的范围,细误差约为粗误差的十六分之一,两结果的差约为细误差的十五倍,所以用 |Sh/2Sh|/15 估计细结果误差。它是依赖这一误差模型的估计,强度不同于已知导数界给出的保证。

程序要把总容差分配到各子区间,并累加误差估计。如果达到最大深度、求值次数或浮点分辨率,应返回停止状态;积分接近零时还需要正的绝对容差。梯形网格减半时可以复用旧节点,只补新中点,减少函数求值成本。FNC:Numerical integration

曲线有尖点,或者只有观测表格时

|x| 在跨越零的区间积分,先在零处分段,每段都是线性的,梯形公式便精确。对 1/x 在零附近积分,不能直接在零采样;可作 x=t2 代换,把 x1/2dx 变成 2dt,先消除端点奇性再求积。

规则采样还可能看不见振荡。例如在 [0,1]n 等分端点采样 sin2(πnx),每个值都是零,梯形结果为零,而真实积分为二分之一。这是特意与网格对齐的反例,说明有限样本在没有更多函数信息时不能保证捕捉所有变化。

若输入是不等时刻 ti 的观测值 yi,分段线性插值给出 Q=iti+1ti2(yi+yi+1). 每段必须用自己的时间宽度。插值生成的新点不等于新增真实观测。若各样本误差绝对值不超过 η,正权重求积的数据误差至多为 ηiwi=(ba)η,这部分不会因插值加密而消失。

能够自由选择节点时,还可使用高斯求积等方法提高多项式精确次数;只能使用固定观测表时,选择受到已有数据限制。函数结构、节点来源和误差目标共同决定适合的公式。

历史

求面积的近似方法早于微积分,插值与微积分的发展进一步给出了系统的求积公式。十八世纪的托马斯·辛普森与一、四、一规则的传播相关,这一规则本身有更早的先例。MacTutor:Thomas Simpson 现代数值积分在这些公式上进一步发展了节点选择、自适应细分和浮点误差控制。

参考资料与知识联系