数值积分：通过复合梯形法计算
在物理学和工程学中，很多问题都可以通过数值积分来求解，特别是当我们无法得到解析解时。数值积分是通过计算积分区间内离散点的函数值来近似积分的结果。在这篇博客中，我将讨论如何使用 复合梯形法 来进行数值积分，并以一个简单的计算为例，展示如何使用该方法，并附带Matlab程序。
数值积分的基本原理
数值积分的基本思想是将一个连续的积分问题转化为离散的求和问题，通常是通过将积分区间划分成若干小的子区间。
复合梯形法
复合梯形法就是其中的一种方法，它的基本原理如下：
- 首先，将积分区间分为 n 个小区间，每个小区间的宽度为 h。
- 然后，利用梯形公式近似每个小区间内函数的积分：
将积分区间
[
a
,
b
]
[a, b]
[a,b] 划分为
n
n
n 等分，则每个小区间的宽度为：
h
=
b
−
a
n
h = \frac{b - a}{n}
h=nb−a
因此，积分可以表示为：
I
=
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
−
1
∫
x
k
x
k
+
1
f
(
x
)
d
x
I = \int_a^b f(x) \, dx = \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} f(x) \, dx
I=∫abf(x)dx=k=0∑n−1∫xkxk+1f(x)dx
每个子区间
[
x
k
,
x
k
+
1
]
[x_k, x_{k+1}]
[xk,xk+1] 上的积分使用梯形公式，得到：
∫
x
k
x
k
+
1
f
(
x
)
d
x
≈
h
2
[
f
(
x
k
)
+
f
(
x
k
+
1
)
]
\int_{x_k}^{x_{k+1}} f(x) \, dx \approx \frac{h}{2} \left[ f(x_k) + f(x_{k+1}) \right]
∫xkxk+1f(x)dx≈2h[f(xk)+f(xk+1)]
因此，整个区间的积分为：
I
=
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
n
−
1
∫
x
k
x
k
+
1
f
(
x
)
d
x
≈
∑
k
=
0
n
−
1
h
2
[
f
(
x
k
)
+
f
(
x
k
+
1
)
]
=
h
2
{
[
f
(
x
0
)
+
f
(
x
1
)
]
+
[
f
(
x
1
)
+
f
(
x
2
)
]
+
…
+
[
f
(
x
n
−
2
)
+
f
(
x
n
−
1
)
]
+
[
f
(
x
n
−
1
)
+
f
(
x
n
)
]
}
=
h
2
[
f
(
x
0
)
+
2
∑
k
=
1
n
−
1
f
(
x
k
)
+
f
(
x
n
)
]
=
h
2
[
f
(
a
)
+
f
(
b
)
]
+
h
⋅
∑
i
=
1
n
−
1
f
(
x
k
)
\begin{aligned} I&= \int_a^b f(x) dx = \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} f(x) \, dx \\ &\approx \sum_{k=0}^{n-1} \frac{h}{2} \left[ f(x_k) + f(x_{k+1}) \right] \\ &= \frac{h}{2} \{ [ f(x_0) + f(x_1)] + [f(x_1) + f(x_2)] +\dots \\ &+ [f(x_{n-2}) + f(x_{n-1})] + [f(x_{n-1}) + f(x_{n})] \} \\ &= \frac{h}{2} \left[ f(x_0) + 2\sum_{k=1}^{n-1} f(x_k) + f(x_n) \right] \\ &= \frac{h}{2} \left[ f(a)+f(b) \right] + h\cdot \sum_{i=1}^{n-1}f(x_k) \\ \end{aligned}
I=∫abf(x)dx=k=0∑n−1∫xkxk+1f(x)dx≈k=0∑n−12h[f(xk)+f(xk+1)]=2h{[f(x0)+f(x1)]+[f(x1)+f(x2)]+…+[f(xn−2)+f(xn−1)]+[f(xn−1)+f(xn)]}=2h[f(x0)+2k=1∑n−1f(xk)+f(xn)]=2h[f(a)+f(b)]+h⋅i=1∑n−1f(xk)
其中：
•
x
i
x_i
xi 是区间内的离散点：
x
i
=
a
+
i
⋅
h
x_i = a + i \cdot h
xi=a+i⋅h 。
•
f
(
x
)
f(x)
f(x) 是被积函数。
最后一项称为复化梯形公式，记为：
T
n
=
h
2
[
f
(
a
)
+
2
∑
k
=
1
n
−
1
f
(
x
k
)
+
f
(
b
)
]
T_n = \frac{h}{2} \left[ f(a) + 2 \sum_{k=1}^{n-1} f(x_k) + f(b) \right]
Tn=2h[f(a)+2k=1∑n−1f(xk)+f(b)]
误差：若 f ( x ) ∈ C 2 [ a , b ] f(x) \in C^2[a, b] f(x)∈C2[a,b]，其求积余项
R n ( f ) = − b − a 12 h 2 f ′ ′ ( η ) , 其中 η ∈ [ a , b ] R_n(f) = -\frac{b-a}{12} h^2 f^{\prime\prime}(\eta), \ 其中 \eta \in [a, b] Rn(f)=−12b−ah2f′′(η), 其中η∈[a,b] 可以看出误差是 h 2 h^2 h2阶，复化梯形公式是收敛的。此外，T_n 的求积系数均为正，由定理 2 知复化梯形公式是稳定的。
复合梯形法的实现步骤
- 计算区间两端的函数值：这是梯形法中的基础项，计算 f ( a ) f(a) f(a) 和 f ( b ) f(b) f(b) 。
- 计算区间中间点的函数值：计算区间内的所有点 f ( x 1 ) , f ( x 2 ) , … , f ( x n − 1 ) f(x_1), f(x_2), \dots, f(x_{n-1}) f(x1),f(x2),…,f(xn−1)。
- 应用梯形公式：将所有小区间的结果加总，得到整个积分的近似值。
复合梯形法的代码实现 (Matlab)
例如，考虑函数
f
(
x
)
=
ln
∣
tanh
(
x
)
∣
2
f(x) = \ln |\tanh(x)|^2
f(x)=ln∣tanh(x)∣2 在区间
[
0
,
2
π
]
[0, 2\pi]
[0,2π] 上的积分。使用复合梯形法来计算该积分：
∫
0
2
π
ln
∣
tanh
(
x
)
∣
2
d
x
≈
h
2
[
ln
∣
tanh
(
0
)
∣
2
+
2
∑
i
=
1
n
−
1
ln
∣
tanh
(
x
i
)
∣
2
+
ln
∣
tanh
(
2
π
)
∣
2
]
\begin{aligned} \int_0^{2\pi} \ln |\tanh(x)|^2 \, dx & \approx \frac{h}{2} \left[ \ln |\tanh(0)|^2 + 2 \sum_{i=1}^{n-1} \ln |\tanh(x_i)|^2 + \ln |\tanh(2\pi)|^2 \right] \end{aligned}
∫02πln∣tanh(x)∣2dx≈2h[ln∣tanh(0)∣2+2i=1∑n−1ln∣tanh(xi)∣2+ln∣tanh(2π)∣2]
下面是使用复合梯形法实现数值积分的 MATLAB 代码：
% 定义积分区间
a = 0; % 积分下限
b = 2*pi; % 积分上限
n = 100; % 将区间分成 100 个小区间
step = (b - a) / n; % 每个小区间的宽度
x = a:step:b; % 生成区间内的离散点
% 定义被积函数
f = @(x) log(abs(tanh(x)))^2; % 这里使用一个简单的示例函数
% 复合梯形法积分
T0 = (step/2) * (f(a) + f(b)); % 计算区间两端的函数值
T1 = 0;
for i = 1:(n-1)
T1 = T1 + (step) * f(x(i+1)); % 求和中间点的函数值
end
integral_result = T0 + T1; % 得到积分结果
% 输出结果
disp(['积分结果：', num2str(integral_result)]);