算法学习：插值型求积公式

牛顿-柯斯特(Newton-Cotes)求积公式

定义

牛顿-柯斯特(Newton-Cotes)求积公式是插值型求积公式的特殊形式
在插值求积公式

∫baf(x)dx≈∫baP(x)dx=∑k=0nAkf(xk)∫abf(x)dx≈∫abP(x)dx=∑k=0nAkf(xk)\int_a^bf(x)dx \approx \int_a^bP(x)dx=\sum\limits_{k=0}^n A_kf(x_k)

中所取节点是等距时称为牛顿-柯斯特公式
其中插值多项式

P(x)=∑k=0nℓk(x)f(xk)P(x)=∑k=0nℓk(x)f(xk)P(x)=\sum\limits_{k=0}^n\ell_k(x)f(x_k)

求积系数

Ak=∫baℓk(x)dxAk=∫abℓk(x)dxA_k=\int_a^b\ell_k(x)dx

这里ℓk(x)ℓk(x)\ell_k(x)指的是插值基函数，即有

Ak=∫baℓk(x)dx=∫ba∏j=0,j≠knx−xjxk−xjdxAk=∫abℓk(x)dx=∫ab∏j=0,j≠knx−xjxk−xjdxA_k=\int_a^b\ell_k(x)dx=\int_a^b\prod\limits_{j=0,j\neq k}^n\frac{x-x_j}{x_k-x_j}dx

推导

在[a,b][a,b][a,b]区间内设置等距的插值基点a=x0<x1⋯<xn=ba=x0<x1⋯<xn=ba=x_0，设节点步长为h=b−an,xk=a+kh,k∈{0,1,⋯,n}h=b−an,xk=a+kh,k∈{0,1,⋯,n}h=\frac{b-a}{n},x_k=a+kh,k\in \{0,1,\cdots ,n\}
积分作变量替换x=a+thx=a+thx=a+th
xk−xj=(k−j)h,x−xj=(t−j)h,dx=hdtxk−xj=(k−j)h,x−xj=(t−j)h,dx=hdtx_k-x_j=(k-j)h,x-x_j=(t-j)h,dx=hdt
(xk−x0)⋯(xk−xk−1)(xk−xk+1)⋯(xk−xn)(xk−x0)⋯(xk−xk−1)(xk−xk+1)⋯(xk−xn)(x_k-x_0)\cdots (x_k-x_{k-1})(x_k-x_{k+1})\cdots(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))hn=(−1)n−kk!(n−k)!hn=k(k-1)\cdots1(-1)(-2)\cdots(-(n-k))h^n=(-1)^{n-k}k!(n-k)!h^n
(x−x0)⋯(x−xk−1)(x−xk+1)⋯(x−xn)(x−x0)⋯(x−xk−1)(x−xk+1)⋯(x−xn)(x-x_0)\cdots (x-x_{k-1})(x-x_{k+1})\cdots(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)hn=t(t-1)\cdots(t-k+1)(t-k-1)\cdots(t-n)h^n
当x=ax=ax=a时，t=0t=0t=0，当x=bx=bx=b时，t=nt=nt=n
于是有Ak=∫baℓk(x)dx=h∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)(−1)n−kk!(n−k)!dtAk=∫abℓk(x)dx=h∫0nt(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)(−1)n−kk!(n−k)!dtA_k=\int_a^b\ell_k(x)dx=h\int_0^n\frac{t(t-1)\cdots(t-k+1)(t-k-1)\cdots(t-n)}{(-1)^{n-k}k!(n-k)!}dt
=h⋅(−1)n−kk!(n−k)!∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt=h⋅(−1)n−kk!(n−k)!∫0nt(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt=h\cdot \frac{(-1)^{n-k}}{k!(n-k)!} \int_0^nt(t-1)\cdots(t-k+1)(t-k-1)\cdots(t-n)dt
=nh⋅(−1)n−kk!(n−k)!n∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt=nh⋅(−1)n−kk!(n−k)!n∫0nt(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt=nh\cdot \frac{(-1)^{n-k}}{k!(n-k)!n}\int_0^nt(t-1)\cdots(t-k+1)(t-k-1)\cdots(t-n)dt
=(b−a)⋅(−1)n−kk!(n−k)!n∫n0t(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt=(b−a)⋅(−1)n−kk!(n−k)!n∫0nt(t−1)⋯(t−k+1)(t−k−1)⋯(t−n)dt=(b-a)\cdot \frac{(-1)^{n-k}}{k!(n-k)!n}\int_0^nt(t-1)\cdots(t-k+1)(t-k-1)\cdots(t-n)dt
=(b−a)⋅(−1)n−kk!(n−k)!n∫n0∏j=0,j≠kn(t−j)dt=(b−a)⋅(−1)n−kk!(n−k)!n∫0n∏j=0,j≠kn(t−j)dt=(b-a)\cdot \frac{(-1)^{n-k}}{k!(n-k)!n}\int_0^n\prod\limits_{j=0,j\neq k}^n(t-j) dt
不难发现，式子中的(−1)n−kk!(n−k)!n∫n0∏j=0,j≠kn(t−j)dt(−1)n−kk!(n−k)!n∫0n∏j=0,j≠kn(t−j)dt \frac{(-1)^{n-k}}{k!(n-k)!n}\int_0^n\prod\limits_{j=0,j\neq k}^n(t-j) dt与步长h无关，所以可以预处理。我们将上述表达式称之为柯斯特求积系数，记为

c(n)k=(−1)n−kk!(n−k)!n∫n0∏j=0,j≠kn(t−j)dtck(n)=(−1)n−kk!(n−k)!n∫0n∏j=0,j≠kn(t−j)dtc_k^{(n)}= \frac{(-1)^{n-k}}{k!(n-k)!n}\int_0^n\prod\limits_{j=0,j\neq k}^n(t-j) dt

那么我们的求积系数

Ak=(b−a)c(n)kAk=(b−a)ck(n)A_k=(b-a)c_k^{(n)}

于是求得牛顿-柯斯特公式

∫baf(x)dx≈(b−a)∑k=0nc(n)kf(xk)∫abf(x)dx≈(b−a)∑k=0nck(n)f(xk)\int_a^bf(x)dx \approx(b-a)\sum\limits_{k=0}^n c_k^{(n)}f(x_k)

~~上面你就当我装了一波逼，其实闷声背公式也是可以滴~~

公式应用：梯形公式

梯形公式

我们将上述式子令n=1，得到柯斯特系数
c(1)0=−∫10(t−1)dt=12c0(1)=−∫01(t−1)dt=12c_0^{(1)}=-\int_0^1(t-1)dt=\frac{1}{2}
c(1)1=−∫10tdt=12c1(1)=−∫01tdt=12c_1^{(1)}=-\int_0^1tdt=\frac{1}{2}
于是得到梯形公式

I1(f)=(b−a)12f(a)+f(b−a)12f(b)I1(f)=(b−a)12f(a)+f(b−a)12f(b)I_1(f)=(b-a)\frac{1}{2}f(a)+f(b-a)\frac{1}{2}f(b)

梯形公式具有一次代数精确度~~（就是没什么卵用）~~

辛普森积分

辛普森积分在立体几何中的公式应用

在立体几何中，辛普森积分用来计算拟柱体体积公式

公式内容

设拟柱体的高（两底面α,βα,β\alpha,\beta间的距离）为H，如果用平行于底面的平面γ去截该图形，所得到的截面面积是平面γγ\gamma与平面αα\alpha之间距离h的不超过3次的函数，那么该拟柱体的体积V为

V=H(S1+4S2+S3)6V=H(S1+4S2+S3)6V=\frac{H(S_1+4S_2+S_3)}{6}

式中，S1S1S_1和S2S2S_2是两底面的面积，S0S0S_0是中截面的面积

公式证明

V=∫f(x)dxV=∫f(x)dxV=\int f(x)dx
=ah44+bh33+ch22+dh=ah44+bh33+ch22+dh=a\frac{h^4}{4}+b\frac{h^3}{3}+c\frac{h^2}{2}+dh
利用公式计算体积，可以得到：
V=H(S1+4S2+S3)6V=H(S1+4S2+S3)6V=\frac{H(S_1+4S_2+S_3)}{6}
=h(f(0)+4f(h2)+f(h))/6=h(f(0)+4f(h2)+f(h))/6=h ( f(0) + 4f(\frac{h}{2}) + f(h) ) /6
=16h[d+4(ah38+bh24+ch2+d)+(ah3+bh2+ch+d)]=16h[d+4(ah38+bh24+ch2+d)+(ah3+bh2+ch+d)]=\frac{1}{6}h [ d + 4 (a\frac{h^3} {8} + b\frac{h^2} {4} + c\frac{h}{2} + d) + (ah^3 + bh^2 + ch + d) ]
=ah44+bh33+ch22+dh=ah44+bh33+ch22+dh=a\frac{h^4}{4}+b\frac{h^3}{3}+c\frac{h^2}{2}+dh
于是原命题得证

辛普森积分拟合目标函数

根据牛顿-柯斯特求积公式，令n=2，带入得到柯斯特系数
c(2)0=14∫20(t−1)(t−2)dt=16c0(2)=14∫02(t−1)(t−2)dt=16c_0^{(2)}=\frac{1}{4}\int_0^2(t-1)(t-2)dt=\frac{1}{6}
c(2)1=−12∫20t(t−2)dt=46c1(2)=−12∫02t(t−2)dt=46c_1^{(2)}=-\frac{1}{2}\int_0^2t(t-2)dt=\frac{4}{6}
c(2)2=14∫20t(t−1)dt=16c2(2)=14∫02t(t−1)dt=16c_2^{(2)}=\frac{1}{4}\int_0^2t(t-1)dt=\frac{1}{6}
于是有

I2(f)=(b−a)(16f(a)+46f(a+b2)+16f(b))I2(f)=(b−a)(16f(a)+46f(a+b2)+16f(b))I_2(f)=(b-a)(\frac{1}{6}f(a)+\frac{4}{6}f(\frac{a+b}{2})+\frac{1}{6}f(b))
=16(b−a)(f(a)+4f(a+b2)+f(b))=16(b−a)(f(a)+4f(a+b2)+f(b))=\frac{1}{6}(b-a)(f(a)+4f(\frac{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)ans=∫xminxmaxf(x)ans=\int_{xmin}^{xmax}f(x)
函数无法直接求出，于是我们采用自适应辛普森积分法即可。
(ps：这题卡精度，小心为妙)

代码

/**************************************************************Problem: 2178User: 2014lvzelongLanguage: C++Result: AcceptedTime:4412 msMemory: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
*/

小结

对于插值型求积公式，在信息学中主要应用自适应辛普森积分算法。此算法的优点在于好写易懂，并有一定的精确度。但是对于某些“二分精度”题，还是小心为妙。