Quy tắc Simpson 1/3 cho tích phân xác định

Giới thiệu

Giống như Quy tắc Thang (Trapezoidal Rule), quy tắc Simpson 1/3 cũng được dùng để tìm giá trị tích phân từ khoảng a đến b. Sự khác biệt chính giữa quy tắc thang và quy tắc Simpson 1/3 là: ở quy tắc thang, toàn bộ đoạn được chia thành các hình thang, nhưng trong trường hợp này, mỗi hình thang cũng được chia thành hai phần.

Đối với quy tắc này, ta sẽ tuân theo công thức sau:

Ở đây h là độ rộng của khoảng chia, và n là số khoảng chia. Ta có thể tìm h bằng cách sử dụng:

Đầu vào và Đầu ra

Đầu vào: Hàm f(x): (x+(1/x)). Giới hạn dưới và trên: 1, 2. Số khoảng chia: 20.
Đầu ra: Kết quả là: 2.19315

Thuật toán

integrateSimpson(a, b, n)
Đầu vào: Giới hạn dưới và trên của tích phân và số khoảng chia n.
Đầu ra: Kết quả sau tích phân.

Begin
    h := (b - a)/n
    res := f(a) + f(b)
    lim := n/2
    for i := 1 to lim, do
        oddSum := oddSum + f(a + (2i - 1)h)
    done
    oddSum := oddSum * 4
    for i := 1 to lim-1, do
        evenSum := evenSum + f(a + 2ih)
    done
    evenSum := evenSum * 2
    res := res + oddSum + evenSum
    res := res * (h/3)
    return res
End

Ví dụ (C++)

#include <iostream>
#include <cmath>
using namespace std;

float mathFunc(float x) {
    return (x+(1/x));  // hàm 1 + 1/x
}

float integrate(float a, float b, int n) {
    float h, res = 0.0, oddSum = 0.0, evenSum = 0.0, lim;
    int i;
    h = (b-a)/n;  // tính khoảng cách giữa hai khoảng chia
    res = (mathFunc(a)+mathFunc(b));  // tổng ban đầu dùng f(a) và f(b)
    lim = n/2;
    for(i = 1; i <= lim; i++)
        oddSum += mathFunc(a+(2*i-1)*h);  // tổng các số ở vị trí lẻ
    oddSum *= 4;  // tổng lẻ nhân 4
    for(i = 1; i < lim; i++)
        evenSum += mathFunc(a+(2*i)*h);  // tổng các số ở vị trí chẵn
    evenSum *= 2;  // tổng chẵn nhân 2
    res += oddSum+evenSum;
    res *= (h/3);
    return res;  // kết quả tích phân
}

main() {
    float result, lowLim, upLim;
    int interval;
    cout << "Enter Lower Limit, Upper Limit and interval: ";
    cin >> lowLim >> upLim >> interval;
    result = integrate(lowLim, upLim, interval);
    cout << "The answer is: " << result;
}

Kết quả

Enter Lower Limit, Upper Limit and interval: 1 2 20
The answer is: 2.19315