日志
算法描述
1 牛顿插值法基本思路
给定插值点序列(xi, f(xi)), i=0,1,…,n。构造牛顿插值多项式Nn(u)。输入要计算的函数点x并计算Nn(x)的值，利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。
2 牛顿插值法计算步骤
1) 输入n的值及 (xi, f(xi)), i=0,1,…,n；要计算的函数点 x。
2) 对给定的x,由
计算Nn(x)的值。
3) 输出Nn(x)。
===== 上面其实是直接抄来的 =====
===== 下面是我的实现 =====
// 依赖于：矩阵和向量类
double newtonInterpolate( NA::VecDP const& vx,
NA::VecDP const& vy,
double x )
{
int const n=vx.size();
NA_ASSERT ( n==vy.size() );
NA::MatDP table(n,n); // 差商表
int i,j;
// 计算差商 {{{
for ( i=0; i<n; ++i ) {
table[i][0]=vy[i];
}
for ( i=1; i<n; ++i ) {
for ( j=1; j<=i; ++j ) {
table[i][j]= ( table[i][j-1]-table[i-1][j-1] ) / ( vx[i]-vx[i-j]);
// 可以看出 table[i-1][j-1] 从此不需要再用到，
// 因此 table 用一维数组表示其实就可以了
}
}
// }}}
double t=1.0;
double ans=table[0][0];
for ( i=1; i<n; ++i ) {
t *= (x-vx[i-1]);
ans += t*table[i][i];
}
return ans;
}
// 空间优化的版本
// 按列计算，从左向右、从下往上
double newtonInterpolate_O1
( NA::VecDP const& vx,
NA::VecDP const& vy,
double x )
{
int const n=vx.size();
NA_ASSERT(n==vy.size());
NA::VecDP table(n);
int i,j;
for ( i=0; i<n; ++i )
table[i]=vy[i];
for ( i=1; i<=n; ++i )
for ( j=n-1; j>=i; --j )
table[j] = (table[j]-table[j-1]) / (vx[j]-vx[j-i]);
double t=1.0;
double ans=table[0];
for ( i=1; i<n; ++i ) {
t *= (x-vx[i-1]);
ans += t*table[i];
}
return ans;
}
OVER
转发至微博
转发至微博
0人 |
分享到：
评论