MATLAB - 梯形法则

MATLAB - 使用梯形法则进行数值积分

数值积分（Numerical integration）是一种近似计算函数定积分的基本技术，尤其在难以或无法找到解析解时非常有用。**梯形法则（trapezoidal rule）**是一种流行的数值方法，它通过将曲线下的面积分割成一系列梯形并对它们的面积求和来近似计算。这种方法通常比使用矩形法更准确。

梯形法则背后的理论

步骤 1：划分区间 [a, b]

首先，我们将积分区间 [a, b] 划分为 n 个更小的子区间，每个子区间宽度相等，为 Δx。

宽度 Δx 计算如下： $$ \Delta x = \frac{b - a}{n} $$

这将创建 n+1 个点：x₀, x₁, x₂, …, xₙ，其中：

$$ x_i = a + i \cdot \Delta x \quad \text{for } i = 0, 1, 2, \dots, n $$

这里，x₀ = a 和 xₙ = b 是区间的端点。

步骤 2：形成梯形并计算它们的面积

对于每个子区间 [xᵢ₋₁, xᵢ]，我们形成一个梯形。梯形的平行边是端点处的函数值 f(xᵢ₋₁) 和 f(xᵢ)，其高为 Δx。

单个梯形的面积由标准公式给出：

$$ \text{Area}i = \frac{1}{2} \times (\text{sum of parallel sides}) \times \text{height} = \frac{\Delta x}{2} [f(x{i-1}) + f(x_{i})] $$

步骤 3：求面积和

总近似积分是所有 n 个梯形面积之和：

$$ \int_{a}^{b} f(x) , dx \approx \sum_{i=1}^{n} \text{Area}_i $$

展开这个求和公式会揭示一个更有效的公式。请注意，所有内部点（f(x₁) 到 f(xₙ₋₁)）都被两个相邻的梯形共享，因此它们被计算了两次。

$$ \int_{a}^{b} f(x) , dx \approx \frac{\Delta x}{2} \left[ f(x_0) + 2f(x_1) + 2f(x_2) + \dots + 2f(x_{n-1}) + f(x_n) \right] $$

在 MATLAB 中实现梯形法则

MATLAB 提供了内置函数 trapz，可以直接对离散数据点执行梯形积分。当您拥有实验数据（例如，随时间变化的传感器读数）而不是数学函数时，这非常有用。

语法

Q = trapz(Y)
Q = trapz(X, Y)
Q = trapz(___, dim)

语法解释

Q = trapz(Y)：计算向量 Y 中数据的积分，假设单位间距（即 X = [1, 2, 3, ...]）。如果 Y 是矩阵，trapz 将独立积分每列。

Q = trapz(X, Y)：计算 Y 相对于 X 中指定坐标的积分。X 和 Y 必须是等长的向量。这是非单位间距最常见的用法。

Q = trapz(___, dim)：沿着指定维度 dim 对多维数组进行积分。例如，trapz(X, Y, 2) 沿着矩阵 Y 的行进行积分。

示例 1：对均匀和非均匀间距的数据进行积分

让我们使用离散点计算函数 y = x² 从 x=1 到 x=5 的积分。

% 定义坐标向量和值向量
X = 1:5; % X = [1, 2, 3, 4, 5]
Y = X.^2; % Y = [1, 4, 9, 16, 25]

% 情况 1：使用 trapz(Y) 假设单位间距（此处正确）
Q1 = trapz(Y);

% 情况 2：使用 trapz(X, Y) 明确提供间距
Q2 = trapz(X, Y);

% 解析解为 ∫x² dx 从 1 到 5 = [x³/3] 从 1 到 5 = 125/3 - 1/3 = 124/3 ≈ 41.33
fprintf('Result with assumed unit spacing (Q1): %.4f\n', Q1);
fprintf('Result with explicit spacing (Q2): %.4f\n', Q2);
fprintf('Analytical Result: %.4f\n', 124/3);

% --- 输出 ---
% Result with assumed unit spacing (Q1): 42.0000
% Result with explicit spacing (Q2): 42.0000
% Analytical Result: 41.3333

数值结果 42.0 是真实值 41.333 的一个很好的近似。精度可以通过增加数据点来提高。

示例 2：沿列和行对矩阵进行积分

% 创建一个表示 3 个不同数据集的 3x3 矩阵
Y = [1 2 3;
     4 5 6;
     7 8 9];

% 沿每列积分（dim=1 是默认值）
Q_cols = trapz(Y);

% 沿每行积分（指定 dim=2）
Q_rows = trapz(Y, 2);

disp('Original Matrix Y:');
disp(Y);
disp('Integral of each column:');
disp(Q_cols);
disp('Integral of each row:');
disp(Q_rows);

% --- 输出 ---
% Original Matrix Y:
%      1     2     3
%      4     5     6
%      7     8     9
%
% Integral of each column:
%      8    10    12
%
% Integral of each row:
%      4
%     10
%     16

实际应用：从速度数据计算距离

假设一个传感器记录车辆在不同时间点的速度。我们可以对速度进行积分以找出总行驶距离。

% 时间（秒，非均匀间距）
time_s = [0, 10, 25, 40, 60];

% 速度（米/秒）
velocity_mps = [0, 22, 28, 30, 25];

% 通过对速度随时间积分来计算总距离
total_distance_m = trapz(time_s, velocity_mps);

fprintf('Total distance traveled: %.2f meters\n', total_distance_m);

% --- 输出 ---
% Total distance traveled: 1475.00 meters

现代替代方案：integral 函数

对于需要积分数学函数的情况，MATLAB 的 integral 函数是更现代、更准确、更受推荐的选择。它使用一种更复杂的技巧，称为自适应正交（adaptive quadrature）。

% 使用函数句柄定义要积分的函数
fun = @(x) x.^2;

% 定义积分区间
a = 1;
b = 5;

% 计算积分
Q = integral(fun, a, b);

fprintf('The result from integral() is: %.4f\n', Q);

% --- 输出 ---
% The result from integral() is: 41.3333

最佳实践和常见陷阱

• trapz 与 integral 对比：当您拥有离散数据点（例如，来自实验）时，使用 trapz。当您拥有一个需要评估的符号函数时，使用 integral。
• 单调坐标：使用 trapz(X,Y) 时，请确保坐标向量 X 是单调的（始终递增或始终递减）。
• 向量长度：X 和 Y 必须具有相同的元素数量。
• 准确性：梯形法则的准确性取决于数据点的数量。更多点（更小的 Δx）通常会带来更准确的结果。