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

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

转换成均插表(或称差商表)形式如下:

 

 http://s9/middle/616e189f49f825966ed78&690代数插值 – 差商表的求法" TITLE="牛顿法 代数插值 – 差商表的求法" />

定义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：可以将插商表转换为如下形式。

http://s16/middle/616e189f49f825dddc8af&690代数插值 – 差商表的求法" TITLE="牛顿法 代数插值 – 差商表的求法" />

 

根据上图，可以给出实际一维数组存储时的序列关系，如下图所示：

http://s2/middle/616e189f07659d6a19921&690代数插值 – 差商表的求法" TITLE="牛顿法 代数插值 – 差商表的求法" />

 

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

};

 

 

 

 

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>

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

    

}

 

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