题1:简单积分

Description

设函数 f(x)=ax2+bx+cf(x)=ax^2+bx+c 求如下积分:

abf(x)dx\int_{a}^{b}f(x)dx

Solution

普通高中课程标准实验教科书数学人教A版选修2-2告诉我们,用牛顿-莱布尼茨公式:

接下来的工作交给人类智慧就能完成辣。

题2:复杂积分

Description

设函数 f(x)=(log7x+1)3f(x)=(log_7x+1)^3 求如下积分:

abf(x)dx\int_{a}^{b}f(x)dx

Solution

做完书上的例题后,不就是复合函数求导吗......

其实...你早该发现这是错的,因为求的是 f(x)dxf(x)dx,而不是 f(x)f^{'}(x) 。 2333.... 那么怎么办呢额。


辛普森公式

介绍一下辛普森(simpson)公式

首先我们知道该积分可以表示为f(x)f(x)[a,b][a,b] 上于 xx轴交的面积(有正负)。

对于函数上的三个点进行积分插值: f(x)=ax2+bx+cf(x)=ax^2+bx+c , 那么

我们简单地推导一下:

对对对√,就是这样。那么对于一个任意的连续函数我们可以将它分为好多个点,然后相邻每三个点类似的求一次面积,然后累加得到总和。显然,分的点的个数越多,精确度越高。

但是精确度越高,所花的时间也就越多了,那我们有没有一种办法可以做到两者间的权衡呢? 有的,介绍一下:

自适应辛普森法

自适应辛普森法

所谓自适应,就是将精确度可能相差较大的区间多分几个点,将精确度可能相差较小的少分几个点。看起来这很ok。那么 可能相差如何定义呢?(定义 c=(a+b)2c=\frac{(a+b)}{2}

我们定义 可能相差较小当且仅当:(以下用S(a,b)S(a,b)表示区间[a,b][a,b]的面积)

S(a,c)+S(c,b)S(a,b)<eps|S(a,c)+S(c,b)-S(a,b)|<eps

此时直接返回结果,否则继续递归调用

abf(x)dx=acf(x)dx+cbf(x)dx\int_{a}^{b}f(x)dx=\int_{a}^{c}f(x)dx+\int_{c}^{b}f(x)dx

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

题目描述

计算积分

LRcx+dax+bdx\int_{L}^{R}\frac{cx+d}{ax+b}dx

结果保留至小数点后66位。

数据保证计算过程中分母不为00且积分能够收敛。

输入输出格式

输入格式: 一行,包含66个实数a,b,c,d,L,Ra,b,c,d,L,R

输出格式: 一行,积分值,保留至小数点后66位。

输入输出样例

输入样例#1:

1 2 3 4 5 6

输出样例#1:

2.732937

说明

      且

Solution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <cstdio>
#include <cmath>
const double eps = 1e-10;
double a,b,c,d,L,R;

inline double F(double x) {
return (c*x+d) / (a*x+b);
}
inline double Simpson(double l, double r) {
return (r-l)/6*(F(l) + 4*F((l+r)/2) + F(r));
}
inline double SOL(double l, double r) {
double mid = (l+r)/2;
if (fabs(Simpson(l, mid) + Simpson(mid, r) - Simpson(l, r)) < eps)
return Simpson(l, r);
return SOL(l, mid) + SOL(mid, r);
}

int main() {
scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &L, &R);
printf("%.6f\n", SOL(L, R));
}

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

题目描述 计算积分

0xaxxdx\int_{0}^{\infty }x^{\frac{a}{x}-x} dx

保留至小数点后55位。若积分发散,请输出"orz"。

输入输出格式

输入格式: 一行,包含一个实数,为aa的值

输出格式: 一行,积分值或 orz

输入输出样例 输入样例#1:

2.33

输出样例#1:

1.51068

说明

请注意时空限制。

Solution

Simpson1

Download
f(x)=xaxx(a=1,5,10,20,50)f(x) = x^{\frac{a}{x}-x}(a=1,5,10,20,50)

Simpson2

Download
f(x)=xaxx(a=1,5,10,20,50)f(x) = x^{\frac{a}{x}-x}(a=-1,-5,-10,-20,-50)

可以注意到当 a<0a < 0 的时候 , limx0f(x)+\lim_{x \to 0} f(x) \to +\infty 积分是不收敛的。当 a>0a > 0 时,积分的收敛边界接近于 88 ,为求稳把边界调大一点。

自适应辛普森法的一个重点是避免特殊点a=0a = 0 时, f(0)=1f(0) = 1a>0a > 0 时, f(0)=0f(0) = 0 所以 x=0x = 0 这个点应该被判掉

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <cstdio>
#include <cmath>
const double eps = 1e-10;
double a;
inline double F(double x) {
return pow(x, a/x-x);
}
inline double Simpson(double l, double r) {
return (r-l)/6*(F(l) + 4*F((l+r)/2) + F(r));
}

inline double SOL(double l, double r, double ans) {
double mid = (l+r)/2, lans=Simpson(l, mid), rans=Simpson(mid, r);
if (fabs(lans + rans - ans) < eps) return ans;
return SOL(l, mid, lans) + SOL(mid, r, rans);
}

int main() {
scanf("%lf", &a);
if(a<0) puts("orz");
else printf("%.5f\n", SOL(eps, 12.0, Simpson(eps, 12.0)));
}

「BZOJ2178」圆的面积并

Time Limit: 20 Sec
Memory Limit: 259 MB

Description

给出N个圆,求其面积并

Input&Output

Input 先给一个数字NN ,     接下来是N行是圆的圆心,半径,其绝对值均为小于 10001000 的整数

Output 面积并,保留三位小数

Sample Input

721
…

Sample Output

12707279.690

Solution

非正解:辛普森积分法 先来画个图:

辛普森积分法1

我们来想想怎么做诶,我们先确定整个图形的 minmin xxmaxmax xx ,然后.....

线 线

诶、做完了?

不!看下图:

辛普森积分法2

那中间的空格怎么办?该怎么减?好像不能这样做辣

我们介绍一种基于 辛普森积分法 求面积并的 新方法

  • 优点:代码复杂度小,容易想到,适用范围广
  • 缺点:误差较大,容易被数据卡掉

在比赛中 骗分 还是挺有用的。

方法介绍:

我们定义函数 f(x1)f(x_1) 为直线 x=x1x=x_1 与图形所交线段的长度,下图蓝线所示:

辛普森积分法3

然后套用辛普森公式,当误差在一定范围内的时候接受得到的近似面积

细节注意:

  1. 需要去除被其它圆所包含的圆,因为没有意义
  2. 计算图形所交线段的长度时要注意两种情况

AC代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
#define re register int
#define rep(i,a,b) for(re i=(a);i<=(b);++i)
#define inf 0x3f3f3f3f
const double eps=1e-12;

int n, skd=0, tag[1020];
int ld=inf, rd=-inf;

struct Circle {
int x, y, r;
friend bool operator < (const Circle& A, const Circle& B) {
return A.r < B.r;
}
} t[1020], c[1020];
struct Segment {
double l, r;
friend bool operator < (const Segment& A, const Segment& B) {
return A.l<B.l || (A.l==B.l&&A.r<B.r);
}
} z[1020];

inline double F(double _x) {
int cnt=0;
rep(i,1,n)
if(_x>=c[i].x-c[i].r && _x<=c[i].x+c[i].r) {
double base = sqrt(c[i].r*c[i].r - (_x-c[i].x)*(_x-c[i].x));
z[++cnt] = (Segment) {c[i].y-base, c[i].y+base};
}
sort(z+1, z+cnt+1);
double ans = 0, h = -inf;
rep(i,1,cnt) {
if(h < z[i].l) ans += z[i].r-z[i].l, h = z[i].r;
else if(z[i].r > h) ans += z[i].r-h, h = z[i].r;
}
return ans;
}

inline double Simpson(double l, double r) {
return (r-l)/6.0 * (F(l) + 4.0*F((l+r)/2) + F(r));
}

inline double SOL(double l, double r, double ans) {
double mid = (l+r)/2, lans = Simpson(l, mid), rans = Simpson(mid, r);
if(fabs(lans + rans - ans) <= eps) return ans;
return SOL(l, mid, lans) + SOL(mid, r, rans);
}

int main() {
scanf("%d", &n);
rep(i,1,n) {
scanf("%d%d%d", &t[i].x, &t[i].y, &t[i].r);
ld = min(ld, t[i].x-t[i].r);
rd = max(rd, t[i].x+t[i].r);
}
sort(t+1, t+n+1);
rep(i,1,n)
rep(j,i+1,n)
if( (t[i].x-t[j].x)*(t[i].x-t[j].x) + (t[i].y-t[j].y)*(t[i].y-t[j].y)
<= (t[j].r-t[i].r)*(t[j].r-t[i].r) ) {tag[i] = 1; break;}
rep(i,1,n)
if(!tag[i]) c[++skd] = t[i];
n = skd;
printf("%.3lf\n", SOL(1.0*ld, 1.0*rd, Simpson(1.0*ld, 1.0*rd)));
}