积分

参考资料

• 数值积分 - OI Wiki
• 辛普森積分法 - 维基百科

辛普森积分法

$$\int_a^b f(x)\mathrm{d}x\approx\frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$$

例题

洛谷 P4525 【模板】自适应辛普森法 1

试计算积分：

$$\displaystyle{\int_L^R\frac{cx+d}{ax+b}\mathrm{d}x}$$

结果保留至小数点后 $6$ 位。

代码（2）

牛顿莱布尼茨公式

#include <bits/stdc++.h>
using namespace std;

double a,b,c,d,l,r;
double f1(double x)
{
	return c/a*x+(a*d-b*c)/(a*a)*log(fabs(a*x+b));
}
double f2(double x)
{
	return c/(b*2)*x*x+d/b*x;
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	cin>>a>>b>>c>>d>>l>>r;
	cout<<fixed<<setprecision(6)<<(a!=0?f1(r)-f1(l):f2(r)-f2(l))<<'\n';
	return 0;
}

自适应辛普森法

#include <bits/stdc++.h>
#define mid ((l+r)/2)
using namespace std;

const double eps=1e-8;
double a,b,c,d;
double F(double x)
{
	return (c*x+d)/(a*x+b);
}
double simp(double l,double r)
{
	return (F(l)+F(r)+F(mid)*4)*(r-l)/6;
}
double asr(double l,double r,double e,double ans)
{
	double L=simp(l,mid),R=simp(mid,r);
	if(fabs(L+R-ans)<=e*15)return L+R+(L+R-ans)/15;
	return asr(l,mid,e/2,L)+asr(mid,r,e/2,R);
}
double calc(double l,double r,double e)
{
	return asr(l,r,e,simp(l,r));
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	double l,r;
	cin>>a>>b>>c>>d>>l>>r;
	cout<<fixed<<setprecision(6)<<calc(l,r,eps)<<'\n';
	return 0;
}

洛谷 P4526 【模板】自适应辛普森法 2

给定实数 $a$，计算积分 $\displaystyle{\int_0^\infty x^{\frac{a}{x}-x}\mathrm{d}x}$。若积分发散，请输出 orz。（$|a|\le 50$）

代码（1）

#include <bits/stdc++.h>
#define mid ((l+r)/2)
using namespace std;

const double eps=1e-8;
double a;
double F(double x)
{
	return pow(x,a/x-x);
}
double simp(double l,double r)
{
	return (F(l)+F(r)+F(mid)*4)*(r-l)/6;
}
double asr(double l,double r,double e,double ans)
{
	double L=simp(l,mid),R=simp(mid,r);
	if(fabs(L+R-ans)<=e*15)return L+R+(L+R-ans)/15;
	return asr(l,mid,e/2,L)+asr(mid,r,e/2,R);
}
double calc(double l,double r,double e)
{
	return asr(l,r,e,simp(l,r));
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	cin>>a;
	if(a<0)cout<<"orz"<<'\n';
	else cout<<fixed<<setprecision(5)<<calc(eps,20,eps)<<'\n';
	return 0;
}