自适应辛普森公式
辛普森公式是数值方法中常用的计算函数定积分的近似方法。
- 计算定积分的方法
- 辛普森公式的推导
- 其他N-C公式
- 自适应方法
计算定积分的方法
- 求原函数
- 直接查表
- 这个求不出来怎么办
∫x0t3et−1dt ∫ 0 x t 3 e t − 1 d t
- 那就不求原函数
- 数值积分
辛普森公式的推导
辛普森(Simpson)公式是牛顿-科特斯公式当n=2时的情形,也称为三点公式。利用区间二等分的三个点来进行积分插值。其科特斯系数分别为1/6,4/6,1/6。 —— [ 百度百科 ]
于是我们引入了辛普森公式
∀f(x)=a1x3+a2x2+a3x+a4
∀
f
(
x
)
=
a
1
x
3
+
a
2
x
2
+
a
3
x
+
a
4
∫baf(x)dx=Δx[(f(a)+4f(a+b2)+f(b)]3
∫
a
b
f
(
x
)
d
x
=
Δ
x
[
(
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
]
3
Δx=a+b2
Δ
x
=
a
+
b
2
怎么推导呢?
配方法!!
∫baf(x)dx=Δx[(f(a)+4f(a+b2)+f(b)]3
∫
a
b
f
(
x
)
d
x
=
Δ
x
[
(
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
]
3
=14a1(a4−b4)+13a2(a3−b3)+12a2(a2−b2)+a3(a−b)
=
1
4
a
1
(
a
4
−
b
4
)
+
1
3
a
2
(
a
3
−
b
3
)
+
1
2
a
2
(
a
2
−
b
2
)
+
a
3
(
a
−
b
)
=b−a2{13(a1a3+a2a2+a3a+a4)+43[a1(a+b2)3+a2(a+b2)2+a3a+b2+a4]+13(a1b3+a2b2+a3b+a4)}
=
b
−
a
2
{
1
3
(
a
1
a
3
+
a
2
a
2
+
a
3
a
+
a
4
)
+
4
3
[
a
1
(
a
+
b
2
)
3
+
a
2
(
a
+
b
2
)
2
+
a
3
a
+
b
2
+
a
4
]
+
1
3
(
a
1
b
3
+
a
2
b
2
+
a
3
b
+
a
4
)
}
=a+b2(f(a)+4f(a+b2)+f(b)3
=
a
+
b
2
(
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
3
这就说明了:
任何一个不大于三次的函数的定积分都可以写成这样的辛普森公式!
有什么用呢?
用来近似!
假装其他函数的某一部分是一个不大于三次的函数
那么它的定积分也就好求了
用辛普森公式就行了
先用Simpson公式近似,万一近似导致误差太大,那么就二分,分成两段区间分别近似。
怎么近似呢?
直接带公式啊
什么叫误差太大呢?
先算区间
A=Simpson(a,b)
A
=
S
i
m
p
s
o
n
(
a
,
b
)
的近似值,然后分别算一下
L=Simpson(a,a+b2)
L
=
S
i
m
p
s
o
n
(
a
,
a
+
b
2
)
和
R=Simpson(a+b2,b)
R
=
S
i
m
p
s
o
n
(
a
+
b
2
,
b
)
的近似值
假如
L+R−A≤15∗σ(σ是预设的精度)
L
+
R
−
A
≤
15
∗
σ
(
σ
是
预
设
的
精
度
)
那么我们就需要分别重新计算
(a,a+b2)(a+b2,b)
(
a
,
a
+
b
2
)
(
a
+
b
2
,
b
)
的积分值
此时精度的限制要缩小一半
然后我用Java实现了一下
发现貌似记得太清楚了
跟lrj老师写的一模一样
class Simpson {
// 三点simpson法
public double simpson(double a, double b) {
double c = a + (b-a)/2;
return (F(a)+4*F(c)+F(b))*(b-a)/6;
}
// 自适应Simpson公式(递归过程)。已知整个区间[a,b]上的三点simpson值A
public double asr(double a, double b, double eps, double A) {
double c = a + (b-a)/2;
double L = simpson(a, c), R = simpson(c, b);
if(Math.abs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15.0;
return asr(a, c, eps/2, L) + asr(c, b, eps/2, R);
}
// 自适应Simpson公式(主过程)
public double asr(double a, double b, double eps) {
return asr(a, b, eps, simpson(a, b));
}
}