\qquad
实际问题当中常常计算积分,有些数值方法如微分方程和积分方程的求解,也都和积分计算相联系。使用牛顿-莱布尼兹(
N
e
w
t
o
n
−
L
e
i
b
n
i
z
Newton-Leibniz
Newton−Leibniz)公式求解往往比较困难,有时原函数不能用初等函数表示或者求解过程十分困难,因此有必要研究积分的数值计算问题。
\qquad
由积分中值定理,在积分区间
[
a
,
b
]
[a,b]
[a,b]中存在一点
ξ
\xi
ξ使得
∫
b
a
f
(
x
)
d
x
=
(
b
−
a
)
f
(
ξ
)
\int_{b}^{a}f(x)dx=(b-a)f(\xi)
∫baf(x)dx=(b−a)f(ξ)
\qquad
但点
ξ
\xi
ξ的具体位置一般难以确定,称
f
(
ξ
)
f(\xi)
f(ξ)为区间
[
a
,
b
]
[a,b]
[a,b]上的平均高度,只要对平均高度
f
(
ξ
)
f(\xi)
f(ξ)提出一种算法,相应的获得一种数值积分方法。
\qquad
如果将两端的高度
f
(
a
)
、
f
(
b
)
f(a)、f(b)
f(a)、f(b)的算术平均值作为平均高度
f
(
ξ
)
f(\xi)
f(ξ)的近似值,则求积公式为
∫
a
b
f
(
x
)
d
x
≈
b
−
a
2
[
f
(
a
)
+
f
(
b
)
]
\int_{a}^{b}f(x)dx\approx\frac{b-a}{2}[f(a)+f(b)]
∫abf(x)dx≈2b−a[f(a)+f(b)]称为梯形公式(n=1)。若使用区间中间点
c
=
a
+
b
2
c=\frac{a+b}{2}
c=2a+b作为
ξ
\xi
ξ,则可得公式为
∫
a
b
f
(
x
)
d
x
≈
(
b
−
a
)
f
(
a
+
b
2
)
\int_{a}^{b}f(x)dx\approx(b-a)f(\frac{a+b}{2})
∫abf(x)dx≈(b−a)f(2a+b)称为中矩形公式或矩形公式。
\qquad
设将积分区间
[
a
,
b
]
[a,b]
[a,b]划分为
n
n
n等份,步长
h
=
b
−
a
n
h=\frac{b-a}{n}
h=nb−a,选去等距节点
x
k
=
a
+
k
h
x_{k}=a+kh
xk=a+kh构造的插值型求积公式为
!
!
!
\color{#F00}{!!!}
!!!
I
n
=
(
b
−
a
)
∑
k
=
0
n
C
k
(
n
)
f
(
x
k
)
I_{n}=(b-a)\sum_{k=0}^{n}C_{k}^{(n)}f(x_{k})
In=(b−a)k=0∑nCk(n)f(xk)称为牛顿-柯特斯(Newton-Cotes)公式,式中
C
k
(
n
)
C_{k}^{(n)}
Ck(n)称为柯特斯系数,引进变换
x
=
a
+
t
h
x=a+th
x=a+th,则有
C
k
(
n
)
=
(
−
1
)
n
−
k
n
k
!
(
n
−
k
)
!
∫
0
n
∏
j
=
0
j
≠
k
n
(
t
−
j
)
d
t
C_{k}^{(n)}=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_{0}^{n}\mathop{\mathop{\prod}\limits_{j=0}}\limits_{j\neq k}^{n}(t-j)dt
Ck(n)=nk!(n−k)!(−1)n−k∫0nj̸=kj=0∏n(t−j)dt
\qquad
辛普森(Simpson)公式(n=2)
I
n
=
(
b
−
a
)
[
1
6
f
(
a
)
+
4
6
f
(
a
+
b
2
)
+
1
6
f
(
b
)
]
I_{n}=(b-a)[\frac{1}{6}f(a)+\frac{4}{6}f(\frac{a+b}{2})+\frac{1}{6}f(b)]
In=(b−a)[61f(a)+64f(2a+b)+61f(b)]代数精确度为3,余项表示为
R
[
f
]
=
K
f
(
4
)
(
η
)
,
η
∈
(
a
,
b
)
R[f]=Kf^{(4)}(\eta),\ \eta\in(a,b)
R[f]=Kf(4)(η), η∈(a,b)
\qquad
更一般的,在区间
[
a
,
b
]
[a,b]
[a,b]上适当选取某些节点
x
k
x_{k}
xk,然后用
f
(
x
k
)
f(x_{k})
f(xk)的加权平均的到平均高度的
f
(
ξ
)
f(\xi)
f(ξ)的近似值,构造的求积公式如下
∫
a
b
f
(
x
)
d
x
≈
∑
k
=
0
n
A
k
f
(
x
k
)
\int_{a}^{b}f(x)dx\approx\sum_{k=0}^{n}A_{k}f(x_{k})
∫abf(x)dx≈k=0∑nAkf(xk)式中
x
k
x_{k}
xk称为求积节点,
A
k
A_{k}
Ak称为求积系数,亦称为伴随节点
x
k
x_{k}
xk的权。权
A
k
A_{k}
Ak仅与节点
x
k
x_{k}
xk的选择有关,不依赖于被积函数
f
(
s
)
f(s)
f(s)的具体形式。
\qquad
代数精度 定义1
!
!
!
\color{#F00}{!!!}
!!!如果某个求积公式对于次数不超过
m
m
m的多项式均能准确的成立,但对于
m
+
1
m+1
m+1次多项式就不准确成立,则称该求积公式具有
m
m
m次代数精度。欲使求积公式具有
m
m
m次代数精度,只要令其对
f
(
x
)
=
1
,
x
,
⋯
 
,
x
m
f(x)=1,x,\cdots,x^{m}
f(x)=1,x,⋯,xm都能准确成立,即
{
∑
A
k
=
b
−
a
∑
A
k
x
k
=
1
2
(
b
2
−
a
2
)
,
⋮
∑
A
k
x
k
m
=
1
m
+
1
(
b
m
+
1
−
a
m
+
1
)
\begin{cases} \sum A_{k}=b-a\\ \sum A_{k}x_{k}=\frac{1}{2}(b^{2}-a^{2}),\\ \vdots \\ \sum A_{k}x_{k}^{m}=\frac{1}{m+1}(b^{m+1}-a^{m+1}) \end{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧∑Ak=b−a∑Akxk=21(b2−a2),⋮∑Akxkm=m+11(bm+1−am+1)
\qquad
求积公式的收敛性与稳定性 定义2在求积公式中,若
l
i
m
n
→
∞
h
→
0
∑
k
=
0
n
A
k
f
(
x
k
)
−
∫
a
b
f
(
x
)
d
x
\mathop{\mathop{lim}\limits_{n\to\infty}}\limits_{h\to 0}\sum_{k=0}^{n}A_{k}f(x_{k})-\int_{a}^{b}f(x)dx
h→0n→∞limk=0∑nAkf(xk)−∫abf(x)dx其中
h
=
m
a
x
1
≤
i
≤
n
{
x
i
−
x
i
−
1
}
h=\mathop{max}\limits_{1\leq i\leq n}\{x_{i}-x_{i-1}\}
h=1≤i≤nmax{xi−xi−1},则称求积公式收敛。在求积公式中,由于计算
f
(
x
k
)
f(x_{k})
f(xk)可能产生误差
δ
k
\delta_{k}
δk实际的得到的为
f
k
~
\tilde{f_{k}}
fk~,即
f
(
x
k
)
=
f
k
~
+
δ
k
f(x_{k})=\tilde{f_{k}}+\delta_{k}
f(xk)=fk~+δk,记
I
n
(
f
)
=
∑
k
=
0
n
A
k
f
(
x
k
)
I
n
(
f
~
)
=
∑
k
=
0
n
A
k
f
k
~
I_{n}(f)=\sum_{k=0}^{n}A_{k}f(x_{k})\quad I_{n}(\tilde{f})=\sum_{k=0}^{n}A_{k}\tilde{f_{k}}
In(f)=k=0∑nAkf(xk)In(f~)=k=0∑nAkfk~如果对任给小正数
ε
>
0
\varepsilon>0
ε>0,只要误差
∣
δ
k
∣
|\delta_{k}|
∣δk∣充分小就有,
∣
I
n
(
f
)
−
I
n
(
f
~
)
∣
=
∣
∑
k
=
0
n
A
k
[
f
(
x
k
)
−
f
k
~
]
∣
≤
ε
|I_{n}(f)-I_{n}(\tilde{f})|=|\sum_{k=0}^{n}A_{k}[f(x_{k})-\tilde{f_{k}}]|\leq\varepsilon
∣In(f)−In(f~)∣=∣k=0∑nAk[f(xk)−fk~]∣≤ε
表
明
求
积
公
式
计
算
是
稳
定
的
\color{#F00}{表明求积公式计算是稳定的}
表明求积公式计算是稳定的。
\qquad
定义3 对任意给定的
ε
>
0
\varepsilon>0
ε>0,若
∃
δ
>
0
\exists\delta>0
∃δ>0,只要
∣
f
(
x
k
)
−
f
k
~
∣
≤
δ
(
k
=
0
,
1
,
⋯
 
,
n
)
|f(x_{k})-\tilde{f_{k}}|\leq\delta(k=0,1,\cdots,n)
∣f(xk)−fk~∣≤δ(k=0,1,⋯,n)就有上式成立,则称求积公式是稳定的。
\qquad
求积公式的余项 若求积公式的代数精度为
m
m
m,则由求积公式余项的表达式可将余项表达为
R
[
f
]
=
∫
a
b
f
(
x
)
d
x
−
∑
k
=
0
n
A
k
f
(
x
k
)
=
K
f
(
m
+
1
)
(
η
)
R[f]=\int_{a}^{b}f(x)dx-\sum_{k=0}^{n}A_{k}f(x_{k})=Kf^{(m+1)}(\eta)
R[f]=∫abf(x)dx−k=0∑nAkf(xk)=Kf(m+1)(η)
\qquad
K
K
K为不依赖
f
(
x
)
f(x)
f(x)的待定参数,
η
∈
(
a
,
b
)
\eta\in(a,b)
η∈(a,b)。结果表明当
f
(
x
)
f(x)
f(x)是次数小于等于
m
m
m的多项式时,由于
f
(
m
+
1
)
(
x
)
=
0
f^{(m+1)}(x)=0
f(m+1)(x)=0,故此时
R
[
f
]
=
0
R[f]=0
R[f]=0,即求积公式精确成立。当
f
(
x
)
=
x
m
+
1
f(x)=x^{m+1}
f(x)=xm+1时,
f
(
m
+
1
)
(
x
)
=
(
m
+
1
)
!
f^{(m+1)}(x)=(m+1)!
f(m+1)(x)=(m+1)!,又
R
n
(
f
)
≠
0
R_{n}(f)\neq 0
Rn(f)̸=0,可求得
K
=
1
(
m
+
1
)
!
[
∫
a
b
x
m
+
1
d
x
−
∑
k
=
0
n
A
k
x
k
m
+
1
]
=
1
(
m
+
1
)
!
[
1
m
+
2
(
b
m
+
2
−
a
m
+
2
)
−
∑
k
=
0
n
A
k
x
k
m
+
1
]
K=\frac{1}{(m+1)!}[\int_{a}^{b}x^{m+1}dx-\sum_{k=0}^{n}A_{k}x_{k}^{m+1}]\\ =\frac{1}{(m+1)!}[\frac{1}{m+2}(b^{m+2}-a^{m+2})-\sum_{k=0}^{n}A_{k}x_{k}^{m+1}]
K=(m+1)!1[∫abxm+1dx−k=0∑nAkxkm+1]=(m+1)!1[m+21(bm+2−am+2)−k=0∑nAkxkm+1]带入余项可以的到具体的余项表达式。梯形公式和中矩形公式的代数精度为1,其余项表达式为
R
[
f
]
=
K
f
′
′
(
η
)
R[f]=Kf^{''}(\eta)
R[f]=Kf′′(η)
复合求积公式
\qquad
为了提高精度通常将积分区域分成若干子区间(通常等分),在每个子区间上用低阶求积公式,称为复合求积法。
\qquad
复合梯形公式 将区间
[
a
,
b
]
[a,b]
[a,b]分为
n
n
n等份,分点
x
k
=
a
+
k
h
,
h
=
b
−
a
n
,
k
=
0
,
1
,
⋯
 
,
n
x_{k}=a+kh,h=\frac{b-a}{n},k=0,1,\cdots,n
xk=a+kh,h=nb−a,k=0,1,⋯,n在每个子区间上
[
x
k
,
x
k
+
1
]
(
k
=
0
,
1
,
⋯
 
,
n
−
1
)
[x_{k},x_{k+1}](k=0,1,\cdots,n-1)
[xk,xk+1](k=0,1,⋯,n−1)上采用梯形公式则得
I
=
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
−
1
∫
x
k
x
k
+
1
f
(
x
)
d
x
=
h
2
∑
k
=
0
n
−
1
[
f
(
x
k
+
f
(
x
k
+
1
)
)
]
+
R
n
(
f
)
I=\int_{a}^{b}f(x)dx=\sum_{k=0}^{n-1}\int_{x_{k}}^{x_{k+1}}f(x)dx=\frac{h}{2}\sum_{k=0}^{n-1}[f(x_{k}+f(x_{k+1}))]+R_{n}(f)
I=∫abf(x)dx=k=0∑n−1∫xkxk+1f(x)dx=2hk=0∑n−1[f(xk+f(xk+1))]+Rn(f)记
T
n
=
h
2
∑
k
=
0
n
−
1
[
f
(
x
k
+
f
(
x
k
+
1
)
)
]
T_{n}=\frac{h}{2}\sum_{k=0}^{n-1}[f(x_{k}+f(x_{k+1}))]
Tn=2hk=0∑n−1[f(xk+f(xk+1))]称为复合梯形公式。其余项表示为
R
n
(
f
)
=
I
−
T
n
=
∑
k
=
0
n
−
1
[
−
h
3
12
f
′
′
(
η
k
)
]
,
η
k
∈
(
x
k
,
x
k
+
1
)
R_{n}(f)=I-T_{n}=\sum_{k=0}^{n-1}[-\frac{h^{3}}{12}f^{''}(\eta_{k})],\ \eta_{k}\in(x_{k},x_{k+1})
Rn(f)=I−Tn=k=0∑n−1[−12h3f′′(ηk)], ηk∈(xk,xk+1)复合辛普森求积公式与此相似,子区间上展开方式变为辛普森公式。
\qquad
定理 设
f
(
x
)
∈
C
∞
[
a
,
b
]
f(x)\in C^{\infty}\ [a,b]
f(x)∈C∞ [a,b],则有
T
(
h
)
=
I
+
α
1
h
2
+
α
2
h
4
+
⋯
+
α
l
h
2
l
+
⋯
T(h)=I+\alpha_{1}h^{2}+\alpha_{2}h^{4}+\cdots+\alpha_{l}h^{2l}+\cdots
T(h)=I+α1h2+α2h4+⋯+αlh2l+⋯其中系数
α
l
(
l
=
1
,
2
,
⋯
 
)
\alpha_{l}(l=1,2,\cdots)
αl(l=1,2,⋯)与
h
h
h无关(由泰勒展开可得)。只要真值与近似值得误差能表示为
h
h
h的幂级数,都能使用外推算法(将步长变为
1
2
h
\frac{1}{2}h
21h),提高精度。外推算法可写为统一形式(龙贝格算法):
T
m
(
h
)
=
4
m
4
m
−
1
T
m
−
1
(
h
2
)
−
1
4
m
−
1
T
m
−
1
(
h
)
T_{m}(h)=\frac{4^{m}}{4^{m}-1}T_{m-1}(\frac{h}{2})-\frac{1}{4^{m}-1}T_{m-1}(h)
Tm(h)=4m−14mTm−1(2h)−4m−11Tm−1(h)
自适应积分方法
\qquad 在复合求积过程中,当被积函数变化剧烈时,使用较小步长,当被积函数变化平缓时,使用较大步长。