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

下面的求插商的方法并不是好的求插商的方式，因为他的效率并不是很高，不论是从空间效率还是时间效率，但是下面主要探讨的是一种将塔形的数据转换成一位数组的方式。实际上求插商仅通过一个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

#include

#include

//在开始的时候要指定变量的维数

template

_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 datas;//计算时使用的一维数组

};

template

_Dims>

Turriform::Turriform(void)

{

T demo;

datas.assign(digits*(digits + 1)/2,demo);//初始化所有变量

}

template

_Dims>

size_t

Turriform::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

_Dims>

Turriform::~Turriform(void)

{

}

template

_Dims>

T

Turriform::getItem(size_t

nLeft, size_t nRight)const

{

return

datas[index(nLeft,nRight)];

}

template

_Dims>

void Turriform::setItem(size_t

nLeft,size_t nRight,const T

& val)

{

datas[index(nLeft,nRight)] = val;

}

说明：template

_Dims>在开始的时候需要指定元素类型和维数大小，维数大小用于确定一维数组实际元素的个数。

#include

#include

#include

#include "Turriform.h"

using std::vector;

using std::cout;

void function(vector & x,

vector & fx,

vector &

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

_Dims>

void Interpolation(const vector & x,

const

vector &

fx,vector &

ret)

{//牛顿插商求解

double val = 0;

Turriform 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 x;

vector fx;

vector ret;

function(x, fx, ret);

Interpolation<5>(x, fx,

ret);//指定了维

std::copy(ret.begin(),ret.end(),

std::ostream_iterator(std::cout," "));//输出结果

std::cout<<:endl>

}

上面的代码并不能说明什么，而且比其他的实现方式代码量和复杂度可能更高些，但是，对于问题的分析过程是很重要的。实际在设计数据结构时，复杂的数据结构往往可以转换成一位数组来进行存储，并可以节省大量的储存空间。特别是这种类二叉树结构，希望大家遇到问题多思考。