C++で実装するベッセル補間：公式の解説とサンプルコード

補間（Interpolation）とは、既知の値の間にある未知の値を推定する手法の一つです。離散的な既知データ点の範囲内に新しいデータ点を構築するプロセスを指します。

補間を使用する主な理由の一つは、計算コストの削減です。ある値を求めるための数式（関数）が複雑すぎたり、計算に多大な時間やリソースが必要だったりする場合、補間が有効な選択肢となります。元の関数を使って少数のデータ点だけを計算し、残りの値は補間によって推定します。完全に正確とは言えませんが、十分に近い値が得られます。

つまり、ここでは計算コストの削減と手軽さが、補間誤差による損失を上回るというトレードオフが成立しているのです。

ベッセルの補間公式

f(u) = {(f(0)+f(1))/2} + {u - 1/2}𝛥f(0) + {u(u-1)/2!}{(𝛥2 f(-1) + 𝛥2 f(0))/2} + {u(u-1)(u - 1/2)/3!}𝛥3f(-1) + {u(u+1)(u-1)(u-2)/4!}{(𝛥4f(-2) + 𝛥4f(-1))/2} + ...

ここで、

• f(0) は原点となる点で、通常は中央の点を指します。
• u = x - f(0) / h（h は差分の間隔）

C++による実装例

以下は、ベッセル補間を実際に動かすためのC++プログラムです。階差表を作成し、公式に沿って順次計算を行います。

#include <iostream>
using namespace std;

// u の多項式部分を計算する関数
calU(float u, int n){
    if (n == 0)
        return 1;
    float result = u;
    for (int i = 1; i <= n / 2; i++)
        result = result*(u - i);
    for (int i = 1; i < n / 2; i++)
        result = result*(u + i);
    return result;
}

// 階乗を計算する関数
int factorial(int n){
    if(n == 1)
        return 1;
    return n * factorial(n-1);
}

int main(){
    int n = 6;
    float x[] = { 50, 51, 52, 53, 54, 55 };
    float y[n][n];
    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;

    // 階差表の作成
    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 value = 53.2;
    float sum = (y[2][0] + y[3][0]) / 2;
    int index;
    if (n % 2)
        index = n/2;
    else
        index = n/2 - 1;

    float u = (value - x[index]) / (x[1] - x[0]);

    // ベッセルの補間公式に基づく合計の計算
    for (int i = 1; i < n; i++) {
        if (i % 2)
            sum += (((u-(0.5))*calU(u, i - 1)*y[index][i])/factorial(i));
        else
            sum += ((calU(u, i)*(y[index][i]+y[-- index][i])/(factorial(i)*2)));
    }

    cout<<"Value at "<<value<<" found using Bessels's interpolation is "<<sum;
    return 0;
}

実行結果

Value at 53.2 found using Bessels's interpolation is 7.54985

このプログラムでは、x = 50〜55 の6つの既知データ点から、x = 53.2 における関数値をベッセル補間によって推定しています。結果として 7.54985 が得られました。このように、少数の既知点から中間の値を効率よく近似できるのが補間法の強みです。