MATLAB - 辛普森法则

辛普森法则是一种用于逼近函数定积分的数值方法。它比梯形法则更精确，尤其适用于相当平滑且表现良好的函数。辛普森法则使用二次多项式来逼近每个子区间上的函数，而不是线性逼近。

辛普森法则通常以两种主要形式应用

辛普森 1/3 法则 - 对于将单个区间 [a,b] 分成 n 个等间距子区间（其中 n 为偶数）的情况，积分的近似值由下式给出。

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

$\mathrm{其中 \: \Delta x\: = \: \frac{b-a}{n} \: 且 \: x_{i} \: = \: a \: + \: i\Delta x \: 对于 \: i\:=\: 0,1,2,\dots , n}$

辛普森 3/8 法则 - 对于将单个区间 [a,b] 分成 n 个等间距子区间（其中 n 为 3 的倍数）的情况，近似值由下式给出。

$$\mathrm{\displaystyle\int_{a}^{b} f(x) \: dx \: \approx \: \frac{3 \Delta x}{8} [f(x_{0}) \:+\: 3f(x_{1}) \:+\: 3f(x_{2}) \:+\: \cdots \:+\: 3f(x_{n-1}) \:+\: f(x_{n})]}$$

$\mathrm{其中 \: \Delta x\: = \: \frac{b-a}{n} \: 且 \: x_{i} \: = \: a \: + \: i\Delta x \: 对于 \: i\:=\: 0,1,2,\dots , n}$

辛普森法则如何工作？

划分区间 - 您要对函数 f(x) 进行积分的区间 [a,b] 被分成 n 个子区间。对于辛普森 1/3 法则，n 必须为偶数。

应用二次逼近 - 在每个子区间对上，辛普森法则通过一个二次多项式来逼近函数 f(x)，该多项式拟合区间端点和中点处的函数值。

求和面积 - 曲线下的面积通过对这些二次逼近下的面积求和来近似。

让我们看一个简单的示例，它展示了辛普森法则的工作原理。

示例

让我们使用辛普森 1/3 法则和 n=4 个子区间来逼近 f(x)=x2 从 a=0 到 b=2 的积分。

1. 计算 Δx

$$\mathrm{\Delta x\: = \: \frac{b \: -\:a}{n} \: = \: \frac{2\: - \:0}{4} \: = \: 0.5}$$

2. 确定 xi 值

$$\mathrm{x_{0} \: = \: 0, \: x_{1} \:= \: 0.5, \: x_{2}\: = \:1, \: x_{3} \: = \: 1.5, \: x_{4} \: =\: 2}$$

3. 计算函数值。

$$\mathrm{f(x_{0})\: =\: 0^{2} \: = \: 0, \: f(x_{1}) \: =\: 0.5^{2} \: = \: 0.25, \: f(x_{2})\: = \:1^{2} \: = \: 1, \: f(x_{3}) \: = \: 1.5^{2} \: = \: 2.25, \: f(x_{4}) \: =\: 2^{2} \:=\: 4}$$

4. 应用辛普森 1/3 法则。

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} [f(0) \:+\: 4f(0.5) \:+\: 2f(1) \:+\: 4f(1.5) \:+\: f(2)]}$$

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} [0 \:+\: 4\cdot 0.25 \:+\: 2 \cdot 1 \:+\: 4 \cdot 2.25 \:+\: 4]}$$

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} [0 \:+\: 1 \:+\: 2 \:+\: 9 \:+\: 4]}$$

$$\mathrm{\displaystyle\int_{0}^{2} x^{2} \: dx \: \approx \: \frac{0.5}{3} \:\cdot\: 16 \:=\: \frac{8}{3} \: \approx \: 2.6667 }$$

x2 从 0 到 2 的精确积分是 23/2 = 8/3 ≈ 2.6667。因此，辛普森法则提供了精确的近似值。

MATLAB 中的辛普森法则

要使用 for 循环在 MATLAB 中实现辛普森法则，您将通过将其划分为子区间并应用辛普森 1/3 法则来计算函数在指定区间上的积分。

我们将使用辛普森 1/3 法则来逼近函数 f(x)=x2 从 a=0 到 b=2 的积分。

% Define the function to integrate
f = @(x) x.^2;

% Define the interval [a, b]
a = 0;
b = 2;

% Number of subintervals (must be even for Simpson's 1/3 Rule)
n = 10;

% Calculate step size
h = (b - a) / n;

% Initialize the sum with the first and last terms
integral = f(a) + f(b);

% Apply Simpson's 1/3 rule using a for loop
for i = 1:n-1
    x = a + i * h;
    
    % Use the appropriate weight for the current point
    if mod(i, 2) == 0
        integral = integral + 2 * f(x);  % even index -> weight is 2
    else
        integral = integral + 4 * f(x);  % odd index -> weight is 4
    end
end

% Multiply by the step size and finalize the result
integral = (h / 3) * integral;

% Display the result
disp(['Approximate integral: ', num2str(integral)]);

在上面的代码中，我们有。

区间 [a,b] 被分成 n 个子区间，因此步长为 h = b-a/n。
求和从端点 a 和 b 处的函数的第一个和最后一个值开始。
循环迭代区间内的内部点。根据索引 i 是奇数还是偶数，相应权重（4 或 2）将应用于该点处的函数值。
循环结束后，结果乘以 h/3 以给出最终的积分近似值。

在 matlab 命令窗口中执行后，我们得到的输出如下。

>> % Define the function to integrate
f = @(x) x.^2;

% Define the interval [a, b]
a = 0;
b = 2;

% Number of subintervals (must be even for Simpson's 1/3 Rule)
n = 10;

% Calculate step size
h = (b - a) / n;

% Initialize the sum with the first and last terms
integral = f(a) + f(b);

% Apply Simpson's 1/3 rule using a for loop
for i = 1:n-1
    x = a + i * h;
    
    % Use the appropriate weight for the current point
    if mod(i, 2) == 0
        integral = integral + 2 * f(x);  % even index -> weight is 2
    else
        integral = integral + 4 * f(x);  % odd index -> weight is 4
    end
end

% Multiply by the step size and finalize the result
integral = (h / 3) * integral;

% Display the result
disp(['Approximate integral: ', num2str(integral)]);

Approximate integral: 2.6667
>>

打印页面