数值积分
定积分的定义
简单来说,函数
在区间
上的定积分
指的是
在区间
中与
轴围成的区域的面积(其中
轴上方的部分为正值,
轴下方的部分为负值)。
很多情况下,我们需要高效,准确地求出一个积分的近似值。下面介绍的 辛普森法,就是这样一种求数值积分的方法。
辛普森法
这个方法的思想是将被积区间分为若干小段,每段套用二次函数的积分公式进行计算。
二次函数积分公式(辛普森公式)
对于一个二次函数
,有:
推导过程: 对于一个二次函数
; 求积分可得
在这里 D 是一个常数,那么

根据这个辛普森公式,我们先介绍一种普通的辛普森积分法。
普通辛普森法
1743 年,这种方法发表于托马斯·辛普森的一篇论文中。
描述
给定一个自然数
,将区间
分成
个等长的区间
。

我们就可以计算每个小区间
,
的积分值,将所有区间的积分值相加即为总积分。
对于
,
的一个区间,选其中的三个点
就可以构成一条抛物线从而得到一个函数
,这个函数存在且唯一。计算原函数在该区间的积分值就变成了计算新的二次函数
在该段区间的积分值。这样我们就可以利用辛普森公式来近似计算它。

将其分段求和即可得到如下结论:

误差
我们直接给出结论,普通辛普森法的误差为:

其中
是位于区间
的某个值。
实现
自适应辛普森法
普通的方法为保证精度在时间方面无疑会受到
的限制,我们应该找一种更加合适的方法。
现在唯一的问题就是如何进行分段。如果段数少了计算误差就大,段数多了时间效率又会低。我们需要找到一个准确度和效率的平衡点。
我们这样考虑:假如有一段图像已经很接近二次函数的话,直接带入公式求积分,得到的值精度就很高了,不需要再继续分割这一段了。
于是我们有了这样一种分割方法:每次判断当前段和二次函数的相似程度,如果足够相似的话就直接代入公式计算,否则将当前段分割成左右两段递归求解。
现在就剩下一个问题了:如果判断每一段和二次函数是否相似?
我们把当前段直接代入公式求积分,再将当前段从中点分割成两段,把这两段再直接代入公式求积分。如果当前段的积分和分割成两段后的积分之和相差很小的话,就可以认为当前段和二次函数很相似了,不用再递归分割了。
上面就是自适应辛普森法的思想。在分治判断的时候,除了判断精度是否正确,一般还要强制执行最少的迭代次数。
参考代码如下:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 | double simpson(double l, double r) {
double mid = (l + r) / 2;
return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6; // 辛普森公式
}
double asr(double l, double r, double eqs, double ans, int step) {
double mid = (l + r) / 2;
double fl = simpson(l, mid), fr = simpson(mid, r);
if (abs(fl + fr - ans) <= 15 * eqs && step < 0)
return fl + fr + (fl + fr - ans) / 15; // 足够相似的话就直接返回
return asr(l, mid, eqs / 2, fl, step - 1) +
asr(mid, r, eqs / 2, fr, step - 1); // 否则分割成两段递归求解
}
double calc(double l, double r, double eps) {
return asr(l, r, eps, simpson(l, r), 12);
}
|
1
2
3
4
5
6
7
8
9
10
11
12 | def simpson(l, r):
mid = (l + r) / 2
return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6 # 辛普森公式
def asr(l, r, eqs, ans, step):
mid = (l + r) / 2
fl = simpson(l, mid); fr = simpson(mid, r)
if abs(fl + fr - ans) <= 15 * eqs and step < 0:
return fl + fr + (fl + fr - ans) / 15 # 足够相似的话就直接返回
return asr(l, mid, eqs / 2, fl, step - 1) + \
asr(mid, r, eqs / 2, fr, step - 1) # 否则分割成两段递归求解
def calc(l, r, eps):
return asr(l, r, eps, simpson(l, r), 12)
|
习题
本页面最近更新:2023/1/5 20:18:43,更新历史
发现错误?想一起完善? 在 GitHub 上编辑此页!
本页面贡献者:H-J-Granger, Chrogeek, countercurrent-time, Enter-tainer, Great-designer, iamtwz, Ir1d, ksyx, mao1t, Menci, NachtgeistW, ShaoChenHeng, StudyingFather, SukkaW, zyj-111
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用