数值积分(梯形法)

                     

贡献者: 待更新

预备知识 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;

                     

© 小时科技 保留一切权利