《数值计算方法》系列总目录
第一章 误差序列实验
第二章 非线性方程f(x)=0求根的数值方法
第三章 CAD模型旋转和AX=B的数值方法
第四章 插值与多项式逼近的数值计算方法
第五章 曲线拟合的数值方法
第六章 数值微分计算方法
第七章 数值积分计算方法
第八章 数值优化方法
第七章
一、算法原理
1、组合梯形公式和辛普森公式原理
在计算函数积分时，传统的方式是将其不定积分求解，再利用牛顿-莱布尼兹公式进行求解，这是求解连续函数定积分的方法。然而在计算机中不能对连续的函数进行求解，因为根据牛顿-莱布尼兹公式推导过程可知，其要求分割的小区间长度趋于无穷小才能得到所求积分，然而在计算机中是无法表达无穷小且有截断误差。因此，遇到此类问题常用的方法就是对函数进行离散化求解。
数值积分的目的是,通过在有限个采样点上计算
f
(
x
)
f(x)
f(x)的值来逼近
f
(
x
)
f(x)
f(x)在区间
[
a
,
b
]
[a,b]
[a,b]上的定积分。因此设
a
=
x
0
<
x
1
<
…
<
X
M
=
b
a= x_0<x_1<…<X_M=b
a=x0<x1<…<XM=b，称形如
Q
[
f
]
=
∑
k
=
0
M
w
k
f
(
x
k
)
=
w
0
f
(
x
0
)
+
w
1
f
(
x
1
)
+
.
.
.
+
w
M
f
(
x
M
)
Q[f] = \sum\limits_{k = 0}^M {{w_k}f({x_k})} = {w_0}f({x_0}) + {w_1}f({x_1}) + ... + {w_M}f({x_M})
Q[f]=k=0∑Mwkf(xk)=w0f(x0)+w1f(x1)+...+wMf(xM) 且具有性质
∫
a
b
f
(
x
)
d
x
=
Q
[
f
]
+
E
[
f
]
\int_a^b {f(x)dx} = Q[f] + E[f]
∫abf(x)dx=Q[f]+E[f] 的公式为数值积分或面积公式。项E[f]称为积分的截断误差,值 称为面积节点, 称为权。根据应用的需要,节点{x}的选择有很多种方法。梯形公式、辛普森公式和布尔公式都选择的是等距节点。通过泰勒多项式的推导可知其误差定义，即设面积公式的精度为正整数n, n使得对所有次数
i
<
=
n
i<=n
i<=n的多项式
P
i
(
x
)
P_i(x)
Pi(x),都满足
E
[
P
i
]
=
0
E[P_i]=0
E[Pi]=0,而对某些次数为
n
+
1
n+1
n+1的多项式
P
n
+
1
(
x
)
P_{n+1}(x)
Pn+1(x)有
E
[
P
n
+
1
]
≠
0
E[P_{n+1}]≠0
E[Pn+1]=0。因此积分公式中的误差项只含有高于
n
n
n阶的多项式。
面积公式的推导有时是基于多项式插值的。前面已经讨论过,过
M
+
1
M +1
M+1个等距点 存在惟一的次数小于等于
M
M
M的多项式
P
M
(
x
)
P_M(x)
PM(x)。当用该多项式来近似
[
a
,
b
]
[a,b]
[a,b]上的
f
(
x
)
f(x)
f(x)时，
P
M
(
x
)
P_M(x)
PM(x)的积分就近似等于
f
(
x
)
f(x)
f(x)的积分,这一结果的公式称为牛顿-科特斯公式。当使用采样点
x
0
=
a
x_0=a
x0=a和
x
M
=
b
x_M=b
xM=b时,称为闭型牛顿–科特斯公式。下面介绍几种常用的积分公式，即
M
=
1
，
2
，
3
M=1，2，3
M=1，2，3和
4
4
4情况下的公式，其中要求结点是等距的，且步长为
h
h
h。
∫
x
0
x
1
f
(
x
)
d
x
≈
h
2
(
f
0
+
f
1
)
\int_{{x_0}}^{{x_1}} {f(x)dx} \approx \frac{h}{2}({f_0} + {f_1})
∫x0x1f(x)dx≈2h(f0+f1)
∫
x
0
x
2
f
(
x
)
d
x
≈
h
3
(
f
0
+
4
f
1
+
f
2
)
\int_{{x_0}}^{{x_2}} {f(x)dx} \approx \frac{h}{3}({f_0} + 4{f_1} + {f_2})
∫x0x2f(x)dx≈3h(f0+4f1+f2)
∫
x
0
x
3
f
(
x
)
d
x
≈
3
h
8
(
f
0
+
3
f
1
+
3
f
2
+
f
3
)
\int_{{x_0}}^{{x_3}} {f(x)dx} \approx \frac{{3h}}{8}({f_0} + 3{f_1} + 3{f_2} + {f_3})
∫x0x3f(x)dx≈83h(f0+3f1+3f2+f3)
∫
x
0
x
4
f
(
x
)
d
x
≈
2
h
45
(
7
f
0
+
32
f
1
+
12
f
2
+
32
f
3
+
7
f
4
)
\int_{{x_0}}^{{x_4}} {f(x)dx} \approx \frac{{2h}}{{45}}(7{f_0} + 32{f_1} + 12{f_2} + 32{f_3} + 7{f_4})
∫x0x4f(x)dx≈452h(7f0+32f1+12f2+32f3+7f4)
其中每一项的精度分别为
n
=
1
,
3
,
3
,
5
n=1,3,3,5
n=1,3,3,5。且误差项分别为
h
3
12
f
(
2
)
(
c
)
,
h
5
90
f
(
4
)
(
c
)
,
3
h
5
80
f
(
4
)
(
c
)
,
8
h
7
945
f
(
6
)
(
c
)
\frac{{{h^3}}}{{12}}{f^{(2)}}(c),\frac{{{h^5}}}{{90}}{f^{(4)}}(c),\frac{{3{h^5}}}{{80}}{f^{(4)}}(c),\frac{{8{h^7}}}{{945}}{f^{(6)}}(c)
12h3f(2)(c),90h5f(4)(c),803h5f(4)(c),9458h7f(6)(c)，可以看出任何一项没有
n
n
n阶以下导数的误差，且当步长
h
h
h越小时，其误差会越小。下面简要给出辛普森公式的证明。
首先，从基于三点的拉格朗日多项式开始证明，其
M
=
2
M=2
M=2，设过
x
0
，
x
1
，
x
2
x_0，x_1，x_2
x0，x1，x2三点的拉格朗日多项式为：
P
2
(
x
)
=
f
0
(
x
−
x
1
)
(
x
−
x
2
)
(
x
0
−
x
1
)
(
x
0
−
x
2
)
+
f
1
(
x
−
x
0
)
(
x
−
x
2
)
(
x
1
−
x
0
)
(
x
1
−
x
2
)
+
f
2
(
x
−
x
0
)
(
x
−
x
1
)
(
x
2
−
x
0
)
(
x
2
−
x
1
)
{P_2}(x) = {f_0}\frac{{(x - {x_1})(x - {x_2})}}{{({x_0} - {x_1})({x_0} - {x_2})}} + {f_1}\frac{{(x - {x_0})(x - {x_2})}}{{({x_1} - {x_0})({x_1} - {x_2})}} + {f_2}\frac{{(x - {x_0})(x - {x_1})}}{{({x_2} - {x_0})({x_2} - {x_1})}}
P2(x)=f0(x0−x1)(x0−x2)(x−x1)(x−x2)+f1(x1−x0)(x1−x2)(x−x0)(x−x2)+f2(x2−x0)(x2−x1)(x−x0)(x−x1) 因此
P
2
（
x
）
P_2（x）
P2（x）是对函数的近似函数，对其的定积分可近似成对原函数的积分，则有
∫
x
0
x
2
f
(
x
)
d
x
=
f
0
∫
x
0
x
2
(
x
−
x
1
)
(
x
−
x
2
)
(
x
0
−
x
1
)
(
x
0
−
x
2
)
d
x
+
f
1
∫
x
0
x
2
(
x
−
x
0
)
(
x
−
x
2
)
(
x
1
−
x
0
)
(
x
1
−
x
2
)
d
x
+
f
2
∫
x
0
x
2
(
x
−
x
0
)
(
x
−
x
1
)
(
x
2
−
x
0
)
(
x
2
−
x
1
)
d
x
\int_{{x_0}}^{{x_2}} {f(x)dx} = {f_0}\int_{{x_0}}^{{x_2}} {\frac{{(x - {x_1})(x - {x_2})}}{{({x_0} - {x_1})({x_0} - {x_2})}}dx} + {f_1}\int_{{x_0}}^{{x_2}} {\frac{{(x - {x_0})(x - {x_2})}}{{({x_1} - {x_0})({x_1} - {x_2})}}dx} + {\rm{ }}{f_2}\int_{{x_0}}^{{x_2}} {\frac{{(x - {x_0})(x - {x_1})}}{{({x_2} - {x_0})({x_2} - {x_1})}}} dx
∫x0x2f(x)dx=f0∫x0x2(x0−x1)(x0−x2)(x−x1)(x−x2)dx+f1∫x0x2(x1−x0)(x1−x2)(x−x0)(x−x2)dx+f2∫x0x2(x2−x0)(x2−x1)(x−x0)(x−x1)dx 现在进行等量代换， ，因此积分限改变，最终可化简为
∫
x
0
x
2
f
(
x
)
d
x
≈
f
0
∫
0
2
h
(
t
−
1
)
h
(
t
−
2
)
(
−
h
)
(
−
2
h
)
h
d
t
+
f
1
∫
0
2
h
(
t
−
0
)
h
(
t
−
2
)
(
h
)
(
−
h
)
h
d
t
+
f
2
∫
0
2
h
(
t
−
0
)
h
(
t
−
1
)
(
2
h
)
(
h
)
h
d
t
=
h
3
(
f
0
+
4
f
1
+
f
2
)
\int_{{x_0}}^{{x_2}} {f(x)dx} \approx {f_0}\int_0^2 {\frac{{h(t - 1)h(t - 2)}}{{( - h)( - 2h)}}hdt} + {f_1}\int_0^2 {\frac{{h(t - 0)h(t - 2)}}{{(h)( - h)}}hdt} + {\rm{ }}{f_2}\int_0^2 {\frac{{h(t - 0)h(t - 1)}}{{(2h)(h)}}hdt} \\ {\rm{ }} = \frac{h}{3}({f_0} + 4{f_1} + {f_2})
∫x0x2f(x)dx≈f0∫02(−h)(−2h)h(t−1)h(t−2)hdt+f1∫02(h)(−h)h(t−0)h(t−2)hdt+f2∫02(2h)(h)h(t−0)h(t−1)hdt=3h(f0+4f1+f2) 最终得到辛普森公式。其误差项可由拉格朗日多项式的插值余项得到。可以看见，使用辛普森公式必须要求其结点等距。而且，当结点数越多，同一个区间上步长就会越小，其拉格朗日多项式逼近程度越好，其积分近似值就越精确。有了对四个牛顿-科特斯公式的推导，就可以推导组合梯形公式和辛普森公式。从梯形公式可以看出，对于一个较大区间进行积分求解，因其步长h较大，会导致特别大的误差，因此可以考虑将大区间分割成一段一段的小区间，再对每一个小区间进行梯形公式，最终得到的公式即为组合梯形公式。因此，组合辛普森公式也是一样的道理。现在给出组合梯形公式和组合辛普森公式的表达式。
设在区间
[
a
,
b
]
[a,b]
[a,b]上有等距节点
[
x
k
,
x
k
+
1
]
[{x_k},{x_{k + 1}}]
[xk,xk+1]。其中
k
=
0
,
1
,
…
M
−
1
k=0,1,…M-1
k=0,1,…M−1，因此共有M个子区间，其步长为
h
=
(
b
−
a
)
/
M
h=(b-a)/M
h=(b−a)/M，其中
M
M
M个子区间的组合梯形公式有3中表达方式，
T
(
f
,
h
)
=
h
2
∑
k
=
1
M
(
f
(
x
k
−
1
)
+
f
(
x
k
)
)
T
(
f
,
h
)
=
h
2
(
f
0
+
2
f
1
+
2
f
2
+
.
.
.
+
2
f
M
−
2
+
2
f
M
−
1
+
f
M
)
T
(
f
,
h
)
=
h
2
(
f
(
a
)
+
f
(
b
)
)
+
h
∑
k
=
1
M
−
1
f
(
x
k
)
T(f,h) = \frac{h}{2}\sum\limits_{k = 1}^M {\left( {f({x_{k - 1}}) + f({x_k})} \right)} \\ T(f,h) = \frac{h}{2}\left( {{f_0} + 2{f_1} + 2{f_2} + ... + 2{f_{M - 2}} + 2{f_{M - 1}} + {f_M}} \right)\\ T(f,h) = \frac{h}{2}\left( {f(a) + f(b)} \right) + h\sum\limits_{k = 1}^{M - 1} {f({x_k})}
T(f,h)=2hk=1∑M(f(xk−1)+f(xk))T(f,h)=2h(f0+2f1+2f2+...+2fM−2+2fM−1+fM)T(f,h)=2h(f(a)+f(b))+hk=1∑M−1f(xk) 其中
f
f
f代表所积分的函数，
h
h
h表示区间的步长。
同样，
M
M
M个子区间的组合辛普森公式有3中表达方式，
S
(
f
,
h
)
=
h
3
∑
k
=
1
M
(
f
(
x
2
k
−
2
)
+
4
f
(
x
2
k
−
1
)
+
f
(
x
2
k
)
)
S
(
f
,
h
)
=
h
3
(
f
0
+
4
f
1
+
2
f
2
+
4
f
3
+
.
.
.
+
2
f
2
M
−
2
+
4
f
2
M
−
1
+
f
2
M
)
S
(
f
,
h
)
=
h
3
(
f
(
a
)
+
f
(
b
)
)
+
2
h
3
∑
k
=
1
M
−
1
f
(
x
2
k
)
+
4
h
3
∑
k
=
1
M
−
1
f
(
x
2
k
−
1
)
S(f,h) = \frac{h}{3}\sum\limits_{k = 1}^M {\left( {f({x_{2k - 2}}) + 4f({x_{2k - 1}}) + f({x_{2k}})} \right)} \\ S(f,h) = \frac{h}{3}\left( {{f_0} + 4{f_1} + 2{f_2} + 4{f_3} + ... + 2{f_{2M - 2}} + 4{f_{2M - 1}} + {f_{2M}}} \right)\\ S(f,h) = \frac{h}{3}\left( {f(a) + f(b)} \right) + \frac{{2h}}{3}\sum\limits_{k = 1}^{M - 1} {f({x_{2k}})} + \frac{{4h}}{3}\sum\limits_{k = 1}^{M - 1} {f({x_{2k - 1}})}
S(f,h)=3hk=1∑M(f(x2k−2)+4f(x2k−1)+f(x2k))S(f,h)=3h(f0+4f1+2f2+4f3+...+2f2M−2+4f2M−1+f2M)S(f,h)=3h(f(a)+f(b))+32hk=1∑M−1f(x2k)+34hk=1∑M−1f(x2k−1) 其中，对
M
M
M个子区间依次进行梯形公式和辛普森公式即可得到最终表达式。利用拉格朗日多项式的余项可以得出，组合梯形公式和组合辛普森公式的误差分别为
E
T
(
f
,
h
)
=
−
(
b
−
a
)
f
(
2
)
(
c
)
h
2
12
E
S
(
f
,
h
)
=
−
(
b
−
a
)
f
(
4
)
(
c
)
h
4
180
{E_T}(f,h) = \frac{{ - (b - a){f^{(2)}}(c){h^2}}}{{12}}\\ {E_S}(f,h) = \frac{{ - (b - a){f^{(4)}}(c){h^4}}}{{180}}
ET(f,h)=12−(b−a)f(2)(c)h2ES(f,h)=180−(b−a)f(4)(c)h4 因此可以看出，在同一区间上，利用辛普森公式可以得到在较大步长时也能得到理想精度，因此其计算量会会大大减少，且精确度也会提高，这点在后面的实例分析中会展示。
2、递归公式和龙贝格积分原理
由第一节可知，使用的子区间数越多,则逼近的精度越高。如何选择子区间的数目?采用2个子区间,4个子区间,…,不断进行试验,直至得到想要的精度,这个过程能够帮助回答这一问题。首先要生成一个梯形公式的序列
T
(
J
)
{T(J)}
T(J),，当子区间数目增加一倍时,函数求值的次数也增加了近一倍,因为必须在所有先前的点和先前区间的中点处对函数进行求值。
设
J
>
=
1
J>=1
J>=1，点
x
k
=
a
+
k
h
{x_k=a+kh}
xk=a+kh将
[
a
,
b
]
[a,b]
[a,b]划分成
2
J
=
2
M
2^J=2M
2J=2M个宽度为(b-a)/2M的子区间，那么梯形公式
T
(
f
,
h
)
T(f,h)
T(f,h)和
T
(
f
,
2
h
)
T(f,2h)
T(f,2h)满足
T
(
f
,
h
)
=
T
(
f
,
2
h
)
2
+
h
∑
k
=
1
M
f
(
x
2
k
−
1
)
T(f,h) = \frac{{T(f,2h)}}{2} + h\sum\limits_{k = 1}^M {f({x_{2k - 1}})}
T(f,h)=2T(f,2h)+hk=1∑Mf(x2k−1) 其中
T
(
f
,
h
)
T(f,h)
T(f,h)是对函数使用步长为h的梯形公式求得的积分近似值，并将其表示为
T
(
J
)
T(J)
T(J)。称式为递归梯形公式。其中
h
=
(
b
−
a
)
/
2
J
h=(b-a)/2^J
h=(b−a)/2J,
x
k
=
a
+
k
h
{xk=a+kh}
xk=a+kh。其证明过程简要为：对偶数结点进行步长为
2
h
2h
2h和
h
h
h的梯形公式，然后利用等量代换法可以得到关于
T
(
J
)
T(J)
T(J)和
T
(
J
−
1
)
T(J-1)
T(J−1)的关系。同样，可以对
T
(
J
)
T(J)
T(J)和
T
(
J
−
1
)
T(J-1)
T(J−1)进行另外的线性组合，从而得到推导出辛普森公式
S
(
f
,
h
)
S(f,h)
S(f,h)的公式，即
S
(
J
)
=
4
T
(
J
)
−
T
(
J
−
1
)
3
,
f
o
r
J
=
1
,
2
,
.
.
.
S(J) = \frac{{4T(J) - T(J - 1)}}{3},{\rm{ }}for{\rm{ }}J = 1,2,...
S(J)=34T(J)−T(J−1),forJ=1,2,... 称式为递归辛普森公式。同理，还可以利用推导出的辛普森公式推导出递归布尔公式，同样是对不同步长的辛普森公式进行线性组合，最终得到递归布尔公式为
B
(
J
)
=
16
S
(
J
)
−
S
(
J
−
1
)
15
,
f
o
r
J
=
2
,
3
,
.
.
.
B(J) = \frac{{16S(J) - S(J - 1)}}{{15}},{\rm{ }}for{\rm{ }}J = 2,3,...
B(J)=1516S(J)−S(J−1),forJ=2,3,... 其中递归布尔公式的推导需要有组合布尔公式的表达式，即将区间分为
4
M
4M
4M时的表达式，从而能够根据递归辛普森公式的线性组合写出递归布尔公式。
不难发现，通过一些公式总能够把低阶的积分公式推导出高阶公式，从而获得更高的精度，因此可以对其推广，从而得到由低阶推导高阶的积分公式的通项公式，即龙贝格积分。设用步长
h
h
h和
2
h
2h
2h得到个逼近公式的两个结果,则两个结果的代数运算将得到改进的答案。每次改进将误差项的阶由
O
(
h
2
N
)
O(h^{2N})
O(h2N)提高到
O
(
h
2
N
+
2
)
O(h^{2N+2})
O(h2N+2)。该过程称为龙贝格积分,它有自己的优点和缺点。
在布尔公式之后很少用到牛顿-科特斯公式。这是因为9点牛顿-科特斯面积公式中有负的权值,而所有超过10个点的公式中都有负的权,这会导致由舍入带来的误差。龙贝格积分的优点在于它的所有的权都是正的,而且它使用等距的节点,易于计算横坐标的值。龙贝格积分的一个缺点是,为了将误差由
O
(
h
2
N
)
O(h^{2N})
O(h2N)降低到
O
(
h
2
N
+
2
)
O(h^{2N+2})
O(h2N+2),函数求值次数增加了一倍。使用连续公式能减少计算量。龙贝格积分的推导基于理论上的假设:如果对所有的
N
N
N有
f
∈
C
N
[
a
,
b
]
f∈CN[a,b]
f∈CN[a,b],则梯形公式的误差项可以表示为一个只包含h的偶数次幂的级数,即，
∫
a
b
f
(
x
)
d
x
=
T
(
f
,
h
)
+
E
T
(
f
,
h
)
,
E
T
(
f
,
h
)
=
a
1
h
2
+
a
2
h
4
+
a
3
h
6
+
.
.
.
\int_a^b {f(x)dx} = T(f,h) + {E_T}(f,h),{E_T}(f,h) = {a_1}{h^2} + {a_2}{h^4} + {a_3}{h^6} + ...
∫abf(x)dx=T(f,h)+ET(f,h),ET(f,h)=a1h2+a2h4+a3h6+... 由于公式(13)中只包含h的偶数次项,可以连续地使用理查森改进,首先消去
a
1
a_1
a1,接着消去
a
2
a_2
a2,然后是
a
3
a_3
a3,依次类推。该过程产生偶数阶次的误差项
O
(
h
4
)
,
O
(
h
6
)
,
O
(
h
8
)
O(h^4),O(h^6),O(h^8)
O(h4),O(h6),O(h8)等。于是可以得到关于
T
(
f
,
2
h
)
T(f,2h)
T(f,2h)和
T
(
f
,
h
)
T(f, h)
T(f,h)的关系式。因此可以得到龙贝格积分的理查森改进，即给定
Q
Q
Q的两个逼近
R
(
2
h
,
K
−
1
)
R(2h,K-1)
R(2h,K−1)和
R
(
h
,
K
−
1
)
R(h,K-1)
R(h,K−1)，满足
Q
=
R
(
h
,
K
−
1
)
+
c
1
h
2
K
+
c
2
h
2
K
+
2
+
.
.
.
Q
=
R
(
2
h
,
K
−
1
)
+
c
1
4
K
h
2
K
+
c
2
4
K
+
1
h
2
K
+
2
+
.
.
.
Q = R(h,K - 1) + {c_1}{h^{2K}} + {c_2}{h^{2K + 2}} + ...\\ Q = R(2h,K - 1) + {c_1}{4^K}{h^{2K}} + {c_2}{4^{K + 1}}{h^{2K + 2}} + ...
Q=R(h,K−1)+c1h2K+c2h2K+2+...Q=R(2h,K−1)+c14Kh2K+c24K+1h2K+2+... 通过消元可以得到关系式
Q
=
4
K
R
(
h
,
K
−
1
)
−
R
(
2
h
,
K
−
1
)
4
K
−
1
+
O
(
h
2
K
+
2
)
Q = \frac{{{4^K}R(h,K - 1) - R(2h,K - 1)}}{{{4^K} - 1}} + O({h^{2K + 2}})
Q=4K−14KR(h,K−1)−R(2h,K−1)+O(h2K+2) 其中，误差相的精度提高了2次，通过此方法，可以通过迭代的方法，从而大大降低误差，简化计算复杂度。其过程类似于求解牛顿多项式系数和理查森外推法一样，构造一张记录每个步长以及每一阶积分值的表格。
3、自适应积分原理
组合积分公式要求使用等距节点。典型的情况是,在整个积分区间使用相同的小步长
h
h
h，以保证整体精度。这并没有考虑到曲线的某些部分变化剧烈,需要比其他部分多加考虑的情况,因此引入一种能够在函数值变化大的区间采用较小步长的方法是很有用的。该技术称为自适应积分,它的基础是辛普森公式。
区间
[
a
k
,
b
k
]
[a_k,b_k]
[ak,bk]的4个子区间上的组合辛普森公式可用另一种方法实现:先将该区间划分为两个相等的子区间
[
a
k
1
,
b
k
1
]
[a_{k1},b_{k1}]
[ak1,bk1]和
[
a
k
2
,
b
k
2
]
[a_{k2},b_{k2}]
[ak2,bk2],并在两段上递归地应用辛普森公式。这只需要增加两个
f
(
x
)
f(x)
f(x)求值计算,其结果为
S
(
a
k
1
,
b
k
1
)
+
S
(
a
k
2
,
b
k
2
)
=
h
6
(
f
(
a
k
1
)
+
4
f
(
c
k
1
)
+
f
(
b
k
1
)
)
+
h
6
(
f
(
a
k
2
)
+
4
f
(
c
k
2
)
+
f
(
b
k
2
)
)
S({a_{k1}},{b_{k1}}) + S({a_{k2}},{b_{k2}}) = \frac{h}{6}\left( {f({a_{k1}}) + 4f({c_{k1}}) + f({b_{k1}})} \right) + \frac{h}{6}\left( {f({a_{k2}}) + 4f({c_{k2}}) + f({b_{k2}})} \right)
S(ak1,bk1)+S(ak2,bk2)=6h(f(ak1)+4f(ck1)+f(bk1))+6h(f(ak2)+4f(ck2)+f(bk2)) 其中,
a
k
1
=
a
k
，
b
k
1
=
a
k
2
=
c
k
，
b
k
2
=
b
k
，
c
k
1
a_{k1}=a_k，b_{k1}=a_{k2}=c_k，b_{k2}=b_k，c_{k1}
ak1=ak，bk1=ak2=ck，bk2=bk，ck1是
[
a
k
1
,
b
k
1
]
[a_{k1},b_{k1}]
[ak1,bk1]的中点，
c
k
1
c_{k1}
ck1是
[
a
k
2
,
b
k
2
]
[a_{k2},b_{k2}]
[ak2,bk2]的中点。因此可以得到式与辛普森公式的关系，即
S
(
a
k
,
b
k
)
−
h
5
f
(
4
)
(
d
2
)
90
=
S
(
a
k
1
,
b
k
1
)
+
S
(
a
k
2
,
b
k
2
)
−
h
5
16
f
(
4
)
(
d
2
)
90
S({a_k},{b_k}) - {h^5}\frac{{{f^{(4)}}({d_2})}}{{90}} = S({a_{k1}},{b_{k1}}) + S({a_{k2}},{b_{k2}}) - \frac{{{h^5}}}{{16}}\frac{{{f^{(4)}}({d_2})}}{{90}}
S(ak,bk)−h590f(4)(d2)=S(ak1,bk1)+S(ak2,bk2)−16h590f(4)(d2) 其中，等式左边为辛普森公式的结果，右边为在区间上左右等分区间上的辛普森公式的结果，其误差项也可以计算出来。合并误差项可以得到
−
h
5
f
(
4
)
(
d
2
)
90
≈
16
15
(
S
(
a
k
1
,
b
k
1
)
+
S
(
a
k
2
,
b
k
2
)
−
S
(
a
k
,
b
k
)
)
- {h^5}\frac{{{f^{(4)}}({d_2})}}{{90}} \approx \frac{{16}}{{15}}\left( {S({a_{k1}},{b_{k1}}) + S({a_{k2}},{b_{k2}}) - S({a_k},{b_k})} \right)
−h590f(4)(d2)≈1516(S(ak1,bk1)+S(ak2,bk2)−S(ak,bk)) 因此将式带入式可以得到误差估计：
∣
∫
a
k
b
k
f
(
x
)
d
x
−
S
(
a
k
1
,
b
k
1
)
−
S
(
a
k
2
,
b
k
2
)
∣
≈
1
15
∣
S
(
a
k
1
,
b
k
1
)
+
S
(
a
k
2
,
b
k
2
)
−
S
(
a
k
,
b
k
)
∣
\left| {\int_{{a_k}}^{{b_k}} {f(x)dx} - S({a_{k1}},{b_{k1}}) - S({a_{k2}},{b_{k2}})} \right| \approx \frac{1}{{15}}\left| {S({a_{k1}},{b_{k1}}) + S({a_{k2}},{b_{k2}}) - S({a_k},{b_k})} \right|
∣∣∣∣∣∫akbkf(x)dx−S(ak1,bk1)−S(ak2,bk2)∣∣∣∣∣≈151∣S(ak1,bk1)+S(ak2,bk2)−S(ak,bk)∣ 将系数换成
1
/
10
1/10
1/10，即可得到精度测试。
对于区间
[
a
k
,
b
k
]
[a_k,b_k]
[ak,bk]指定容差
ε
k
>
0
ε_k>0
εk>0。如果
1
10
∣
S
(
a
k
1
,
b
k
1
)
+
S
(
a
k
2
,
b
k
2
)
−
S
(
a
k
,
b
k
)
∣
<
ε
k
\frac{1}{{10}}\left| {S({a_{k1}},{b_{k1}}) + S({a_{k2}},{b_{k2}}) - S({a_k},{b_k})} \right| < {\varepsilon _k}
101∣S(ak1,bk1)+S(ak2,bk2)−S(ak,bk)∣<εk 则推断
∣
∫
a
k
b
k
f
(
x
)
d
x
−
S
(
a
k
1
,
b
k
1
)
−
S
(
a
k
2
,
b
k
2
)
∣
<
ε
k
\left| {\int_{{a_k}}^{{b_k}} {f(x)dx} - S({a_{k1}},{b_{k1}}) - S({a_{k2}},{b_{k2}})} \right| < {\varepsilon _k}
∣∣∣∣∣∫akbkf(x)dx−S(ak1,bk1)−S(ak2,bk2)∣∣∣∣∣<εk 即认为利用公式逼近函数的积分是在精度范围内的。因为推推导可知，式左边是逼近精度很好的近似，因此只要满足了式，就有理由认为近似值满足精度要求。其中对于每个区间都有一个误差限
ε
k
ε_k
εk。
通过应用辛普森公式实现自适应积分。从
{
[
a
0
,
b
0
]
,
ε
0
}
\{[a_0,b_0],ε_0\}
{[a0,b0],ε0}开始,其中
ε
0
ε_0
ε0是区间
[
a
0
,
b
0
]
[a_0,b_0]
[a0,b0]上数值积分的容差。该区间细分为两个子区间,分别记为KaTeX parse error: Expected '}', got 'EOF' at end of input: [a_{01],b_{01}]和
[
a
02
,
b
02
]
[a_{02},b_{02}]
[a02,b02]。如果通过了精度测试,则在区间
[
a
0
,
b
0
]
[a_0,b_0]
[a0,b0]上应用枳分公式,过在测试;否则,如果没有通过测试,则将两个子区间重新标记为
[
a
1
,
b
1
]
[a_1,b_1]
[a1,b1]和
[
a
2
,
b
2
]
[a_2,b_2]
[a2,b2],其上分别采用容差
ε
1
=
ε
0
/
2
ε_1=ε_0/2
ε1=ε0/2和
ε
2
=
ε
0
/
2
ε_2=ε_0/2
ε2=ε0/2。
这样就得到了两个子区间及其相应的容差,对它们进一步细分和测试:
{
[
a
1
,
b
1
]
,
ε
1
}
\{[a_1,b_1],ε_1\}
{[a1,b1],ε1}和
{
[
a
2
,
b
2
]
,
ε
2
}
\{[a_2,b_2],ε_2\}
{[a2,b2],ε2},其中
ε
1
+
ε
2
=
ε
0
ε_1+ε_2=ε_0
ε1+ε2=ε0。如果必须继续进行自适应积分,则必须进一步细分子区间并进行测试,每个子区间都有相应的容差。
在第2步中首先考虑
{
[
a
1
,
b
1
]
,
ε
1
}
\{[a_1,b_1],ε_1\}
{[a1,b1],ε1},并将区间
[
a
1
,
b
1
]
[a_1,b_1]
[a1,b1]细分为
[
a
11
,
b
11
]
[a_{11},b_{11}]
[a11,b11]和
[
a
12
,
b
12
]
[a_{12},b_{12}]
[a12,b12]。如果它们按照容差
ε
1
ε_1
ε1通过了精度测试,则在区间
[
a
1
,
b
1
]
[a_1,b_1]
[a1,b1]上应用公式能保证此区间上的精度;否则,如果不能以容差
ε
1
ε_1
ε1通过精度测试,则第3步必须对两个子区间
[
a
11
,
b
11
]
[a_{11},b_{11}]
[a11,b11]和
[
a
12
,
b
12
]
[a_{12},b_{12}]
[a12,b12]细分,并以减小的容差进行测试。在第2步中还要考察
{
[
a
2
,
b
2
]
,
ε
2
}
\{[a_2,b_2],ε_2\}
{[a2,b2],ε2},并将
[
a
2
,
b
2
]
[a_2,b_2]
[a2,b2]细分为
[
a
21
,
b
21
]
[a_{21},b_{21}]
[a21,b21]和
[
a
22
,
b
22
]
[a_{22},b_{22}]
[a22,b22],如果它们按照容差ε2通过了测试,则在区间
[
a
2
,
b
2
]
[a_2,b_2]
[a2,b2]上应用公式能保证此区间上的精度;否则,如果不能以容差
ε
2
ε_2
ε2通过精度测试,则第3步必须对两个子区间
[
a
21
,
b
21
]
[a_{21},b_{21}]
[a21,b21]和
[
a
22
,
b
22
]
[a_{22},b_{22}]
[a22,b22]细分,并以减小的容差
ε
2
ε_2
ε2进行测试。
如果必须继续进行自适应积分,则必须以各自相应的容差测试较小的区间。公式中的误差项显示,每一对一个较小的区间进行细分,误差的衰减因子大约为
1
/
16
1/16
1/16。在有限步之后停止。该方法的实现过程中需要一个信号变量,用以指示每个子区间是否通过了精度测试。为避免不必要的
f
(
x
)
f(x)
f(x)求值计算,函数值可以在对应于每个子区间的一个数据表中给出。
自适应积分主要是通过分析误差项得出来的判定依据，当不满足精度测试时，我们便为在这一段区间上需要更加细分的区间及其上的辛普森公式才能得到精度范围内的解。基于这个特性，很容易设计出递归函数算法计算积分，通过设置初始误差，即可得到在误差范围内的积分值，因为每个小区间的误差限的和即为初始误差，则每个区间上的误差和也不会超过初始误差，否则就需要继续细分迭代。从而最终可以得到在精度范围内的积分值。
二、实验内容及核心算法代码
1、组合梯形公式和辛普森公式原理实现
由第一节原理可知，在计算机中不能对连续的函数进行求解。遇到此类问题常用的方法就是对函数进行离散化求解。面积公式的推导有时是基于多项式插值的。前面已经讨论过,过
M
+
1
M +1
M+1个等距点 存在惟一的次数小于等于
M
M
M的多项式
P
M
(
x
)
P_M(x)
PM(x)。当用该多项式来近似
[
a
,
b
]
[a,b]
[a,b]上的
f
(
x
)
f(x)
f(x)时，
P
M
(
x
)
P_M(x)
PM(x)的积分就近似等于
f
(
x
)
f(x)
f(x)的积分,这一结果的公式称为牛顿-科特斯公式。当使用采样点
x
0
=
a
x_0=a
x0=a和
x
M
=
b
x_M=b
xM=b时,称为闭型牛顿–科特斯公式。
利用闭型牛顿–科特斯公式求解数值积分的过程较为简单，唯一需要注意的地方就是选取结点的过程，即在考虑计算精度时还要考虑计算复杂度，因此要考虑选取选择区间的个数。
因此，利用组合梯形公式和组合辛普森公式求解数值积分的基本过程为：
- 确定所积分的区间和所分割区间的个数。其中组合梯形公式所分割的区间为 M M M，组合辛普森公式所分割区间个数为 2 M 2M 2M。
- 求解所有结点的横坐标，并求解其上的函数值。
- 利用组合梯形公式和组合辛普森公式，根据已（注意结点的函数值和的系数）。
因此，通过以下实例对实验结果进行测试。
- 修改组合梯形公式,使之可以求只有若干点函数值已知的函数积分。编写程序求区间 [ a , b ] [a,b] [a,b]上过 M M M个给定点的函数 f ( x ) f(x) f(x)的积分逼近。注意节点不需要等距。利用该程序求过点 的函数的积分逼近。
当选取节点不等距时，就不能利用等距节点推导出的组合梯形公式和组合辛普森公式，只能对每相邻节点使用梯形公式或辛普森公式求解分段的积分值，在对所有区间进行求和。但是直接对每个区间进行梯形公式得到的结果精确度较低，因此另外一种方法是先对离散的结点进行拉格朗日插值拟合，再对拟合后的曲线进行区间上组合梯形公式和组合辛普森公式，但是这种方法的计算量相对较大。
其核心函数算法代码为：
//use the combination trapezuim formula to cmp the integration
float CmpIntegrationvalbyComTrap(int N)
{
float h, xf, xb, yf, yb;
float integ = 0.;
xb= Cmp_abscissa(0.);
for (int i = 1; i <= N; ++i)
{
xf = xb;
xb = Cmp_abscissa(float(i));
yf = Cmp_ordinate(xf);
yb = Cmp_ordinate(xb);
h = xb - xf;
integ += h * (yb + yf) / 2.;
}
return integ;
}
//cmp the integration of the function which only several knots has been known,the Fval is the value of each knots
//the numofkont is num of knots,konts_abscissa is the abscissa of discrete knots.
//use the Lagerange interpolation to process the discrete knots,and get the val in the xxs,have N+1 xxs,and the step-length is h of xxs
float CmpUnequaldisIntergrationbyComTrape(float* konts_abscissa, float* Fval, const int numofkont, float* xxs, const int N, float h)
{
//Lagerange interpolation to process firstly
float* ypt;
ypt = MallocArr(N + 1);
for (int i = 0; i <= N; ++i)
{
ypt[i] = LagrangePolyapproximate(konts_abscissa, Fval, xxs[i], numofkont);
}
//use the combination trapezium formula to cmp the integration of xxs,ypt
float integ = h * (ypt[0] + ypt[N]) / 2.;
float sumF = 0.;
for (int i = 1; i < N; ++i)
{
sumF += ypt[i];
}
integ += h * sumF;
return integ;
}
//cmp the integration by combination simpson formula,the arguments info is the same as theCmpUnequaldisIntergrationbyComTrape
float CmpUnequaldisIntergrationbyComSim(float* konts_abscissa, float* Fval, const int numofkont, float* xxs, const int N, float h)
{
//Lagerange interpolation to process firstly
float* ypt;
ypt = MallocArr(N + 1);
for (int i = 0; i <= N; ++i)
{
ypt[i] = LagrangePolyapproximate(konts_abscissa, Fval, xxs[i], numofkont);
}
//use the combination trapezium formula to cmp the integration of xxs,ypt
float integ = h * (ypt[0] + ypt[N]) / 3.;
float sumF1 = 0., sumF2 = 0.;
int M = N / 2;//divede into 2*M segments
for (int i = 1; i < M; ++i)
{
sumF1 += ypt[i * 2];//cmp the sum of f(x|2k)
}
for (int i = 1; i <= M; ++i)
{
sumF2 += ypt[i * 2 - 1];//cmp the sum of f(x|2k-1)
}
integ += (2 * h * sumF1 + 4 * h * sumF2) / 3.;
return integ;
}--
2、递归公式和龙贝格积分原理实现
由第一节原理可知，采用2个子区间,4个子区间,…,不断进行试验,直至得到想要的精度,生成一个梯形公式的序列
{
T
(
J
)
}
\{T(J)\}
{T(J)},，当子区间数目增加一倍时,函数求值的次数也增加了近一倍,因为必须在所有先前的点和先前区间的中点处对函数进行求值。除此之外，通过一些公式总能够把低阶的积分公式推导出高阶公式，从而获得更高的精度，因此可以对其推广，从而得到由低阶推导高阶的积分公式的通项公式，即龙贝格积分。
因此，在利用龙贝格积分求解数值积分时，会生成一个积分数值表，其中第一列即为递归梯形公式不同步长所求得的积分值，第二列为递归辛普森公式所求的值，第三列为递归布尔公式所求得到的积分值，后面即为更高阶公式所求得的积分值。因此，求解龙贝格积分值表即可求得其任意阶任意步长的积分值。使用龙贝格积分求解数值积分的流程为：
- 确定所积函数以及积分区间，利用梯形公式求解龙贝格积分值表第一列的数值。
- 从第二行开始，利用龙贝格积分的理查森改进递推公式求解每一行上的元素，其中第i行就有i个元素。
- 取最后一行的最后一列的元素作为最后的函数积分值。
下面通过以下实例进行实验。
- 正态概率密度函数为 f ( t ) = ( 1 / 2 π ) e − t 2 / 2 f(t) = (1/\sqrt {2\pi } ){e^{ - {t^2}/2}} f(t)=(1/2π)e−t2/2,而累积分布为由积分 f ( t ) = ( 1 / 2 π ) e − t 2 / 2 f(t) = (1/\sqrt {2\pi } ){e^{ - {t^2}/2}} f(t)=(1/2π)e−t2/2定义的函数。计算有8位有效数字的 Φ ( 0.5 ) \Phi (0.5) Φ(0.5), Φ ( 1.0 ) \Phi (1.0) Φ(1.0) , Φ ( 1.5 ) \Phi (1.5) Φ(1.5) Φ ( 2.0 ) \Phi (2.0) Φ(2.0), Φ ( 2.5 ) \Phi (2.5) Φ(2.5), Φ ( 3.0 ) \Phi (3.0) Φ(3.0), Φ ( 3.5 ) \Phi (3.5) Φ(3.5), Φ ( 4.0 ) \Phi (4.0) Φ(4.0)的值。
计算每个
x
x
x的函数积分值，仅仅需要对正态概率密度函数
f
(
t
)
f(t)
f(t)进行
0
0
0到
x
x
x上的积分，可以用龙贝格积分公式求解，然后乘以系数和加上常数项即可得到最终的结果。
其核心函数算法代码为：
//cmp the first col of integtab,and R(J,0) is trapezuim formula,N is the 6length of the first col
void CmpIntegrationbyConTrape(float (*fun)(float), float* T, float xs, float xe, const int N)
{
float h = xe - xs;//the initial step-length
int segment = 1, clock = 0;
T[0] = h * (fun(xs) + fun(xe)) / 2.;
do
{
++clock;
h /= 2.;
segment *= 2;
int M = segment / 2;
float sumf = 0., tx;
for (int i = 1; i <= M; ++i)
{
tx = xs + (2 * i - 1) * h;
sumf += fun(tx);
}
T[clock] = T[clock - 1] / 2. + h * sumf;
} while (clock < N - 1);
}
//use the Romberg to cmp integration,the shape of the integtab is N*N
float CmpIntegrationbyRomberg(float (*fun)(float), float** integtab, float xs, float xe, const int N)
{
//cmp the 1st col,the continuous trapezuim formula
float h = xe - xs;
float* T;
T = MallocArr(N);
CmpIntegrationbyConTrape(fun, T, xs, xe, N);
for (int i = 0; i < N; ++i)
{
integtab[i][0] = T[i];
}
FreeArr(T);
//cmp the rest element of the integtab
float val;
//from the 2nd col
for (int i = 1; i < N; ++i)
{
//for every row
for (int j = 1; j <= i; ++j)
{
val = integtab[i][j - 1] - integtab[i - 1][j - 1];
integtab[i][j] = integtab[i][j - 1] + val / (powf(4, j) - 1);
}
}
return .5 + integtab[N - 1][N - 1];
}
3、自适应积分原理原理实现
由第一节原理可知，组合积分公式要求使用等距节点。典型的情况是,在整个积分区间使用相同的小步长
h
h
h，以保证整体精度。这并没有考虑到曲线的某些部分变化剧烈,需要比其他部分多加考虑的情况,因此引入一种能够在函数值变化大的区间采用较小步长的方法是很有用的。该技术称为自适应积分,它的基础是辛普森公式。
自适应积分的原理推导较为简单，然而计算过程非常繁琐，主要是要一层一层的考虑每个区间的精度测试，然后对没有通过测试的区间进行进一步的细分，最终将所有区间的积分值相加得到最终的结果。
因此，接可以达到在函数变化率较大的地方实现更精细的区间划分，因为函数变化较大的地方，其辛普森公式的逼近效果并不理想，因此一般不会通过精度测试，这就提示我们可以对其进行更加细分的划分。
其核心算法为递归算法，因此其基本流程为：
- 确定所积函数以及积分区间，设定最终结果的误差上限，即初始容差。
- 对区间 [ a k , b k ] [a_k,b_k] [ak,bk]进行组合辛普森公式求值，得到整个区间和两个小区间的积分近似值。
- 对所求得的函数积分近似值进行精度测试。若未通过，则对两个小区间分别进行步骤（2），直到通过精度测试；若通过，则以两个区间所求的组合辛普森公式的积分近似值的和作为函数积分值。
- 直到所有区间都通过了精度测试，将所有区间的积分值相加得到最终的结果。
下面通过以下示例来对自适应积分的实用性进行测试。
- 求以下定积分的近似值,起始容差ε0= 0.00001.
对以上每个定积分,绘制一个包含结点显示和函数曲线的图。
对以上每个函数可知，每个函数都有变化比较大的地方，因此适合用自适应积分进行求值。求解过程千篇一律，仅仅是换了积分函数和积分区间。记录其结点位置，使用matlab进行绘图。
其核心函数算法代码为：
//cmp the integration value by simpson in [a,b],and save S(ak1,bk1) S(ak2,bk2) S(ak,bk) to Svec
void CmpSimpsonval(float (*fun)(float), float a, float b, float& Sk, float& Sk1, float& Sk2)
{
float c = (a + b) / 2.; //the middle of [a,b]
float h = (b - a) / 2.; //the step [a,b]
float c1 = (c + a) / 2., c2 = (b + c) / 2.; //the middle of[a,c],[c,b]
Sk1 = h * (fun(a) + 4 * fun(c1) + fun(c)) / 6.; //S(ak1,bk1)
Sk2 = h * (fun(c) + 4 * fun(c2) + fun(b)) / 6.; //S(ak2,bk2)
Sk = h * (fun(a) + 4 * fun(c) + fun(b)) / 3.; //S(ak,bk)
}
//cmp the adaptive integration and save it to the integtab,and note the iterate time iter and the sum of error
//a0:the initial left abscissa
//b0:the initial right abscissa
//eps:the initial tolerance
//when not meet the tolerance,adapt to the left,so every time note the smaller side
void CmpAdaptiveIntegbySimpson(float (*fun)(float), float a0, float b0, float eps, float& integ, int& iter, float** integtab, float* abscissa_knot, float& sumErr)
{
float Sk, Sk1, Sk2;
CmpSimpsonval(fun, a0, b0, Sk, Sk1, Sk2);//cmp the S(ak1,bk1), S(ak2,bk2), S(ak,bk)
if (abs(Sk1 + Sk2 - Sk) / 10. < eps)
{
float sint = Sk1 + Sk2;
integ += sint;
sumErr = abs(sint - Sk);
integtab[iter][0] = a0; //note the data and save to the integtab,left knot
integtab[iter][1] = b0; //right side
integtab[iter][2] = sint; //the S(ak1,bk1)+S(ak2,bk2)
integtab[iter][3] = sumErr; //sum of err
integtab[iter][4] = eps; //every step's eps
abscissa_knot[iter] = a0; //save the every knot
++iter;
return;
}
else
{
eps /= 2.;
float ck = (a0 + b0) / 2.;
float ak1 = a0, bk1 = ck, ak2 = ck, bk2 = b0;
//cmp from the smaller side,so can knot the abscissa from the lowest to highest
CmpAdaptiveIntegbySimpson(fun, ak1, bk1, eps, integ, iter, integtab, abscissa_knot, sumErr);
CmpAdaptiveIntegbySimpson(fun, ak2, bk2, eps, integ, iter, integtab, abscissa_knot, sumErr);
}
}
三、结果分析
1、组合梯形公式和辛普森公式实例分析
由实验结果可以看出，当使用插值后再使用组合梯形公式和辛普森公式求解的数值积分结果比直接使用组合梯形公式所得到的结果要更小，并且更接近真实结果。原因可能是直接使用组合梯形公式所得到的结果不具有步长合适的特点，就会导致当步长过大导致结果严重偏离的结果。而拉格朗日插值多项式是对原函数较好的逼近，因此使用拉格朗日插值后所得到的结果更精确。并且，使用插值后的组合辛普森公式比组合梯形公式结果更加精确。但因为拉格朗日插值具有计算量大的特点，因此也会增加计算的难度与复杂度，因此在实际使用过程中应当根据所给点的步长来判断是否使用拉格朗日插值先进行近似逼近。
2、递归公式和龙贝格积分实例分析
由表格可以看出，使用龙贝格积分求解的积分近似值在前两列就基本已经收敛了，说明递归梯形公式和递归辛普森公式已经能够满足精度需求，但是为了方便起见，我们选取最后一行最后一列的值作为最终的结果。当 x x x取 0.5 , 1.0 , 1.5 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0 0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0时，其误差如表（10）所示，因此可以看出使用龙贝格积分的近似程度是特别好的。通过选取阶数M的大小，可以求得任意阶任意步长的积分近似值，然而 M M M并不需要选择太大，因为由表格可以看出龙贝格积分的收敛速度是非常快的，因此可以适当减小M的值，以此来简化计算复杂度提高算力，并且将误差控制在一定范围内。
3、自适应积分实例分析
由六组实验结果可以看出，当使用自适应积分时，可以在不必要使用精细步长的地方增大步长，从而减少计算量；还可以在变化较大的地方适当减小步长以获得更高的精度。最后将所有积分值相加得到最终的积分结果，并且可以看出所有误差和也小于初始误差，因此可以通过设置初始误差来获得任意精度的积分值。初始误差越小，获得的积分值越精确，计算次数更多；反之，初始误差越大，获得的积分值越粗糙，计算次数更少。通过6幅实验结果的图表可以更加直观地看出在函数斜率变化较大地方，其区间会划分的更加细，从而看出自适应积分法的智能性以及实用性。
四、结论
通过分析三种基本类型的数值积分算法，我们可以分析其各自的优缺点以及适用场景。对于组合梯形公式和辛普森公式，其中适合于可以等距取样的情况下，且辛普森公式比梯形公式结果更加精确，步长越小，结果越精确，但是两种方法的精度在较为苛刻的情况下就并不适用了；对于递归公式和龙贝格公式，其适用于几乎所有函数，但是对于算力的消耗比较大，但是其结果一般特别准确，且收敛速度很快，因此可以针对不同情况设置阶数
M
M
M从而尽可能的减少算力、优化算法；对于自适应积分，其适用于变化特别剧烈的函数，如高次幂函数、三角函数、指对数函数等等，适用范围比较广泛，但是对于一般比较平缓的函数其作用并不大。其次，还能通过设置初始误差来得到想要的精度内的结果。
综上所述，三种数值积分方法都有各自的优缺点和具体适用范围，但是其核心思想都是将函数离散化求解积分，并且三种方法都是基于四个基本牛顿-科特斯公式推导出来的，而牛顿-科特斯公式则是由拉格朗日插值多项式推导出来的，因此三种方法都有共同的源头。
声明：本系列《数值计算方法》所有章节源代码均为作者编写，大家可以自行下载，仅供学习交流使用，欢迎大家评论留言或者私信交流，请勿用于任何商业行为。其中，各章实验数据结果也在资源链接中，大家可以自行实验对照，欢迎大家批评指正！文件中每个算法均由作者进行初步测试，部分算法可直接用于工程开发作为模块参考。
2142
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
