定积分的计算

定积分的计算有两种方法，利用原函数或者数值计算。但是原函数不一定求得出来，此时需要数值积分。数值积分有好几种算法，需要视具体情况选择精度与时间都能满足需要的算法。各种数值积分的算法可以参考相关的教材，这篇文章的目录也比较全，可以参考一下。一般而言，如果直接使用矩形求和（也称为黎曼和，对计算机而言相当于根据定义计算，只需将$dx$替换为$eps$即可）会超时，则需要考虑自适应辛普森积分。

辛普森积分公式与自适应辛普森积分

$${\int }_a^{b}f(x)dx\approx \frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b))$$
其来源是使用过三点的二次曲线去逼近源曲线，从而得到一个数值结果。更广义上辛普森公式是牛顿-柯斯特公式的一个特例。
自适应辛普森积分则是一种递归调用：对任意区间，将其分成2个子区间，分别使用辛普森积分再求和。如果子区间之和与全区间直接辛普森积分的差别在控制范围以内，则不需要再递归了。否则继续递归。

double simpson(double(*f)(double),double a,double b){double mid = ( a + b ) * 0.5;return (b-a)*(f(a)+f(b)+4.0*f(mid)) / 6.0;
}
//自适应辛普森积分的递归调用
double asr(double(*f)(double),double a,double b,double eps,double ans){double mid = ( a + b ) * 0.5;double lans = simpson(f,a,mid);double rans = simpson(f,mid,b);//精度够直接返回if(fabs(lans+rans-ans)<=eps) return ans;//精度不够递归调用return asr(f,a,mid,eps*0.5,lans) + asr(f,mid,b,eps*0.5,rans);
}
//自适应辛普森积分，区间[a, b]，精度eps，原函数f
double asr(double (*f)(double),double a,double b,double eps){return asr(f,a,b,eps,simpson(f,a,b));
}

洛谷自适应辛普森法1模板题，要求计算
$${\int }_L^R\frac{cx+d}{ax+b}dx$$
到小数点后6位。题目保证收敛。
直接计算的话，通过调整精度，可以得到一个几十分的代码。下面是一个60分的。

#include <bits/stdc++.h>
using namespace std;double const EPS = 6E-8;
double A,B,C,D,L,R;int main(){cin>>A>>B>>C>>D>>L>>R;double sum = 0;for(double x=L;x<=R;x+=EPS){sum += (C*x+D) / (A*x+B);}printf("%.6f\n",EPS * sum);return 0;
}

总之，精度和时间有矛盾，可以考虑使用自适应辛普森方法。

#include <bits/stdc++.h>
using namespace std;double simpson(double(*f)(double),double a,double b){double mid = ( a + b ) * 0.5;return (b-a)*(f(a)+f(b)+4.0*f(mid)) / 6.0;
}
//自适应辛普森积分的递归调用
double asr(double(*f)(double),double a,double b,double eps,double ans){double mid = ( a + b ) * 0.5;double lans = simpson(f,a,mid);double rans = simpson(f,mid,b);//精度够直接返回if(fabs(lans+rans-ans)<=eps) return ans;//精度不够递归调用return asr(f,a,mid,eps*0.5,lans) + asr(f,mid,b,eps*0.5,rans);
}
//自适应辛普森积分，区间[a, b]，精度eps，原函数f
double asr(double (*f)(double),double a,double b,double eps){return asr(f,a,b,eps,simpson(f,a,b));
}double a,b,c,d,L,R;
double f(double x){return (c*x+d)/(a*x+b);
}
int main(){scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);printf("%.6f\n",asr(f,L,R,1E-8));return 0;
}

洛谷自适应辛普森法2模板题要求计算如下积分：
$${\int }_0^{\infty }{x}^{\frac{a}{x}-x}dx$$
保留到小数点后5位，如果发散输出$orz$。可以证明$a<0$时积分发散。其次可知当$x$稍微大一点函数值几乎为0，因此积分区域不用取太大。

#include <bits/stdc++.h>
using namespace std;double const EPS = 1E-6;double simpson(double(*f)(double),double a,double b){double mid = ( a + b ) * 0.5;return (b-a)*(f(a)+f(b)+4.0*f(mid)) / 6.0;
}
//自适应辛普森积分的递归调用
double asr(double(*f)(double),double a,double b,double eps,double ans){double mid = ( a + b ) * 0.5;double lans = simpson(f,a,mid);double rans = simpson(f,mid,b);//精度够直接返回if(fabs(lans+rans-ans)<=eps) return ans;//精度不够递归调用return asr(f,a,mid,eps*0.5,lans) + asr(f,mid,b,eps*0.5,rans);
}
//自适应辛普森积分，区间[a, b]，精度eps，原函数f
double asr(double (*f)(double),double a,double b,double eps){return asr(f,a,b,eps,simpson(f,a,b));
}double a;
double f(double x){return pow(x,a/x-x);
}
int main(){scanf("%lf",&a);if(a<0)puts("orz");else printf("%.5f\n",asr(f,EPS,30.0,EPS));return 0;
}

龙贝格积分(Romberg积分)

龙贝格积分可以称得上是一种有效的加速计算积分的方法。首先引入一个记号$T(h)$表示积分区间长度为$h$，使用梯形去近似定积分。即：
$${\int }_L^{R}f(x)dx\approx \frac{h}{2}(f(L)+f(R))=T(h),其中h=R-L$$
此时$T$可以称为步长为$h$的一阶近似，因为其误差为$O(h)$。由于后续会采用2的幂次，因此使用脚标$0$来记录，记作：$T_0(h)$。
接着可以把全区间划分为更多的梯形，分别是$2,4,8,...$即$2$的幂序列，用这些梯形面积去近似定积分，于是得到了一个 $T$序列，分别是$T_0,T_1,T_2,T_3,...$
其中，
$$T_1=\frac{1}{2}{T}_0+\frac{h}{2}f(L+\frac{h}{2})$$
$$T_2=\frac{1}{2}{T}_1+\frac{h}{4}(f(L+\frac{h}{4})+f(L+\frac{3h}{4}))$$
$$T_n=\frac{1}{2}{T}_{n-1}+\frac{h}{2^n}\sum _{i=1,i取奇数}^{2^n}f(L+\frac{i}{2^n}h)$$
有了$T$序列，可以构造出辛普森序列，省略推导过程，直接给出结论：
$$S_n=\frac{1}{3}(4T_{n+1}-T_n)$$
再利用辛普森序列构造出Cotes序列：
$$C_n=\frac{1}{4^2-1}(4^2{S}_{n+1}-S_n)$$
最后利用$C$序列构造出Romberg序列：
$$R_n=\frac{1}{4^3-1}(4^3{C}_{n+1}-C_n)$$
计算$R_n$到满足一定精度，即可停止。
尽管数学过程并不一目了然，但是把$4$个序列写出来并观察其计算依赖关系，立刻可知，用$4$个一维数组或者一个二维数组依次计算即可。
$$\begin{array}{cccc}{T}_0& & & \\ T_1& S_0& & \\ T_2& S_1& C_0& \\ T_3& S_2& C_1& R_0\\ T_4& S_3& C_2& R_1\\ \cdots & & & \end{array}$$
仍然是洛谷自适应辛普森法1模板题，

#include <bits/stdc++.h>
using namespace std;double const EPS = 1E-8;
inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}//Romberg积分,求f(x)在[x1,x2]上的定积分
//r[0]表示T序列，r[1]表示S序列,r[2]表示C序列,r[3]表示R序列
double Romberg(double(*f)(double), double x1, double x2,double r[][4]){double h = x2 - x1;int m = 1;r[0][3] = 0.0;//首先计算T0r[0][0] = 0.5 * h * (f(x1)+f(x2));for(int i=1;;++i){double *now = *(r + i);double *prev = *(r + i - 1);double h2 = h;h *= 0.5, m <<= 1;//计算Tidouble x = x1 + h;double& sum = now[0];sum = 0.0;for(int j=1;j<m;j+=2){sum += f(x);x += h2;}sum = 0.5 * prev[0] + h * sum;//计算Sinow[1] = (4.0*now[0]-prev[0]) / 3.0;//计算Cinow[2] = (16.0*now[1]-prev[1]) / 15.0;//计算Rinow[3] = (64.0*now[2]-prev[2]) / 63.0;//cout<<i<<" "<<r[i][0]<<" "<<r[i][3]<<endl;//判断if(isEq(prev[3],now[3])) return now[3];}throw runtime_error("XXXXXX");
}//*/double T[10010][4];
double A,B,C,D,L,R;
double f(double x){return (C*x+D)/(A*x+B);
}
int main(){//freopen("1.txt","r",stdin);cin>>A>>B>>C>>D>>L>>R;printf("%.6f\n",Romberg(f,L,R,T));return 0;
}

洛谷自适应辛普森法2模板题用Romberg方法计算。

#include <bits/stdc++.h>
using namespace std;double const EPS = 1E-7;
double const PI = acos(-1.0);inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isLE(double x,double y){return sgn(x-y)<=0;};
inline bool isEq(double x,double y){return is0(x-y);}
inline bool isGE(double x,double y){return sgn(x-y)>=0;}//Romberg积分,求f(x)在[x1,x2]上的定积分
//r[0]表示T序列，r[1]表示S序列,r[2]表示C序列,r[3]表示R序列
double Romberg(double(*f)(double), double x1, double x2,double r[][4]){double h = x2 - x1;int m = 1;r[0][3] = 0.0;//首先计算T0r[0][0] = 0.5 * h * (f(x1)+f(x2));for(int i=1;;++i){double *now = *(r + i);double *prev = *(r + i - 1);double h2 = h;h *= 0.5, m <<= 1;//计算Tidouble x = x1 + h;double& sum = now[0];sum = 0.0;for(int j=1;j<m;j+=2){sum += f(x);x += h2;}sum = 0.5 * prev[0] + h * sum;//计算Sinow[1] = (4.0*now[0]-prev[0]) / 3.0;//计算Cinow[2] = (16.0*now[1]-prev[1]) / 15.0;//计算Rinow[3] = (64.0*now[2]-prev[2]) / 63.0;//cout<<i<<" "<<r[i][0]<<" "<<r[i][3]<<endl;//判断if(isEq(prev[3],now[3])) return now[3];}throw runtime_error("XXXXXX");
}//*/double a,Tmp[1010][4];
double f(double x){return pow(x,a/x-x);
}
int main(){scanf("%lf",&a);if(a<0)puts("orz");else printf("%.5f\n",Romberg(f,EPS,12.0,Tmp));return 0;
}

这个题目用Romberg方法时还要注意积分区域不能选的过大。因为如果选的过大的话，第一次计算时，起点、中点和终点的函数值恰好都是0（精度意义下），则直接退出循环，认为结果就是0。

面积或者体积的计算

hdu1724

hdu1724要求计算椭圆被两竖直线截取的面积。分别用自适应辛普森方法和龙贝格积分完成。

//自适应辛普森积分
#include <bits/stdc++.h>
using namespace std;
double const EPS = 1E-6;
inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}
double simpson(double(*f)(double),double a,double b){double mid = ( a + b ) * 0.5;return (b-a)*(f(a)+f(b)+4.0*f(mid)) / 6.0;
}
//自适应辛普森积分的递归调用
double asr(double(*f)(double),double a,double b,double eps,double ans){double mid = ( a + b ) * 0.5;double lans = simpson(f,a,mid);double rans = simpson(f,mid,b);//精度够直接返回if(fabs(lans+rans-ans)<=eps) return ans;//精度不够递归调用return asr(f,a,mid,eps*0.5,lans) + asr(f,mid,b,eps*0.5,rans);
}
//自适应辛普森积分，区间[a, b]，精度eps，原函数f
double asr(double (*f)(double),double a,double b,double eps){return asr(f,a,b,eps,simpson(f,a,b));
}
double A,B,L,R;double f(double x){return B * sqrt(1-x*x/(A*A));
}
int main(){//freopen("1.txt","r",stdin);int nofkase;scanf("%d",&nofkase);while(nofkase--){scanf("%lf%lf%lf%lf",&A,&B,&L,&R);printf("%.3f\n",asr(f,L,R,EPS)*2.0);}return 0;
}

//龙贝格积分
#include <bits/stdc++.h>
using namespace std;
double const EPS = 1E-5;
inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}//Romberg积分,求f(x)在[x1,x2]上的定积分
//r[0]表示T序列，r[1]表示S序列,r[2]表示C序列,r[3]表示R序列
double Romberg(double(*f)(double), double x1, double x2,double r[][4]){double h = x2 - x1;int m = 1;r[0][3] = 0.0;//首先计算T0r[0][0] = 0.5 * h * (f(x1)+f(x2));for(int i=1;;++i){double *now = *(r + i);double *prev = *(r + i - 1);double h2 = h;h *= 0.5, m <<= 1;//计算Tidouble x = x1 + h;double& sum = now[0];sum = 0.0;for(int j=1;j<m;j+=2){sum += f(x);x += h2;}sum = 0.5 * prev[0] + h * sum;//计算Sinow[1] = (4.0*now[0]-prev[0]) / 3.0;//计算Cinow[2] = (16.0*now[1]-prev[1]) / 15.0;//计算Rinow[3] = (64.0*now[2]-prev[2]) / 63.0;//cout<<i<<" "<<r[i][0]<<" "<<r[i][3]<<endl;//判断if(isEq(prev[3],now[3])) return now[3];}throw runtime_error("XXXXXX");
}//*/double T[22][4];
double A,B,L,R;double f(double x){return B * sqrt(1-x*x/(A*A));
}
int main(){//freopen("1.txt","r",stdin);int nofkase;scanf("%d",&nofkase);while(nofkase--){scanf("%lf%lf%lf%lf",&A,&B,&L,&R);printf("%.3f\n",Romberg(f,L,R,T)*2.0);}return 0;
}

NEERC 2014 D题 Damage Assessment