数值积分

数学 / 数值算法

本地源文件:docs/math__numerical__integral.md

数值积分

定积分的定义

简单来说,函数 𝑓(𝑥)f(x) 在区间 [𝑙,𝑟][l,r] 上的定积分 ∫𝑟𝑙𝑓(𝑥)d𝑥∫lrf(x)dx 指的是 𝑓(𝑥)f(x) 在区间 [𝑙,𝑟][l,r] 中与 𝑥x 轴围成的区域的面积(其中 𝑥x 轴上方的部分为正值,𝑥x 轴下方的部分为负值).

很多情况下,我们需要高效,准确地求出一个积分的近似值.下面介绍的 辛普森法 ,就是这样一种求数值积分的方法.

辛普森法

这个方法的思想是将被积区间分为若干小段,每段套用二次函数的积分公式进行计算.

二次函数积分公式(辛普森公式)

对于一个二次函数 𝑓(𝑥) =𝑎𝑥2 +𝑏𝑥 +𝑐f(x)=ax2+bx+c,有:

∫𝑟𝑙𝑓(𝑥)d𝑥=(𝑟−𝑙)(𝑓(𝑙)+𝑓(𝑟)+4𝑓(𝑙+𝑟2))6∫lrf(x)dx=(r−l)(f(l)+f(r)+4f(l+r2))6

推导过程: 对于一个二次函数 𝑓(𝑥) =𝑎𝑥2 +𝑏𝑥 +𝑐f(x)=ax2+bx+c; 求积分可得 𝐹(𝑥) =∫𝑥0𝑓(𝑥)d𝑥 =𝑎3𝑥3 +𝑏2𝑥2 +𝑐𝑥 +𝐷F(x)=∫0xf(x)dx=a3x3+b2x2+cx+D 在这里 D 是一个常数,那么

∫𝑟𝑙𝑓(𝑥)d𝑥=𝐹(𝑟)−𝐹(𝑙)=𝑎3(𝑟3−𝑙3)+𝑏2(𝑟2−𝑙2)+𝑐(𝑟−𝑙)=(𝑟−𝑙)(𝑎3(𝑙2+𝑟2+𝑙𝑟)+𝑏2(𝑙+𝑟)+𝑐)=𝑟−𝑙6(2𝑎𝑙2+2𝑎𝑟2+2𝑎𝑙𝑟+3𝑏𝑙+3𝑏𝑟+6𝑐)=𝑟−𝑙6((𝑎𝑙2+𝑏𝑙+𝑐)+(𝑎𝑟2+𝑏𝑟+𝑐)+4(𝑎(𝑙+𝑟2)2+𝑏(𝑙+𝑟2)+𝑐))=𝑟−𝑙6(𝑓(𝑙)+𝑓(𝑟)+4𝑓(𝑙+𝑟2))∫lrf(x)dx=F(r)−F(l)=a3(r3−l3)+b2(r2−l2)+c(r−l)=(r−l)(a3(l2+r2+lr)+b2(l+r)+c)=r−l6(2al2+2ar2+2alr+3bl+3br+6c)=r−l6((al2+bl+c)+(ar2+br+c)+4(a(l+r2)2+b(l+r2)+c))=r−l6(f(l)+f(r)+4f(l+r2))

根据这个辛普森公式,我们先介绍一种普通的辛普森积分法.

普通辛普森法

1743 年,这种方法发表于托马斯·辛普森的一篇论文中.

描述

给定一个自然数 𝑛n,将区间 [𝑙,𝑟][l,r] 分成 2𝑛2n 个等长的区间 𝑥x

𝑥𝑖 =𝑙 +𝑖ℎ, 𝑖 =0…2𝑛,xi=l+ih, i=0…2n, ℎ =𝑟−𝑙2𝑛.h=r−l2n.

我们就可以计算每个小区间 [𝑥2𝑖−2,𝑥2𝑖][x2i−2,x2i],𝑖 =1…𝑛i=1…n 的积分值,将所有区间的积分值相加即为总积分.

对于 [𝑥2𝑖−2,𝑥2𝑖][x2i−2,x2i],𝑖 =1…𝑛i=1…n 的一个区间,选其中的三个点 (𝑥2𝑖−2,𝑥2𝑖−1,𝑥2𝑖)(x2i−2,x2i−1,x2i) 就可以构成一条抛物线从而得到一个函数 𝑃(𝑥)P(x),这个函数存在且唯一.计算原函数在该区间的积分值就变成了计算新的二次函数 𝑃(𝑥)P(x) 在该段区间的积分值.这样我们就可以利用辛普森公式来近似计算它.

∫𝑥2𝑖𝑥2𝑖−2𝑓(𝑥) 𝑑𝑥 ≈∫𝑥2𝑖𝑥2𝑖−2𝑃(𝑥) 𝑑𝑥 =(𝑓(𝑥2𝑖−2)+4𝑓(𝑥2𝑖−1)+(𝑓(𝑥2𝑖))ℎ3∫x2i−2x2if(x) dx≈∫x2i−2x2iP(x) dx=(f(x2i−2)+4f(x2i−1)+(f(x2i))h3

将其分段求和即可得到如下结论:

∫𝑟𝑙𝑓(𝑥)𝑑𝑥 ≈(𝑓(𝑥0)+4𝑓(𝑥1)+2𝑓(𝑥2)+4𝑓(𝑥3)+2𝑓(𝑥4)+…+4𝑓(𝑥2𝑁−1)+𝑓(𝑥2𝑁))ℎ3∫lrf(x)dx≈(f(x0)+4f(x1)+2f(x2)+4f(x3)+2f(x4)+…+4f(x2N−1)+f(x2N))h3

误差

我们直接给出结论,普通辛普森法的误差为:

−190(𝑟−𝑙2)5𝑓(4)(𝜉)−190(r−l2)5f(4)(ξ)

其中 𝜉ξ 是位于区间 [𝑙,𝑟][l,r] 的某个值.

实现

C++Python

---|---

---|---

自适应辛普森法

普通的方法为保证精度在时间方面无疑会受到 𝑛n 的限制,我们应该找一种更加合适的方法.

现在唯一的问题就是如何进行分段.如果段数少了计算误差就大,段数多了时间效率又会低.我们需要找到一个准确度和效率的平衡点.

我们这样考虑:假如有一段图像已经很接近二次函数的话,直接带入公式求积分,得到的值精度就很高了,不需要再继续分割这一段了.

于是我们有了这样一种分割方法:每次判断当前段和二次函数的相似程度,如果足够相似的话就直接代入公式计算,否则将当前段分割成左右两段递归求解.

现在就剩下一个问题了:如果判断每一段和二次函数是否相似?

我们把当前段直接代入公式求积分,再将当前段从中点分割成两段,把这两段再直接代入公式求积分.如果当前段的积分和分割成两段后的积分之和相差很小的话,就可以认为当前段和二次函数很相似了,不用再递归分割了.

上面就是自适应辛普森法的思想.在分治判断的时候,除了判断精度是否正确,一般还要强制执行最少的迭代次数.

参考代码如下:

C++Python

---|---

---|---

习题

参考资料

<https://doi.org/10.1145/321526.321537>:该文章讨论了自适应 Simpson 法的改进方案,其中详细论述了上文代码中的常数 15 的由来与优势.

本页面最近更新: 2026/1/7 08:56:54,更新历史 发现错误?想一起完善?在 GitHub 上编辑此页! 本页面贡献者:Ir1d, H-J-Granger, StudyingFather, Enter-tainer, Tiphereth-A, countercurrent-time, NachtgeistW, ksyx, ShaoChenHeng, AngelKitty, CCXXXI, Chrogeek, cjsoft, diauweb, Early0v0, ezoixx130, GekkaSaori, Great-designer, iamtwz, Konano, LovelyBuggies, Makkiy, mao1t, Menci, mgt, minghu6, Nanarikom, P-Y-Y, PotassiumWings, SamZhangQingChuan, sshwy, SukkaW, Suyun514, weiyong1024, alphagocc, GavinZhengOI, Gesrua, kxccc, lychees, Peanut-Tang, shawlleyw, Xeonacid, Yanjun-Zhao, zyj-111 本页面的全部内容在CC BY-SA 4.0SATA 协议之条款下提供,附加条款亦可能应用