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
>>

列印頁面