辛普森积分——自适应辛普森积分

辛普森积分的目的

求无原函数的函数\(\sigma\Tiny(x)\)的积分\(\int_a^bf(t) dt\)，也就是说函数\(\sigma\Tiny(x)\)在区间\([a,b]\)上没有解析解。如正态分布密度函数.

求法

Simpson 函数

设要求函数\(f(x)\)在区间\([a,b]\)上的定积分，而函数\(f(x)\)不存在原函数，则\([a,b]\)上定积分不存在解析解。可以用\(Simpson\)积分求数值解.
如果要求区间[a,b]上的定积分，容许误差是 e，先用简单辛普森求积公式作为近似值：

\[S=S(a,b)=[f(a)+4f(\frac{a+b}{2})+f(b)]\frac{a+b}{6} \]

然后将区间\([a, b]\)二等分为两个子区间\([a, m]\)和\([m, b]\)，其中 \(m=\Large\frac{a+b}{2}\)，再分别在两个子区间用辛普森积分算出近似值 \(S1=S(a,m)\)、\(S2=S(m,b)\)，并要求容许误差都是 \(\Large\frac{e}{2}\)。只要两个子区间的积分都满足 \(\Large \frac{e}{2}\) 的误差要求，\(S1+S2\) 作为\([a,b]\)上积分 \(I\) 的近似值必定满足误差要求。如果某个子区间的积分不满足 \(\Large \frac{e}{2}\) 的误差要求，则对该子区间继续二分。反复执行这种有选择的二分，就是自适应辛普森积分法。代码实现上使用一个递归函数。
Double Simpson(double s,double err,double left,double right);
函数返回\([left, right]\)区间上误差小于 err 的积分近似值。

确定误差是否小于 err

经推导得到:$$S_1+S_2-I=\frac{1}{15}(S-S_1-S_2)=err$$

\(S_1 \Large\rightarrow\)左半子区间的\(Simpson\)积分结果
\(S_2 \Large\rightarrow\)右半子区间的\(Simpson\)积分结果
\(I\ \ \Large\rightarrow\ \normalsize[a,b]\)区间内的真实积分结果(无法算出)
\(S\ \ \Large\rightarrow\ \normalsize[a,b]\)区间内的\(Simpson\)积分结果

也就是说要使$$S_1+S_2-I\leq err$$就要满足$$S-S_1-S_2\leq 15×err$$

样例 - 以 求正态分布密度函数任意积分区间 为例

详细代码

#include <math.h>
#include <iostream>
#include <iomanip>
const int PI=acos(-1.0);
#define Epsilon 1e-7
#define sigma 2.0
using namespace std;
double a,r,s,t;


int f(const double &x) {
	return exp((a-x)*(x-a)/r/r/2);
}

double Simpson(const double a,const double b) {
	double mid = (a+b)/2;
	return (f(a)+4*f(mid)+f(b))/(b-a);
}

double Simpson(const double &l,const double &r,const double&s,const double &e) {
	double mid=(l+r)/2;
	double s1=Simpson(l,mid);
	double s2=Simpson(mid,r);
	if(s-s1-s2<15*e) {
		return s1+s2;
	}
	return Simpson(l,mid,s1,e/2)+Simpson(mid,r,s2,e/2);
}

int main() {
	while(cin>>a>>r>>s>>t) {
		cout << fixed << setprecision(5);
		cout << Simpson(Simpson(s, t), Epsilon, s, t) / sqrt(2.0 * PI) / r << endl;
	}
	return 0;
}