zoj 3898 Stean 高等数学 数值积分

题目

题目链接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3898

题目来源：ZOJ 月赛 给跪

简要题意：一个桶底和盖的$Z$坐标为$Z_1,Z_2$，桶的半径$R = 2 + \cos Z$，求体积和内表面积。

数据范围：$1 \leqslant T \leqslant 300; \quad |Z_1|, |Z_2| \leqslant 1000.00$

题解

这是一道纯高数题。

首先对于体积的话非常好搞，学过微积分可以直接算出来：

$$\begin{align*} V= \int_{Z_1}^{Z_2}\pi (2+\cos z)^2 \mathrm{d}z&= \pi \int_{Z_1}^{Z_2} (4+4\cos z+\cos^2 z)\ \mathrm{d}z\\ &= \pi \int_{Z_1}^{Z_2} \left (4+4\cos z+\frac{1+\cos 2z}{2} \right )\ \mathrm{d}z\\ &= \pi \left ( \frac{9z}{2}+4\sin z + \frac{\sin 2z}{4} \right )_{\bigg|_{Z_1}^{Z_2}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \end{align*}$$

然后表面积的话首先有底为$\pi Z_1^2$。

然后就剩下要求侧面积了。

一个错误的想法是去求$\int_{Z_1}^{Z_2}2\pi (2+\cos z)\ dz$，因为就算分得很小点也没法代替弧长。

实际上我们求的时候需要再乘上弧长，弧长的微元有$\mathrm{d}s=\sqrt{1+{y'}^2}\ \mathrm{d}x$

于是有我们要求的就是：

$$S=\pi Z_1^2+2\pi \int_{Z_1}^{Z_2}{(2+\cos z)\sqrt{1+\sin^2 z}}\ \mathrm{d}z\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad$$

后边那坨积分算出来的东西是个椭圆积分，没法用初等函数表示。

然后就需要使用复合$Simpson$来搞。

不会的可以参考：http://blog.csdn.net/xc19952007/article/details/48472033

实现

直接用复合$Simpson$精度实际上是不够的。

考虑三角函数的特殊性我们可以知道$2\pi$区间内的积分值是相同的。

这样我们就可以处理一段$[0,2\pi]$的和砍掉一堆周期之后剩下的长度小于$2\pi$的一段就行了。

曾经有过一个错误的想法是砍成一堆$\pi$，由于想着是个偶函数就这么干了。

实际上由于它的积分是个奇函数，而且并非三角函数，所以能确定的周期只有$2\pi$，$\pi$作为周期是错的。

这个利用$z+2\pi$的代入很好证明。

代码

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <stack>
#include <queue>
#include <string>
#include <vector>
#include <set>
#include <map>

#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define sz(x) ((int)(x).size())
#define fi first
#define se second
#define getmid(l,r) ((l) + ((r) - (l)) / 2)
using namespace std;
typedef long long LL;
typedef vector<int> VI;
typedef pair<int,int> PII;
LL powmod(LL a,LL b, LL MOD) {LL res=1;a%=MOD;for(;b;b>>=1){if(b&1)res=res*a%MOD;a=a*a%MOD;}return res;}
// head
const double PI = acos(-1.0);
const double EPS = 1e-5;

double calV(double l, double r) {
    return (PI*9*(r-l)/2+4*PI*(sin(r)-sin(l))+PI*(sin(r) * cos(r) - sin(l) * cos(l)) / 2);
}

inline double getR(double x) {
    return 2+cos(x);
}
inline double sqr(double x) {
    return x*x;
}
inline double f(double x) {
    return 2*PI*sqrt(1+sqr(sin(x)))*getR(x);
}

inline double getAppr(double fl, double fm, double fr, double l, double r) {
    return (fl+4*fm+fr)*(r-l)/6.0;
}

double Simpson(double l, double r, double fl, double fr) {
    double m = (l+r)/2, lm = (l+m)/2, rm = (r+m)/2;
    double fm = f(m), flm = f(lm), frm = f(rm);
    double vlr = getAppr(fl, fm, fr, l, r);
    double vlm = getAppr(fl, flm, fm, l, m);
    double vrm = getAppr(fm, frm, fr, m, r);
    return fabs(vlr-vlm-vrm) < EPS ? vlr : Simpson(l, m, fl, fm)+Simpson(m, r, fm, fr);
}

double calS(double l, double r) {
    double PIS = Simpson(0, 2*PI, f(0), f(2*PI));
    double p = floor((r-l)/PI/2);
    return PIS*p + Simpson(l+p*PI*2, r, f(l+p*PI*2), f(r));
}

int main() {
    int t;
    double l, r;
    scanf("%d", &t);
    while (t--) {
        scanf("%lf%lf",&l, &r);
        printf("%.2f %.2f\n", calV(l, r), PI*sqr(getR(l))+calS(l, r));
    }
    return 0;
}