一部分函数的积分是可以通过Newton-Leibniz formula计算的。 然而另一部分就不行，我们就要通过数值积分的方式计算。 我们首先看一下 拉格朗日插值法，有了函数的值，通过拉格朗日插值法就可以计算出函数式。 我们设每个多项式是pj(x)，我们通过构造，得出pj(x)的式子。我们要保证对于pj(x),任取一个i属于Ij，pj(xi)=0且pj(xj)=1; 最后我们就可以推出拉格朗日插值法的公式 我们继续推导 我们对两边积分。 我们让 让 那么我们就得到了数值积分的公式。 我们继续推导。 cni叫做柯特斯系数。 那么我们有 这是牛顿-柯特斯求积公式 我们取3个点x1=a；x2=a+b>>1;x3=b; 那么我们有 这就是辛普森公式。 hdu 1724
#include<cstdio> #include<cmath> using namespace std; typedef long long ll; const double eps=1e-6; double a,b,l,r; double f(double x) { return sqrt((b*b)*(1-(x*x)/(a*a))); } double simpson(double a,double b) { double c=a+(b-a)/2; return (f(a)+4*f(c)+f(b))*(b-a)/6; } double asr(double a,double b,double eps,double A) { double c=a+(b-a)/2; double L=simpson(a,c); double R=simpson(c,b); if(fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0; return asr(a,c,eps/2,L)+asr(c,b,eps/2,R); } double asr(double a,double b,double eps) { return asr(a,b,eps,simpson(a,b)); } int main() { int T; scanf("%d",&T); while(T--) { scanf("%lf%lf%lf%lf",&a,&b,&l,&r); printf("%.3f\n",2*asr(l,r,eps)); } return 0; }我们来计算一波圆周率。
#include<cstdio> #include<cmath> using namespace std; typedef long long ll; double a,b,l,r; double f(double x) { return sqrt(1-x*x); } double simpson(double a,double b) { double c=a+(b-a)/2; return (f(a)+4*f(c)+f(b))*(b-a)/6; } double asr(double a,double b,double eps,double A) { double c=a+(b-a)/2; double L=simpson(a,c); double R=simpson(c,b); if(fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0; return asr(a,c,eps/2,L)+asr(c,b,eps/2,R); } double asr(double a,double b,double eps) { return asr(a,b,eps,simpson(a,b)); } int main() { int T; printf("%.10f",asr(0,1,1e-11)*4); return 0; }