跳转至

# 数值积分

## 定积分的定义

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

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

## 辛普森法

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

二次函数积分公式（辛普森公式）

对于一个二次函数 𝑓(𝑥) =𝑎 𝑥 2 +𝑏 𝑥 +𝑐 f(x)=a x 2+b x+c，有：

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

推导过程： 对于一个二次函数 𝑓(𝑥) =𝑎 𝑥 2 +𝑏 𝑥 +𝑐 f(x)=a x 2+b x+c； 求积分可得 𝐹(𝑥) =∫𝑥 0 𝑓(𝑥)d 𝑥 =𝑎 3 𝑥 3 +𝑏 2 𝑥 2 +𝑐 𝑥 +𝐷 F(x)=∫0 x f(x)d x=a 3 x 3+b 2 x 2+c x+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))∫l r f(x)d x=F(r)−F(l)=a 3(r 3−l 3)+b 2(r 2−l 2)+c(r−l)=(r−l)(a 3(l 2+r 2+l r)+b 2(l+r)+c)=r−l 6(2 a l 2+2 a r 2+2 a l r+3 b l+3 b r+6 c)=r−l 6((a l 2+b l+c)+(a r 2+b r+c)+4(a(l+r 2)2+b(l+r 2)+c))=r−l 6(f(l)+f(r)+4 f(l+r 2))

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

### 普通辛普森法

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

#### 描述

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

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

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

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

∫𝑥 2 𝑖 𝑥 2 𝑖−2 𝑓(𝑥) 𝑑 𝑥 ≈∫𝑥 2 𝑖 𝑥 2 𝑖−2 𝑃(𝑥) 𝑑 𝑥 =(𝑓(𝑥 2 𝑖−2)+4 𝑓(𝑥 2 𝑖−1)+(𝑓(𝑥 2 𝑖))ℎ 3∫x 2 i−2 x 2 i f(x) d x≈∫x 2 i−2 x 2 i P(x) d x=(f(x 2 i−2)+4 f(x 2 i−1)+(f(x 2 i))h 3

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

∫𝑟 𝑙 𝑓(𝑥)𝑑 𝑥 ≈(𝑓(𝑥 0)+4 𝑓(𝑥 1)+2 𝑓(𝑥 2)+4 𝑓(𝑥 3)+2 𝑓(𝑥 4)+…+4 𝑓(𝑥 2 𝑁−1)+𝑓(𝑥 2 𝑁))ℎ 3∫l r f(x)d x≈(f(x 0)+4 f(x 1)+2 f(x 2)+4 f(x 3)+2 f(x 4)+…+4 f(x 2 N−1)+f(x 2 N))h 3

#### 误差

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

−1 9 0(𝑟−𝑙 2)5 𝑓(4)(𝜉)−1 90(r−l 2)5 f(4)(ξ)

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

#### 实现

C++Python

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    constexpr int N = 1000 * 1000;

    double simpson_integration(double a, double b) {
      double h = (b - a) / N;
      double s = f(a) + f(b);
      for (int i = 1; i <= N - 1; ++i) {
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
      }
      s *= h / 3;
      return s;
    }

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    N = 1000 * 1000

    def simpson_integration(a, b):
        h = (b - a) / N
        s = f(a) + f(b)
        for i in range(1, N):
            x = a + h * i
            if i & 1:
                s = s + f(x) * 4
            else:
                s = s + f(x) * 2
        s = s * (h / 3)
        return s

### 自适应辛普森法

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

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

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

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

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

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

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

参考代码如下：

C++Python

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    double simpson(double l, double r) {
      double mid = (l + r) / 2;
      return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;  // 辛普森公式
    }

    double asr(double l, double r, double eps, double ans, int step) {
      double mid = (l + r) / 2;
      double fl = simpson(l, mid), fr = simpson(mid, r);
      if (abs(fl + fr - ans) <= 15 * eps && step < 0)
        return fl + fr + (fl + fr - ans) / 15;  // 足够相似的话就直接返回
      return asr(l, mid, eps / 2, fl, step - 1) +
             asr(mid, r, eps / 2, fr, step - 1);  // 否则分割成两段递归求解
    }

    double calc(double l, double r, double eps) {
      return asr(l, r, eps, simpson(l, r), 12);
    }

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    def simpson(l, r):
        mid = (l + r) / 2
        return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6  # 辛普森公式

    def asr(l, r, eps, ans, step):
        mid = (l + r) / 2
        fl = simpson(l, mid)
        fr = simpson(mid, r)
        if abs(fl + fr - ans) <= 15 * eps and step < 0:
            return fl + fr + (fl + fr - ans) / 15  # 足够相似的话就直接返回
        return asr(l, mid, eps / 2, fl, step - 1) + asr(
            mid, r, eps / 2, fr, step - 1
        )  # 否则分割成两段递归求解

    def calc(l, r, eps):
        return asr(l, r, eps, simpson(l, r), 12)

## 习题

- Luogu4525【模板】自适应辛普森法 1
- HDU1724 Ellipse
- NOI2005 月下柠檬树

## 参考资料

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.0 和 SATA 协议之条款下提供，附加条款亦可能应用