昨天
A
C
M
A
C
M
模拟的时候遇到了一道
S
i
m
p
s
o
n
S
i
m
p
s
o
n
积分相关的题
完全不知道怎么求，我们组
F
i
s
h
m
e
n
F
i
s
h
m
e
n
被
B
y
m
B
y
m
嘲讽了很久
于是今天下午结合各种资料还是看了一下
这个东西不要觉得它看上去讲什么积分很高级 实际上认真推导也不是很难
首先
S
i
m
p
s
o
n
S
i
m
p
s
o
n
积分的式子是这个东西
积分号的下标代表区间左端点，上标代表区间右端点
这个值的意思就是函数在
[
a
,
b
]
[
a
,
b
]
区间与
x
x
轴，
，
x
=
b
x
=
b
围成的面积
首先引入一些东西，
S
i
m
p
s
o
n
S
i
m
p
s
o
n
积分就是用一个二次曲线来无限逼近原曲线
来达到求得近似值的效果
首先引入一个东西 对于一个二次函数
f
(
x
)
=
a
x
2
+
b
x
+
c
f
(
x
)
=
a
x
2
+
b
x
+
c
求积分可得
F
(
x
)
=
∫
x
0
f
(
x
)
d
x
=
a
3
x
3
+
b
2
x
2
+
c
x
+
D
F
(
x
)
=
∫
0
x
f
(
x
)
d
x
=
a
3
x
3
+
b
2
x
2
+
c
x
+
D
在这里
D
D
是一个常数
这个东西具体怎么证等我问了物理组的同学再来填坑
那么
∫
R
L
f
(
x
)
d
x
=
a
3
(
R
3
−
L
3
)
+
b
2
(
R
2
−
L
2
)
+
c
(
R
−
L
)
∫
L
R
f
(
x
)
d
x
=
a
3
(
R
3
−
L
3
)
+
b
2
(
R
2
−
L
2
)
+
c
(
R
−
L
)
∫
R
L
f
(
x
)
d
x
=
(
R
−
L
)
[
a
3
(
L
2
+
R
2
+
L
R
)
+
b
2
(
L
+
R
)
+
c
]
∫
L
R
f
(
x
)
d
x
=
(
R
−
L
)
[
a
3
(
L
2
+
R
2
+
L
R
)
+
b
2
(
L
+
R
)
+
c
]
∫
R
L
f
(
x
)
d
x
=
b
−
a
6
(
2
a
L
2
+
2
a
R
2
+
2
a
L
R
+
3
b
L
+
3
b
R
+
6
c
)
∫
L
R
f
(
x
)
d
x
=
b
−
a
6
(
2
a
L
2
+
2
a
R
2
+
2
a
L
R
+
3
b
L
+
3
b
R
+
6
c
)
∫
R
L
f
(
x
)
d
x
=
b
−
a
6
{
(
a
L
2
+
b
L
+
c
)
+
(
a
R
2
+
b
R
+
c
)
+
4
[
a
(
L
+
R
2
)
2
+
b
(
L
+
R
2
)
+
c
]
}
∫
L
R
f
(
x
)
d
x
=
b
−
a
6
{
(
a
L
2
+
b
L
+
c
)
+
(
a
R
2
+
b
R
+
c
)
+
4
[
a
(
L
+
R
2
)
2
+
b
(
L
+
R
2
)
+
c
]
}
∫
R
L
f
(
x
)
d
x
=
b
−
a
6
[
f
(
L
)
+
f
(
R
)
+
4
f
(
L
+
R
2
)
]
∫
L
R
f
(
x
)
d
x
=
b
−
a
6
[
f
(
L
)
+
f
(
R
)
+
4
f
(
L
+
R
2
)
]
就这样证完了… 然后我们来考虑这个代码怎么实现
double Simpson(double a, double b) {
int mid = (a + b) / 2;
return (b - a) / 6 * (f(a) + f(b) + 4 * f(mid));
}
好了写完了
这样子直接计算我们发现误差非常大，毕竟原图像可能不能很好的用二次函数逼近
自适应Simpson积分
那怎么办呢，我们可以考虑模拟微分的过程，把区间分成无穷小份再把得到的答案积分就是最后的答案了
但是一般我们只要求一个近似值，所以我们就递归到满足正常误差即可
还有一个问题 每次递归下去都有精度误差，那么累加起来可能就和正确答案相差甚远了
这个时候我们可以在递归的时候提高精度要求，尽量让小的答案精确这样就可以了
在刘汝佳的代码中判断是这样写的
double _asr(double l, double r, double v, double eps) {
double mid = (l + r) / 2;
double v1 = simpson(l, mid), v2 = simpson(mid, r);
if(fabs(v1+v2-v) <= 15*eps) return v1+v2 + (v1+v2-v)/15;
else return _asr(l, mid, v1, eps/2) + _asr(mid, r, v2, eps/2);
}大概意思就是把多出的精度损失也加回去，这样子能够比较好的保证精确度
这个就不用讲了 套板子即可
Codes
#include<bits/stdc++.h>
using namespace std;
typedef double db;
const db eps = 1e-6;
db a, b, c, d;
db f(db x) {
return (c * x + d) / (a * x + b);
}
db integral(db x1, db x2) {
db mid = (x1 + x2) / 2;
return (x2 - x1) / 6 * (f(x1) + f(x2) + f(mid) * 4);
}
db Simpson(db x1, db x2, db eps) {
db mid = (x1 + x2) / 2, res1, res2, res;
res = integral(x1, x2); res1 = integral(x1, mid); res2 = integral(mid, x2);
if(fabs(res1 + res2 - res) < eps * 15) return res1 + res2 + (res1 + res2 - res) / 15;
return Simpson(x1, mid, eps / 2) + Simpson(mid, x2, eps / 2);
}
int main() {
#ifndef ONLINE_JUDGE
freopen("4525.in", "r", stdin);
freopen("4525.out", "w", stdout);
#endif
db L, R;
cin >> a >> b >> c >> d >> L >> R;
printf("%.6lf", Simpson(L, R, eps));
return 0;
}