牛顿插值法算法

牛顿插值法(算法)

Newton Interpolation

牛顿插值算法是根据n + 1个点x0, x1, ... x n(x0 < x1 < ... x n)与函数値f(x0), f(x1) , ... , f(x n)求出n次多項式p(x)。再通过次多項式p(x)求出任意的点x所对应的f(x)的算法。

１．求n阶差分商f[x0, x1, ... 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;

}