算法学习：插值型求积公式
牛顿-柯斯特(Newton-Cotes)求积公式
定义
牛顿-柯斯特(Newton-Cotes)求积公式是插值型求积公式的特殊形式
在插值求积公式
∫baf(x)dx≈∫baP(x)dx=∑k=0nAkf(xk)
∫
a
b
f
(
x
)
d
x
≈
∫
a
b
P
(
x
)
d
x
=
∑
k
=
0
n
A
k
f
(
x
k
)
中所取节点是等距时称为牛顿-柯斯特公式
其中插值多项式
P(x)=∑k=0nℓk(x)f(xk)
P
(
x
)
=
∑
k
=
0
n
ℓ
k
(
x
)
f
(
x
k
)
求积系数
Ak=∫baℓk(x)dx
A
k
=
∫
a
b
ℓ
k
(
x
)
d
x
这里
ℓk(x)
ℓ
k
(
x
)
指的是插值基函数，即有
Ak=∫baℓk(x)dx=∫ba∏j=0,j≠knx−xjxk−xjdx
A
k
=
∫
a
b
ℓ
k
(
x
)
d
x
=
∫
a
b
∏
j
=
0
,
j
≠
k
n
x
−
x
j
x
k
−
x
j
d
x
推导
在
[a,b]
[
a
,
b
]
区间内设置等距的插值基点
a=x0<x1⋯<xn=b
a
=
x
0
<
x
1
⋯
<
x
n
=
b
，设节点步长为
h=b−an,xk=a+kh,k∈{0,1,⋯,n}
h
=
b
−
a
n
,
x
k
=
a
+
k
h
,
k
∈
{
0
,
1
,
⋯
,
n
}
积分作变量替换
x=a+th
x
=
a
+
t
h
xk−xj=(k−j)h,x−xj=(t−j)h,dx=hdt
x
k
−
x
j
=
(
k
−
j
)
h
,
x
−
x
j
=
(
t
−
j
)
h
,
d
x
=
h
d
t
(xk−x0)⋯(xk−xk−1)(xk−xk+1)⋯(xk−xn)
(
x
k
−
x
0
)
⋯
(
x
k
−
x
k
−
1
)
(
x
k
−
x
k
+
1
)
⋯
(
x
k
−
x
n
)
=k(k−1)⋯1(−1)(−2)⋯(−(n−k))hn=(−1)n−kk!(n−k)!hn
=
k
(
k
−
1
)
⋯
1
(
−
1
)
(
−
2
)
⋯
(
−
(
n
−
k
)
)
h
n
=
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
h
n
(x−x0)⋯(x−xk−1)(x−xk+1)⋯(x−xn)
(
x
−
x
0
)
⋯
(
x
−
x
k
−
1
)
(
x
−
x
k
+
1
)
⋯
(
x
−
x
n
)
=t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)hn
=
t
(
t
−
1
)
⋯
(
t
−
k
+
1
)
(
t
−
k
−
1
)
⋯
(
t
−
n
)
h
n
当
x=a
x
=
a
时，
t=0
t
=
0
，当
x=b
x
=
b
时，
t=n
t
=
n
于是有
Ak=∫baℓk(x)dx=h∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)(−1)n−kk!(n−k)!dt
A
k
=
∫
a
b
ℓ
k
(
x
)
d
x
=
h
∫
0
n
t
(
t
−
1
)
⋯
(
t
−
k
+
1
)
(
t
−
k
−
1
)
⋯
(
t
−
n
)
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
d
t
=h⋅(−1)n−kk!(n−k)!∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt
=
h
⋅
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
∫
0
n
t
(
t
−
1
)
⋯
(
t
−
k
+
1
)
(
t
−
k
−
1
)
⋯
(
t
−
n
)
d
t
=nh⋅(−1)n−kk!(n−k)!n∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt
=
n
h
⋅
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
n
∫
0
n
t
(
t
−
1
)
⋯
(
t
−
k
+
1
)
(
t
−
k
−
1
)
⋯
(
t
−
n
)
d
t
=(b−a)⋅(−1)n−kk!(n−k)!n∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt
=
(
b
−
a
)
⋅
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
n
∫
0
n
t
(
t
−
1
)
⋯
(
t
−
k
+
1
)
(
t
−
k
−
1
)
⋯
(
t
−
n
)
d
t
=(b−a)⋅(−1)n−kk!(n−k)!n∫n0∏j=0,j≠kn(t−j)dt
=
(
b
−
a
)
⋅
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
n
∫
0
n
∏
j
=
0
,
j
≠
k
n
(
t
−
j
)
d
t
不难发现，式子中的
(−1)n−kk!(n−k)!n∫n0∏j=0,j≠kn(t−j)dt
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
n
∫
0
n
∏
j
=
0
,
j
≠
k
n
(
t
−
j
)
d
t
与步长h无关，所以可以预处理。我们将上述表达式称之为柯斯特求积系数，记为
c(n)k=(−1)n−kk!(n−k)!n∫n0∏j=0,j≠kn(t−j)dt
c
k
(
n
)
=
(
−
1
)
n
−
k
k
!
(
n
−
k
)
!
n
∫
0
n
∏
j
=
0
,
j
≠
k
n
(
t
−
j
)
d
t
那么我们的求积系数
Ak=(b−a)c(n)k
A
k
=
(
b
−
a
)
c
k
(
n
)
于是求得牛顿-柯斯特公式
∫baf(x)dx≈(b−a)∑k=0nc(n)kf(xk)
∫
a
b
f
(
x
)
d
x
≈
(
b
−
a
)
∑
k
=
0
n
c
k
(
n
)
f
(
x
k
)
上面你就当我装了一波逼，其实闷声背公式也是可以滴
公式应用：梯形公式
梯形公式
我们将上述式子令n=1，得到柯斯特系数
c(1)0=−∫10(t−1)dt=12
c
0
(
1
)
=
−
∫
0
1
(
t
−
1
)
d
t
=
1
2
c(1)1=−∫10tdt=12
c
1
(
1
)
=
−
∫
0
1
t
d
t
=
1
2
于是得到梯形公式
I1(f)=(b−a)12f(a)+f(b−a)12f(b)
I
1
(
f
)
=
(
b
−
a
)
1
2
f
(
a
)
+
f
(
b
−
a
)
1
2
f
(
b
)
梯形公式具有一次代数精确度（就是没什么卵用）
辛普森积分
辛普森积分在立体几何中的公式应用
在立体几何中，辛普森积分用来计算拟柱体体积公式
公式内容
设拟柱体的高（两底面
α,β
α
,
β
间的距离）为H，如果用平行于底面的平面γ去截该图形，所得到的截面面积是平面
γ
γ
与平面
α
α
之间距离h的不超过3次的函数，那么该拟柱体的体积V为
V=H(S1+4S2+S3)6
V
=
H
(
S
1
+
4
S
2
+
S
3
)
6
式中，
S1
S
1
和
S2
S
2
是两底面的面积，
S0
S
0
是中截面的面积
公式证明
V=∫f(x)dx
V
=
∫
f
(
x
)
d
x
=ah44+bh33+ch22+dh
=
a
h
4
4
+
b
h
3
3
+
c
h
2
2
+
d
h
利用公式计算体积，可以得到：
V=H(S1+4S2+S3)6
V
=
H
(
S
1
+
4
S
2
+
S
3
)
6
=h(f(0)+4f(h2)+f(h))/6
=
h
(
f
(
0
)
+
4
f
(
h
2
)
+
f
(
h
)
)
/
6
=16h[d+4(ah38+bh24+ch2+d)+(ah3+bh2+ch+d)]
=
1
6
h
[
d
+
4
(
a
h
3
8
+
b
h
2
4
+
c
h
2
+
d
)
+
(
a
h
3
+
b
h
2
+
c
h
+
d
)
]
=ah44+bh33+ch22+dh
=
a
h
4
4
+
b
h
3
3
+
c
h
2
2
+
d
h
于是原命题得证
辛普森积分拟合目标函数
根据牛顿-柯斯特求积公式，令n=2，带入得到柯斯特系数
c(2)0=14∫20(t−1)(t−2)dt=16
c
0
(
2
)
=
1
4
∫
0
2
(
t
−
1
)
(
t
−
2
)
d
t
=
1
6
c(2)1=−12∫20t(t−2)dt=46
c
1
(
2
)
=
−
1
2
∫
0
2
t
(
t
−
2
)
d
t
=
4
6
c(2)2=14∫20t(t−1)dt=16
c
2
(
2
)
=
1
4
∫
0
2
t
(
t
−
1
)
d
t
=
1
6
于是有
I2(f)=(b−a)(16f(a)+46f(a+b2)+16f(b))
I
2
(
f
)
=
(
b
−
a
)
(
1
6
f
(
a
)
+
4
6
f
(
a
+
b
2
)
+
1
6
f
(
b
)
)
=16(b−a)(f(a)+4f(a+b2)+f(b))
=
1
6
(
b
−
a
)
(
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
)
这就是辛普森积分公式，说明了辛普森积分公式就是牛顿-柯斯特求积公式的特殊情况。
辛普森积分具有良好的稳定性与收敛性，在数值逼近领域有很重要的应用。
自适应辛普森积分算法
引例
bzoj2178圆的面积并
算法流程
辛普森积分不论具有多么优秀的稳定性与收敛性，在大数据范围内仍会出现较大误差。因此自适应辛普森积分算法是辛普森积分的衍生算法。其基本步骤如下
- 对于区间[l,r]上的函数求出其辛普森积分的答案ans
- 分别求出区间[l,m],[m,r]两个子区间的辛普森积分的答案之和ret
- 若ans和ret的误差范围不超过eps，返回ret的值，否则递归左右子区间。
在对于圆滑曲线的近似积分中，自适应辛普森积分具有非常优秀的表现。
引例分析
对于平面上的圆的面积并，我们首先排除内含的所有圆，然后建立函数模型f(a)表示在圆所覆盖直线x=a的线段长度。那么就有
ans=∫xmaxxminf(x)
a
n
s
=
∫
x
m
i
n
x
m
a
x
f
(
x
)
函数无法直接求出，于是我们采用自适应辛普森积分法即可。
(ps：这题卡精度，小心为妙)
代码
/**************************************************************
Problem: 2178
User: 2014lvzelong
Language: C++
Result: Accepted
Time:4412 ms
Memory:1340 kb
****************************************************************/
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<cmath>
#include<string>
using namespace std;
const int N = 1101;
const double eps = 1e-6;
int read() {
char ch = getchar(); int x = 0, f = 1;
while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
while(ch >= '0' && ch <= '9') {x = (x << 1) + (x << 3) - '0' + ch; ch = getchar();}
return x * f;
}
struct C{
int x, y, r;
C(int xx = 0, int yy = 0, int rr = 0) : x(xx), y(yy), r(rr) {}
C operator - (C &a) {return C(x - a.x, y - a.y);}
int operator ^ (C &a) {return x * a.x + y * a.y;}
}c[N];
int n, st, ed, tot, L, R;
bool del[N];
struct Li{double l, r;}p[N];
int sqr(double x) {return x * x;}
int sqr(C a) {return a ^ a;}
bool cmp1(C a, C b) {return a.r < b.r;}
bool cmp2(C a, C b) {return a.x - a.r < b.x - b.r;}
bool cmp3(Li a, Li b) {return a.l < b.l;}
void Del() {
sort(c + 1, c + n + 1, cmp1);
for(int i = 1;i <= n; ++i)
for(int j = i + 1;j <= n; ++j)
if(sqr(c[i] - c[j]) <= sqr(c[j].r - c[i].r))
{del[i] = true; break;}
int top = 0;
for(int i = 1;i <= n; ++i) if(!del[i]) c[++top] = c[i];
n = top;
sort(c + 1, c + n + 1, cmp2);
}
double F(double x) {
double ret = 0; tot = 0;
for(int i = 1;i <= n; ++i)
if(fabs(c[i].x - x) <= c[i].r - eps) {
double dis = sqrt(sqr(c[i].r) - (x - c[i].x) * (x - c[i].x));
p[++tot].l = c[i].y - dis; p[tot].r = c[i].y + dis;
}
if(!tot) return 0;
sort(p + 1, p + tot + 1, cmp3);
double h = -2200;
for(int i = 1;i <= tot; ++i) {
if(p[i].l > h) ret += p[i].r - p[i].l, h = p[i].r;
else if(p[i].r > h) ret += p[i].r - h, h = p[i].r;
}
return ret;
}
double calc(double l, double r) {return (r - l) / 6 * (F(l) + 4 * F((l + r) / 2.0) + F(r)) ;}
double Simpson(double l, double r) {
double m = (l + r) / 2.0, s1 = calc(l, m) + calc(m, r), s2 = calc(l, r);
if(fabs(s1 - s2) < eps) return s1;
else return Simpson(l, m) + Simpson(m, r);
}
int main() {
n = read();
for(int i = 1;i <= n; ++i) {
c[i].x = read(); c[i].y = read(); c[i].r = read();
L = min(L, c[i].x - c[i].r); R = max(R, c[i].x + c[i].r);
}
Del();
printf("%.3lf\n", Simpson(-2000.0, 2000.0));
return 0;
}
/*
3
0 0 1
1 0 1
0 1 1
*/
小结
对于插值型求积公式，在信息学中主要应用自适应辛普森积分算法。此算法的优点在于好写易懂，并有一定的精确度。但是对于某些“二分精度”题，还是小心为妙。
