复化求积法是一种数值积分的方法, 它通过将积分区间进行分段, 并在各分段子区间上采用低阶的Newton-Cotes求积公式, 对各个小区间上的积分值进行近似, 最后再累加起来. 这种方法避免了随着插值点增多导致的高阶Newton-Cotes公式不稳定的缺陷, 常用的复化求积法有复化梯形求积和复化Simpson求积.
复化梯形公式
复化梯形公式是一种常用的复化求积法, 它通过将积分区间分成若干小区间, 并在每个小区间上使用梯形公式进行近似积分, 最后将所有小区间的积分值相加得到近似积分值.
具体地, 对于一个在 [ a , b ] [a,b] [a,b]上的定积分 ∫ a b f ( x ) d x \displaystyle\int_{a}^{b}f(x){\rm d}x ∫abf(x)dx, 我们可以将其分成 n n n个小区间 [ x i − 1 , x i ] [x_{i-1},x_i] [xi−1,xi], 其中 i = 1 , 2 , . . . , n i=1,2,...,n i=1,2,...,n, 并设 f ( x ) f(x) f(x)在每个小区间的端点处取值分别为 f ( x i − 1 ) f(x_{i-1}) f(xi−1)和 f ( x i ) f(x_i) f(xi). 在每个小区间上, 我们使用梯形公式进行近似积分, 得到 Δ x i \Delta x_i Δxi和 1 2 [ f ( x i − 1 ) + f ( x i ) ] \dfrac12[f(x_{i-1})+f(x_i)] 21[f(xi−1)+f(xi)], 其中 Δ x i = x i − x i − 1 \Delta x_i=x_i-x_{i-1} Δxi=xi−xi−1. 然后我们将所有小区间的积分值相加, 得到近似积分值
∫ a b f ( x ) d x ≈ ∑ i = 1 n 1 2 [ f ( x i − 1 ) + f ( x i ) ] Δ x i \int_{a}^{b}f(x)dx\approx\sum_{i=1}^{n}\frac12[f(x_{i-1})+f(x_i)]\Delta x_i ∫abf(x)dx≈i=1∑n21[f(xi−1)+f(xi)]Δxi
若记
T
n
(
f
)
T_n(f)
Tn(f)为将区间
[
a
,
b
]
[a,b]
[a,b]进行
n
n
n等分割所得的积分值, 则上式可进一步化简得
T
n
(
f
)
=
∑
k
=
0
n
−
1
h
2
[
f
(
x
k
)
+
f
(
x
k
+
1
)
]
=
h
2
[
f
(
a
)
+
f
(
b
)
+
2
∑
k
=
1
n
−
1
f
(
x
k
)
]
T_n(f)=\sum_{k=0}^{n-1}\frac{h}{2}\Big[f(x_k)+f(x_{k+1})\Big]=\frac{h}{2}\Bigg[f(a)+f(b)+2\sum_{k=1}^{n-1}f(x_k)\Bigg]
Tn(f)=k=0∑n−12h[f(xk)+f(xk+1)]=2h[f(a)+f(b)+2k=1∑n−1f(xk)]
其中,
h
:
=
b
−
a
n
h:=\dfrac{b-a}n
h:=nb−a为每个小区间的测度. 特别地, 根据上式, 我们可以得到如下递推式:
T
2
n
(
f
)
=
1
2
T
n
(
f
)
+
h
2
∑
k
=
0
n
−
1
f
(
x
k
+
1
2
)
T_{2n}(f)=\frac{1}{2}T_n(f)+\frac{h}{2}\sum\limits_{k=0}^{n-1}f(x_{k+\frac{1}{2}})
T2n(f)=21Tn(f)+2hk=0∑n−1f(xk+21)
根据中值定理, 复化梯形公式的余项可表示如下:
∃
η
∈
(
a
,
b
)
,
R
T
n
(
f
)
=
I
−
T
n
(
f
)
=
−
b
−
a
12
h
2
f
′
′
(
η
)
\exists\eta\in(a,b),R_{T_n}(f)=I-T_n(f)=-\dfrac{b-a}{12}h^2f^{\prime\prime}(\eta)
∃η∈(a,b),RTn(f)=I−Tn(f)=−12b−ah2f′′(η)
当
f
′
′
(
x
)
f^{\prime\prime}(x)
f′′(x)在区间
[
a
,
b
]
[a,b]
[a,b]上连续, 且其函数值变化不大时, 即有
f
′
′
(
η
n
)
≈
f
′
′
(
η
2
n
)
f^{\prime\prime}(\eta_n)\approx f^{\prime\prime}(\eta_{2n})
f′′(ηn)≈f′′(η2n), 则有
I
−
T
n
(
f
)
I
−
T
2
n
(
f
)
≈
4
\frac{I-T_n(f)}{I-T_{2n}(f)}\approx4
I−T2n(f)I−Tn(f)≈4
故
I
−
T
2
n
(
f
)
≈
1
3
(
T
2
n
(
f
)
−
T
n
(
f
)
)
I-T_{2n}(f)\approx\frac13(T_{2n}(f)-T_n(f))
I−T2n(f)≈31(T2n(f)−Tn(f))
因此, 对于某个允许误差 ε \varepsilon ε, 可以使用 ∣ T 2 n ( f ) − T n ( f ) ∣ < ε |T_{2n}(f)-T_n(f)|<\varepsilon ∣T2n(f)−Tn(f)∣<ε来判断 T 2 n ( f ) T_{2n}(f) T2n(f)是否达到了精度要求.
当被积函数的函数值变化较大时, 为了提高计算精度, 可以在每个小区间的中点处再构造一个梯形公式进行近似积分, 即利用递推式将分割逐渐加细, 从而得到更高精度的积分.
复化Simpson公式
复化Simpson公式与复化梯形公式类似, 但每个小区间上要取3个点, 因此观测点数量需为大于1的奇数. 下面仅考虑分割为等分分割的情况.
将积分区间分为
n
n
n等分, 分点为
x
k
=
a
+
k
h
,
k
=
0
,
1
,
⋯
,
n
x_k=a+kh,k=0,1,\cdots,n
xk=a+kh,k=0,1,⋯,n, 其中,
h
:
=
b
−
a
n
h:=\dfrac{b-a}n
h:=nb−a. 记区间
[
x
k
,
x
k
+
1
]
[x_k,x_{k+1}]
[xk,xk+1]的中点为
x
k
+
1
2
x_{k+\frac12}
xk+21, 在每个小区间上应用Simpson公式, 有
S
n
(
f
)
=
h
6
[
f
(
a
)
+
2
∑
k
=
1
n
−
1
f
(
x
k
)
+
4
∑
k
=
0
n
−
1
f
(
x
k
+
1
2
)
+
f
(
b
)
]
S_n(f)=\frac{h}{6}\bigg[f(a)+2\sum_{k=1}^{n-1}f(x_k)+4\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})+f(b)\bigg]
Sn(f)=6h[f(a)+2k=1∑n−1f(xk)+4k=0∑n−1f(xk+21)+f(b)]
与复化梯形公式类似, 复化Simpson公式也有如下递推公式, 余项估计及允许误差估计如下:
S
2
n
(
f
)
=
1
2
S
n
(
f
)
+
h
6
∑
k
=
0
n
−
1
[
4
f
(
x
k
+
1
4
)
−
2
f
(
x
k
+
1
2
)
+
4
f
(
x
k
+
3
4
)
]
S_{2n}(f)=\frac12S_n(f)+\frac h6\sum_{k=0}^{n-1}\left[4f(x_{k+\frac14})-2f(x_{k+\frac12})+4f(x_{k+\frac34})\right]
S2n(f)=21Sn(f)+6hk=0∑n−1[4f(xk+41)−2f(xk+21)+4f(xk+43)]
∃
η
∈
(
a
,
b
)
,
R
S
n
(
f
)
=
I
−
S
n
(
f
)
=
−
b
−
a
2880
h
4
f
(
4
)
(
η
)
\exists\eta\in(a,b),R_{S_n}(f)=I-S_n(f)=-\dfrac{b-a}{2880}h^4f^{(4)}(\eta)
∃η∈(a,b),RSn(f)=I−Sn(f)=−2880b−ah4f(4)(η)
I
−
S
2
n
(
f
)
≈
1
15
(
S
2
n
(
f
)
−
S
n
(
f
)
)
I-S_{2n}(f)\approx\frac1{15}(S_{2n}(f)-S_n(f))
I−S2n(f)≈151(S2n(f)−Sn(f))
因此同样可以利用上式逐步对区间的分割进行加细, 从而得到更高精度的积分.
算法实现
首先需要包含头文件:
#include <functional>
#include <math.h>
using namespace std;
复化梯形公式
首先实现了一般积分的复化梯形公式如下:
/*
* 复化梯形公式
* f :被积函数
* a :积分下限
* b :积分上限
* e :精度
*/
double composite_trapezoid(function<double(double)> f, const double &a, const double &b, const double &e = 1e-6)
{
if (a >= b)
throw "区间上限小于区间下限！";
double I(f(a));
if (isnan(I) || isnan(I += f(b)))
throw "区间内存在奇点！";
double step(b - a), pre_I(I *= step /= 2);
if (isnan((I /= 2) += step * f(a + step)))
throw "区间内存在奇点！";
while (abs(I - pre_I) > e)
{
double h(step), x(a + (step /= 2)), s(f(x));
if (isnan(s))
throw "区间内存在奇点！";
while ((x += h) < b)
if (isnan(s += f(x)))
throw "区间内存在奇点！";
pre_I = I;
(I /= 2) += step * s;
}
return I;
}
考虑到某些积分上下限可能为奇点, 在上面程序的基础上提供了手动输入边界函数值的实现:
/*
* 复化梯形公式(端点可能为奇点)
* f :被积函数
* a :积分下限
* b :积分上限
* fa:积分下限处的函数值
* fb:积分上限处的函数值
* e :精度
*/
double composite_trapezoid(function<double(double)> f, const double &a, const double &b, const double &fa, const double &fb, const double &e = 1e-6)
{
if (a >= b)
throw "区间上限小于区间下限！";
double step(b - a), I(fa + fb), pre_I(I *= step /= 2);
if (isnan((I /= 2) += step * f(a + step)))
throw "区间内存在奇点！";
while (abs(I - pre_I) > e)
{
double h(step), x(a + (step /= 2)), s(f(x));
if (isnan(s))
throw "区间内存在奇点！";
while ((x += h) < b)
if (isnan(s += f(x)))
throw "区间内存在奇点！";
pre_I = I;
(I /= 2) += step * s;
}
return I;
}
复化Simpson公式
/*
* 复化Simpson公式(端点可能为奇点)
* f :被积函数
* a :积分下限
* b :积分上限
* fa:积分下限处的函数值
* fb:积分上限处的函数值
* e :精度
*/
double composite_Simpson(function<double(double)> f, const double &a, const double &b, const double &fa, const double &fb, const double &e = 1e-6)
{
if (a >= b)
throw "区间上限小于区间下限！";
double step(b - a), m(a + (step /= 2)), fx(f(m));
if (isnan(fx))
throw "区间内存在奇点！";
double I((fx * 4 + fa + fb) / 6), pre_I(I), pre_fx(fx), h(step);
fx = f(a + (step /= 2));
if (isnan((I /= 2) += ((fx += f(b - step)) * 2 - pre_fx) * h / 3))
throw "区间内存在奇点！";
while (abs(I - pre_I) > e)
{
pre_fx = fx, pre_I = I, h = step;
if (isnan(fx = f(m = a + (step /= 2))))
throw "区间内存在奇点！";
while ((m += h) < b)
if (isnan(fx += f(m)))
throw "区间内存在奇点！";
(I /= 2) += (fx * 2 - pre_fx) * h / 3;
}
return I;
}
/*
* 复化Simpson公式
* f :被积函数
* a :积分下限
* b :积分上限
* e :精度
*/
inline double composite_Simpson(function<double(double)> f, const double &a, const double &b, const double &e = 1e-6)
{
double fa(f(a));
if (isnan(fa))
throw "区间内存在奇点！";
double fb(f(b));
if (isnan(fb))
throw "区间内存在奇点！";
return composite_Simpson(f, a, b, fa, fb, e);
}
测试案例
用复化求积法计算 I = ∫ 0 1 sin x x d x \displaystyle I=\int_0^1\frac{\sin x}x{\rm d}x I=∫01xsinxdx的近似值.
#include <stdio.h>
// 算法实现代码略
int main()
{
auto f = [](double x) -> double
{ return sin(x) / x; };
try
{
double I1(composite_trapezoid(f, 0., 1., 1., sin(1.), 1e-8)), I2(composite_Simpson(f, 0., 1., 1., sin(1.), 1e-8));
printf("%f\t%f", I1, I2);
}
catch (const char *s)
{
printf(s);
}
return 0;
}
最终求得 I ≈ 0.946083 I\approx0.946083 I≈0.946083. 另外, 对程序运行进行计时得复化梯形求积法耗时0.172913 μ \mu μs, 而复化Simpson求积法耗时仅为0.00302124 μ \mu μs, 可见复化Simpson求积法在本例中效率较高.
615
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
