Nội suy Bessel (Bessel's Interpolation) trong C++

Nội suy (Interpolation) là kỹ thuật ước lượng giá trị chưa biết nằm giữa các giá trị đã biết. Giảm chi phí tính toán khi hàm gốc phức tạp.

Công thức nội suy Bessel

f(u) = {(f(0)+f(1))/2} + {u - 1/2}Δf(0) 
     + {u(u-1)/2!}{(Δ²f(-1) + Δ²f(0))/2}
     + {u(u-1)(u-1/2)/3!}Δ³f(-1) + ...

Trong đó:

• f(0) là gốc (thường là điểm giữa)
• u = (x - f(0)) / h, h là khoảng cách giữa các điểm
• Δ là toán tử sai phân

Ví dụ mã nguồn

#include <iostream>
using namespace std;

float tinhU(float u, int n) {
    if (n == 0) return 1;
    float res = u;
    for (int i = 1; i <= n/2; i++) res *= (u - i);
    for (int i = 1; i < n/2; i++) res *= (u + i);
    return res;
}

int giaiThua(int n) {
    if (n == 1) return 1;
    return n * giaiThua(n-1);
}

int main() {
    int n = 6;
    float x[] = {50, 51, 52, 53, 54, 55};
    float y[6][6];
    y[0][0] = 8.000; y[1][0] = 7.746; y[2][0] = 7.674;
    y[3][0] = 7.571; y[4][0] = 7.469; y[5][0] = 7.231;
    
    // Bảng sai phân
    for (int i = 1; i < n; i++)
        for (int j = 0; j < n-i; j++)
            y[j][i] = y[j+1][i-1] - y[j][i-1];
    
    float xTim = 53.2;
    float tong = (y[2][0] + y[3][0]) / 2;
    int idx = (n % 2) ? n/2 : n/2 - 1;
    float u = (xTim - x[idx]) / (x[1] - x[0]);
    
    for (int i = 1; i < n; i++) {
        if (i % 2)
            tong += ((u - 0.5) * tinhU(u, i-1) * y[idx][i]) / giaiThua(i);
        else
            tong += (tinhU(u, i) * (y[idx][i] + y[--idx][i])) / (giaiThua(i) * 2);
    }
    
    cout << "Giá trị tại " << xTim << " bằng nội suy Bessel = " << tong;
    return 0;
}

Kết quả

Giá trị tại 53.2 bằng nội suy Bessel = 7.54985