ニュートン補間(アルゴリズム)

Newton Interpolation

ニュートン補間とは、n + 1個の点x₀, x₁, ... x_n(x₀ < x₁ < ... x_n)と、関数値f(x₀), f(x₁) , ... , f(x_n)に基づいて、n次多項式p(x)を、下記の式で求めて、任意のxからf(x)を算出する補間手法です。

１．n階差分商というf[x₀, x₁, ... x_n]のプログラム。再帰呼び出し。

#define N 20
typedef struct TagXYVALUE
{
    double x;
    double y;
} XYVALUE;

XYVALUE val[N+1];

//階差分商(Divided Differences)
double f(int n0, int ni)
{
    if (n0 == ni)
        return val[n0].y;
    if (n0 + 1== ni)
        return (val[ni].y - val[n0].y) / (val[ni].x - val[n0].x);
    else
        return (f(n0+1, ni) - f(n0, ni-1)) / (val[ni].x - val[n0].x); 
}

２．ニュートン補間のメイン関数で、p(x_n)を求めるプログラム。

double NewtonInterpolation(double x)
{
    double t = 1.0;
    double ft;
    double p = val[0].y; //P(0) = f[0]

    for(int i = 1; i <= N; i++)
    {
        t = t * (x - val[i-1].x);
        ft = f(0, i) * t;
        p = p + ft;
    }
    return p;
}

３．テスト。遅いですので、プログラムの改善必要があると思います。

void CNewtonInterpolationTestView::OnDraw(CDC* pDC)
{
    for (int i = 0; i <= N; i ++)
    {
        val[i].x = i * 15 * atan(1.0) / 45.0 * 2;
        val[i].y = sin(val[i].x);
        pDC->Rectangle((int)(val[i].x * 20) - 2, 150 - (int)(val[i].y * 50) - 2,
                                 (int)(val[i].x * 20) + 2, 150 - (int)(val[i].y * 50) + 2);
    }
    for (int j = 0; j <= N*15; j += 5)
    {
        double x = j * atan(1.0) / 45.0 * 2;
        double y = NewtonInterpolation(x);
        pDC->SetPixel((int)(x * 20) - 2, 150 - (int)(y * 50) - 2, 0x000000ff);
    }
}

４．階差分商の計算方法の改善。上の再帰呼び出しの方法が分かりやすいですが、遅くて実用できません。階差分商は、xやf(x)に無関係なので、事前に算出しておけば、処理速度を上げられます。事前に計算するので、再帰呼び出し方法でも、非再帰（ループ）の方法でもOKですが、下記の非再帰の方法もメモします。今度は、早いです。

double dd[N+1];
void CallDividedDifferences()
{
    for (int n = 0; n <= N; n ++)
    dd[n] = val[n].y;

    for (int i = 1; i <= N; i ++)
    {
        for (int j = N; j >= i; j --)
            dd[j] = (dd[j] - dd[j - 1]) / (val[j].x - val[j-i].x);
    }
}
double NewtonInterpolation_DD(double x)
{
    double t = 1.0;
    double ft;
    double p = val[0].y; //P(0) = f[0]

    for(int i = 1; i <= N; i++)
    {
        t = t * (x - val[i-1].x);
        ft = dd[i] * t;
        p = p + ft;
    }
    return p;
}