MATLAB - 辛普森法则

MATLAB - 使用辛普森法则进行数值积分

辛普森法则（Simpson’s Rule）是一种强大的数值方法，用于逼近定积分的值。它通过使用二次多项式来逼近每个子区间内的函数，从而比梯形法则等更简单的方法提供更高的精度。这使得它对于平滑且行为良好的函数特别有效。

该法则最常以辛普森 1/3 法则（Simpson’s 1/3 Rule）的形式应用，它要求将积分区间 [a, b] 分成偶数 (n) 个子区间。

辛普森 1/3 法则的公式为：

$$\int_{a}^{b} f(x) , dx \approx \frac{\Delta x}{3} [f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + \dots + 4f(x_{n-1}) + f(x_{n})]$$

其中 $$\Delta x = \frac{b-a}{n}$$ 和 $$x_i = a + i\Delta x$$。注意权重的模式是：1, 4, 2, 4, ..., 2, 4, 1。

辛普森法则的工作原理

1. 划分区间：将总区间 [a, b] 分成 n 个等长子区间，其中 n 必须为偶数。
2. 拟合二次多项式：在每对相邻的子区间上，拟合一个唯一的二次多项式（即抛物线），使其通过函数的三个点。
3. 求面积和：通过对这些抛物线下的精确面积求和来近似积分。

MATLAB 中的实现

虽然您可以从头开始实现辛普森法则以理解算法，但现代 MATLAB 提供了更高效的方法。我们将探讨三种方法：基本循环实现、向量化实现以及使用 MATLAB 的内置函数。

方法 1：使用 for 循环实现（用于学习）

这种方法直观明了，清晰地展示了算法的逻辑。我们将逼近函数 $$f(x) = x^2$$ 从 a=0 到 b=2 的积分。

% 1. 定义函数、区间和子区间数量
f = @(x) x.^2;
a = 0;
b = 2;
n = 10; % 必须是偶数

% 2. 检查 n 是否为偶数
if mod(n, 2) ~= 0
    error('Number of subintervals (n) must be even for Simpson''s 1/3 rule.');
end

% 3. 计算步长
h = (b - a) / n;

% 4. 用首项和末项初始化和
integral_sum = f(a) + f(b);

% 5. 遍历内部点以添加加权项
for i = 1:(n - 1)
    x_i = a + i * h;
    if mod(i, 2) == 0 % 偶数索引项权重为 2
        integral_sum = integral_sum + 2 * f(x_i);
    else % 奇数索引项权重为 4
        integral_sum = integral_sum + 4 * f(x_i);
    end
end

% 6. 完成计算
approx_integral = (h / 3) * integral_sum;

fprintf('Approximate integral (loop): %.6f\n', approx_integral);

% 为了比较，精确值为 8/3 ≈ 2.666667

输出： Approximate integral (loop): 2.666667

方法 2：向量化实现（最佳实践）

一种更高效、更具“MATLAB风格”的方法是使用向量运算来避免循环。对于较大的 n，这会显著加快速度。

% 1. 定义函数、区间和 n（同前）
f = @(x) x.^2; a = 0; b = 2; n = 10;
if mod(n, 2) ~= 0, error('n must be even'); end

% 2. 创建所有 x 点的向量
x = linspace(a, b, n + 1);

% 3. 计算对应的 y 值
y = f(x);

% 4. 使用向量求和应用辛普森法则
% 对奇数索引的内部点求和（y(2), y(4), ...）
% 对偶数索引的内部点求和（y(3), y(5), ...）
h = (b - a) / n;
approx_integral_vec = (h/3) * (y(1) + 4*sum(y(2:2:n)) + 2*sum(y(3:2:n-1)) + y(n+1));

fprintf('Approximate integral (vectorized): %.6f\n', approx_integral_vec);

输出： Approximate integral (vectorized): 2.666667

方法 3：使用 integral 的现代 MATLAB 方法（推荐）

对于大多数实际应用，您应该使用 MATLAB 的内置 integral 函数。它采用全局自适应求积（globally adaptive quadrature），比辛普森法则等固定步长规则更鲁棒和准确。它自动处理复杂性并提供误差估计。

% 1. 定义函数和区间（同前）
f = @(x) x.^2;
a = 0;
b = 2;

% 2. 调用 integral 函数
matlab_integral = integral(f, a, b);

fprintf('Integral using built-in function: %.6f\n', matlab_integral);

输出： Integral using built-in function: 2.666667

结论：何时使用每种方法

• 循环实现：用于学习算法或需要检查中间步骤以进行教学时。
• 向量化实现：当您需要专门实现辛普森法则（例如，用于课堂作业）并希望获得高效、编写良好的脚本时使用。
• integral 函数：用于所有实际和专业工作。它经过高度优化，更准确，并且可以处理更广泛的函数，包括那些具有奇异点的函数。