牛顿法代数插值ndash差商表的求法

牛顿法代数插值ndash 差商表的求
法
原文地址：牛顿法代数插值–差商表的求法作者：大关牛顿法代数插值–
差商表的求法
下面的求插商的方法并不是好的求插商的方式，因为他的效率并不是很高，不论是从空间效率还是时间效率，但是下面主要探讨的是一种将塔形的数据转
换成一位数组的方式。

实际上求插商仅通过一个n个元素的一位数组就能解决，但本文强调的是一种思路，希望对大家有所借鉴。

牛顿插商公式：
f[xi,xj]=(f(xj)– f(xi))/(xj– xi)
f[xi,xj,xk]=(f[xj,xk]– f[xi,xj])/(xk– xi)
….
f[x0,x1,x2…,xn]=(f[x1,x2,…,xn]– f[x0,x1,…,xn-1])/(xn– x0)
转换成均插表(或称差商表)形式如下：
定义1：f[xi,xi+1,…xj]简记为f(i,j)其中i=0&&i=n&&j=0&&j=n&&i j；
记f(xi)为f[xi,xi]即f(i,i)
根据定义1可以推出：f[x0,x1]=f(0,1),f[x0,x1…xn]=f(0,n)….
根据定义1：可以将插商表转换为如下形式。

根据上图，可以给出实际一维数组存储时的序列关系，如下图所示：
此时f(0,0)位置是数组下标0，f(1,1)是数组下标为1….
这样，我们从中找出相应的规律。

推论1：已知f(i,j)，n为变量的数目，令k=j– i。

当k不等于0时，f(i,j)在数组中的下标通过计算得：
Index=k*n–((k-1)*k)/2+i
当k等于0时
Index=i。

推论1很容易证明(实际就是一个等差数列求和问题)这里证明略。

推论2：n为变量的数目，则一维数组的长度可以计算得((1+n)*n)/2
推论2可以通过等差数列求和得以证明。

证明略。

推论3：各阶插商就是f(0,k)k=1,2….n.
推论3：根据插商的定义和定义1可以直接推出。

下面先给出一位数组的存储结构代码：
#pragma once
#include vector
#include iostream
#include cassert
//在开始的时候要指定变量的维数
template class T,size_t _Dims class Turriform
{
enum{digits=_Dims}；//记录变量的维数
public：
Turriform(void)；
~Turriform(void)；
//获取f(nleft，nright)处的元素值
T getItem(size_t nLeft,size_t nRight)const；
//设置f(nleft，nright)处的元素值为val void setItem(size_t nLeft,size_t nRight,const T&val)；
private：
//计算f(nleft，nright)的坐标值
size_t index(size_t nLeft,size_t nRight)const；
private：
std：vector Tdatas；//计算时使用的一维数组
}；
template class T,size_t _Dims Turriform T,_Dims：Turriform(void)
{
T demo；
datas.assign(digits*(digits+1)/2,demo)；//初始化所有变量
}
template class T,size_t _Dims size_t Turriform T,_Dims：
index(size_t nLeft,size_t nRight)const
{
if(nLeft nRight)//确保左边的比右边的小
{
size_t tmp=nLeft；
nLeft=nRight；
nRight=nLeft；
}
//计算k值
size_t k=nRight-nLeft；
if(k==0)//k=0的情况
return nLeft；
//k！=0的情况
return k*digits-((k-1)*k)/2+nLeft；
}
template class T,size_t _Dims Turriform T,_Dims：~Turriform(void) {
}
template typename T,size_t _Dims TTurriform T,_Dims：
getItem(size_t nLeft,size_t nRight)const
{
return datas[index(nLeft,nRight)]；
}
template typename T,size_t _Dims void Turriform T,_Dims：
setItem(size_t nLeft,size_t nRight,const T&val)
{
datas[index(nLeft,nRight)]=val；
}
说明：template typename T,size_t _Dims在开始的时候需要指定元素类型和维数大小，维数大小用于确定一维数组实际元素的个数。

#include iostream
#include vector
#include algorithm
#include"Turriform.h"
using std：vector；
using std：cout；
void function(vector double&x,vector double&fx,vector double&ret)
{//放入一个测试的数据
x.reserve(5)；
x.push_back(1)；
x.push_back(2)；
x.push_back(3)；
x.push_back(4)；
x.push_back(5)；
fx.reserve(5)；
fx.push_back(1)；
fx.push_back(4)；
fx.push_back(7)；
fx.push_back(8)；
fx.push_back(6)；
ret.assign(x.size(),0)；
}
template size_t _Dims void Interpolation(const vector
double&x,const vector double&fx,vector double&ret)
{//牛顿插商求解
double val=0；
Turriform double,_Dims tf；
for(size_t i=0；i fx.size()；i++)//先将i-j=0的情况放入插商表{
tf.setItem(i,i,fx[i])；
}
size_t k；
size_t j；
for(size_t n=1；n fx.size()；++n)//计算其他维中的元素{
k=0；
j=k+n；
for(；j fx.size()；++j,++k)
{
val=tf.getItem(k+1,j)-tf.getItem(k,j-1)；//插值公式val/=x[j]-x[k]；
tf.setItem(k,j,val)；
}
}
for(k=0；k x.size()；++k)//获取插值
{
ret[k]=tf.getItem(0,k)；
}
}
int main()
{
vector double x；
vector double fx；
vector double ret；
function(x,fx,ret)；
Interpolation 5(x,fx,ret)；//指定了维
std：copy(ret.begin(),ret.end(),
std：ostream_iterator double(std：cout,""))；//输出结果
std：cout std：endl；
}
上面的代码并不能说明什么，而且比其他的实现方式代码量和复杂度可能更高些，但是，对于问题的分析过程是很重要的。

实际在设计数据结构时，复杂的数据结构往往可以转换成一位数组来进行存储，并可以节省大量的储存空间。

特别是这种类二叉树结构，希望大家遇到问题多思考。

MSN空间完美搬家到新浪博客！。