SOFTIST 编程方法筆記 →目録

牛顿插值法(算法)

Newton Interpolation

牛顿插值算法是根据 n + 1个点x₀, x₁, ... x_n(x₀ < x₁ < ... x_n)与函数値f(x₀), f(x₁) , ... , f(x_n)求出n次多項式p(x)。再通过次多項式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;
}

３．测试。用正弦波的20个采样点，还原出正弦波曲线。计算速度很慢，需要改进程序。

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);
    }
}

４．n阶差分商的计算方法的改善。上面的递归算法虽然很好读，但速度太慢不能实用。n阶差分商与、x和f(x)没有关系，所以可以事先把它一次算好，以提高整体処理速度。因为是事先计算，所以用递归算法，或迭代算法都无所谓。不过为了记录算法，用迭代算法计算。这回速度提高了。

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