数值积分（梯形法）

贡献者：待更新

预备知识　Matlab 的函数，定积分

　　本文介绍一种简单的梯形算法计算数值积分。如图 1 若要用梯形法计算定积分 $\int_a^b f(x) \,\mathrm{d}{x} $，则可将区间 $[a, b]$ 划分为 $N$ 个长度为 $h = (b-a)/N$ 的等长的小区间，区间端点从 $a$ 到 $b$ 分别为 $x_1 = a, x_2 = a + h, \dots, x_{N+1} = b$。

图 1：梯形法数值积分

　　接下来将每个区间的被曲线围出的面积用梯形来计算，第 $i$ 个梯形面积为（两底和乘以高除以二） $[f(x_{i+1}) + f(x_i)]h/2$。定积分约等于所有梯形面积

\begin{equation} \int_a^b f(x) \,\mathrm{d}{x} \approx \sum_i^N \frac h2 [f(x_{i+1}) + f(x_i)] = h \left[\frac12 f(x_1) + \sum_{i = 2}^N f(x_i) + \frac12 f(x_{N+1}) \right] ~. \end{equation}

显然，当 $N$ 取越大时，右边的求和就越接近定积分。

　　这里给出 Matlab 代码

代码 1：trapezoidInt.m

% 梯形法数值积分
% f 是被积函数的函数句柄
% [a, b] 为积分区间
% N 为子区间个数
function I = trapezoidInt(f, a, b, N)
x = linspace(a, b, N+1);
y = arrayfun(f, x); % 求所有 y(i) = f(x(i))
h = (b - a)/N;
I = h*(0.5*y(1) + 0.5*y(N+1) + sum(y(2:N)));
end

　　下面来看两个例子，由算法可知，当被积函数具有 $c_1 x + c_2$ 的形式时，函数曲线围成的面积可由梯形精确解算，所以即使 N 取很小也可以精确到 double 类型的精度（16 位左右）。

>> format long;
>> trapezoidInt(@(x)x, 0, 1, 4)
ans = 0.500000000000000

对其他一些函数，需要较大的 $N$ 才能获得较高的精度。我们已知 $\sin x$ 在 $[0, \pi]$ 内的积分等于 $2$，再看数值积分

>> trapezoidInt(@sin, 0, pi, 10)
ans = 1.983523537509455
>> trapezoidInt(@sin, 0, pi, 100)
ans = 1.999835503887444
>> trapezoidInt(@sin, 0, pi, 1000)
ans = 1.999998355065662
>> format short;

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。