辛普森积分学习笔记
定积分
定积分的定义
设函数 \(f(x)\) 在区间 \([a,b]\) 上有界，在 \([a,b]\) 中插入若干个分点
\[a=x_0<x_1<x_2<\cdots<x_{n-1}<x_n=b \]把区间 \([a,b]\) 分成 \(n\) 个小区间，各小区间的长度依次为 \(\Delta x_i=x_i-x_{i-1}(1\le i\le n)\)。在各小区间上任取一点 \(\xi_i\)（\(\xi\in[x_{i-1},x_i]\)），作乘积 \(f(\xi_i)\Delta x_i\) 并求和
\[S=\sum_{i=1}^{n}f(\xi_i)\Delta_i \]记 \(\lambda=\max(\Delta x_1,\Delta x_2,\cdots,\Delta x_n)\)，如果不论对 \([a,b]\) 怎样分，也不论小区间 \([x_{i-1},x_i]\) 上点 \(\xi_i\) 怎样的取法，只要当 \(\lambda\to 0\) 时，和 \(S\) 总趋于确定的极限 \(I\)，我们称这个极限 \(I\) 为函数 \(f(x)\) 在区间 \([a,b]\) 上的定积分，记作 \(\int_{a}^{b}f(x)dx\)。即
\[\int_{a}^{b}f(x)dx=I=\lim_{\lambda\to 0}\sum_{i=1}^{n}f(\xi_i)\Delta x_i \]其中 \(f(x)\) 叫做被积函数，\(f(x)dx\) 叫做被积表达式， \(x\) 叫做积分变量，\(a\) 叫做积分下限，\(b\) 叫做积分上限，\([a,b]\) 叫做积分区间。积分值仅与被积函数与积分区间有关，与积分变量用什么字母无关，即
\[\int_{a}^{b}f(x)dx=\int_{a}^{b}f(t)dt=\int_{a}^{b}f(u)du \]同时定义
- \(\int_{a}^{a}f(x)dx=0\)
- \(\int_{a}^{b}f(x)dx=-\int_{b}^{a}f(x)dx\)
定积分的几何意义即曲线 \(y=f(x)\) 和直线 \(x=a,x=b\) 以及 \(x\) 轴围成各部分中，\(x\) 轴上方图形面积之和减去 \(x\) 轴下方图形面积之和。
可积的定义
当函数 \(f(x)\) 在区间 \([a,b]\) 上的定积分存在时，称 \(f(x)\) 在区间 \([a,b]\) 上可积。若 \(f(x)\) 在区间 \([a,b]\) 上可积，则 \(f(x)\) 在 \([a,b],(a,b],[a,b),(a,b)\) 上均可积且定积分相等。
充分条件：当函数 \(f(x)\) 在区间 \([a,b]\) 上连续时，则 \(f(x)\) 在区间 \([a,b]\) 上可积；当函数 \(f(x)\) 在区间 \([a,b]\) 上有界，且只有有限个间断点，则 \(f(x)\) 在区间 \([a,b]\) 上可积。
必要条件：若函数 \(f(x)\) 在区间 \([a,b]\) 上可积，则 \(f(x)\) 在 \([a,b]\) 上必有界。
性质
- 定积分对函数有线性性，即若 \(\alpha,\beta\) 为常数，则
- 对区间有可加性，即对于任意三个常数 \(a,b,c\)，都有
- 同时具有保号性，即对于 \(a<b\)，如果在区间 \([a,b]\) 上 \(f(x)\ge 0\)，则
- 如果 \(a<b\) 且 \(M,m\) 分别是函数 \(f(x)\) 在 \([a,b]\) 上的最大值即最小值，则
- 设 \(f(x)\) 在 \([a,b]\) 上连续，则至少存在一点 \(\xi\in[a,b]\) 使得
其中 \(f(\xi)=\frac{1}{b-a}\int_{a}^{b}f(x)dx\) 称为函数 \(f(x)\) 在区间 \([a,b]\) 上的平均值。
辛普森积分
对于朴素的积分计算，有牛顿-莱布尼茨公式，即
\[\int_{a}^{b}f(x)dx=F(b)-F(a) \]其中 \(F(x)\) 是 \(f(x)\) 的原函数，即 \(f(x)=F'(x)\)。
但对于没有原函数或不好求解原函数的函数，我们可以采取梯形法则，即考虑用 \([a,b]\) 两点连线得到的直线代替曲线 \(y=f(x)\) 从而得到一个近似值，也叫做一阶牛顿-柯特斯闭型积分公式
\[I=\int_{a}^{b}f(x)dx\cong \int_{a}^{b}f_n(x)dx \]其中 \(f_n(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}+a_nx^n\)。一阶指的就是 \(n=1\) 的情况。而一个 \(n\) 次多项式的原函数是好求的，\(f_n(x)\) 对应的原函数即为 \(F_n(x)=C+a_0x+\frac{a_1}{2}x^2+\cdots+\frac{a_{n-1}}{n}x^n+\frac{a_n}{n+1}x^{n+1}\)，其中 \(C\) 为常数。那么在这个近似于原函数的积分下，我们可以得到原函数积分的近似值。
而辛普森积分法就是二阶牛顿-柯特斯闭型积分公式，即利用一个二次函数去适配原函数，但众所周知，三个点才能确定一个二次函数，那么我们需要在区间内再取一个点，我们取 \(c=\frac{a+b}{2}\) 作为第三个点，那么有
\[\int_{a}^{b}f(x)dx\cong\int_{a}^{b}f_2(x)dx=\frac{b-a}{6}[f(a)+4f(c)+f(b)] \]下面给出证明：
设由 \((a,f(a)),(b,f(b)),(c,f(c))\) 三点得到的二次函数为 \(f_2(x)=Ax^2+Bx+C\)，所以 \(F_2(x)=\frac{A}{3}x^3+\frac{B}{2}x^2+Cx+\alpha\)，由牛顿-莱布尼茨公式得
\[\begin{aligned}\int_{a}^{b}f(x)dx&\cong\int_{a}^{b}f_2(x)dx\\&=F_2(b)-F_2(a)\\&=\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(b-a)\\&=\frac{b-a}{6}[2A(a^2+ab+b^2)+3B(a+b)+6C]\\&=\frac{b-a}{6}[(Aa^2+Ba+C)+4(A(\frac{a+b}{2})^2+B(\frac{a+b}{2})+C)+(Ab^2+Bb+C)]\\&=\frac{b-a}{6}[(Aa^2+Ba+C)+4(Ac^2+Bc+C)+(Ab^2+Bb+C)]\\&=\frac{b-a}{6}[f(a)+4f(c)+f(b)] \end{aligned} \]那么，在已知 \(f(x)\) 的情况下我们可以得到原函数定积分的近似值。
自适应辛普森积分法
但是，我们考虑对于一个不解近二次函数的函数图像，进行这样的逼近还是过于粗糙了。为了避免误差过大，也是为了利用计算机实现，我们需要学习自适应辛普森积分法。
为了减少误差，我们显然可以将区间不断细化，这样一定可以逼近准确值，但是，计算机的运行速度限制了我们，我们不能将原区间拆分得过于细致，并且我们不能接受较大的误差，那么我们不妨对于区间 \([a,b]\) 先行给出一个估算值 \(A\)，接着对区间 \([a,c],[c,b]\) 分别算出其估算值 \(A_1,A_2\)。由积分的区间可加性，那么 \(I\cong A_1+A_2\cong A\)，只要 \(A\) 与 \(A_1+A_2\) 的差值不至于过大，那么我们就可以认为当前区间的求解已经足够。反之则利用区间的可加性具体求解两侧区间的积分后相加得到答案。
但是，不难考虑到如果 \(A,A_1,A_2\) 的误差均较大，但其差值恰好较小，我们依旧会认为其合理，为了避免这种情况，我们可以人为设定一个递归层数，只有超过这个层数我们才会考虑是否直接退出的问题。同时为了防止小误差相加得到大误差，我们在递归时需要将误差逐次减小。
并且我们发现，我们可能会重复调用许多次 f 函数，因此可以采用记忆化的方式进行卡常。
代码
double f(double x){
//利用高超的算法求解已知的 f
}
double Simpson(double l,double r){
double mid=(l+r)/2;
return (f(l)+f(r)+4*f(mid))*(r-l)/6;//辛普森积分
}
double asr(double l,double r,double ans,double eps){
double mid=(l+r)/2;
double l_=Simpson(l,mid),r_=Simpson(mid,r);//估算左右两个区间的积分
if(fabs(l_+r_-ans)<=15*eps)return l_+r_+(l_+r_-ans)/15;//判断误差
return asr(l,mid,l_,eps/2)+asr(mid,r,r_,eps/2);//递归求解
//事实上，这里的 15,2 都是可以自己改变的参数，根据经验和题目限制选择合适的参数即可
}
double asr(double l,double r,double eps){
return asr(l,r,Simpson(l,r),eps);
}
int main(){
double ans=asr(L,R,1e-6);
//这里的 L,R 是积分区间，但是对于一个收敛于 0 的函数可以自己二分出一个合理的上限避免积分区间过大
//对代码的时间复杂度造成影响的因素有：区间 [L,R]，误差，误差参数。
return 0;
}