牛顿插值法——用Python进行数值计算

拉格朗日插值法的最大毛病就是每次引入一个新的插值节点，基函数都要发生变化，这在一些实际生产环境中是不合适的，有时候会不断的有新的测量数据加入插值节点集，

因此，通过寻找n个插值节点构造的的插值函数与n+1个插值节点构造的插值函数之间的关系，形成了牛顿插值法。推演牛顿插值法的方式是归纳法，也就是计算Ln（x）- Ln+1（x），并且从n=1开始不断的迭代来计算n+1时的插值函数。

 

　　牛顿插值法的公式是：

　　注意：在程序中我用W 代替 

　　计算牛顿插值函数关键是要计算差商，n阶差商的表示方式如下：

　　　　　　　　　　　　　　　　　　　　　　　　

 

　   关于差商我在这里并不讨论

　　计算n阶差商的公式是这样：

　　很明显计算n阶差商需要利用到两个n-1阶差商，这样在编程的时候很容易想到利用递归来实现计算n阶差商，不过需要注意的是递归有栈溢出的潜在危险，在计算差商的时候

更是如此，每一层递归都会包含两个递归，递归的总次数呈满二叉树分布：

　　　　

　　这意味着递归次数会急剧增加:(。所以在具体的应用中还需要根据应用来改变思路或者优化代码。

 

　　废话少说，放码过来。

　　首先写最关键的一步，也就是计算n阶差商：

1 2 3 4 5 6 7 8 9 10 11 12 13

""" @brief:   计算n阶差商 f[x0, x1, x2 ... xn] @param:   xi   所有插值节点的横坐标集合                                                        o @param:   fi   所有插值节点的纵坐标集合                                                      /   \ @return:  返回xi的i阶差商(i为xi长度减1)                                                     o     o @notice:  a. 必须确保xi与fi长度相等                                                        / \   / \            b. 由于用到了递归，所以留意不要爆栈了.                                           o   o o   o            c. 递归减递归(每层递归包含两个递归函数), 每层递归次数呈二次幂增长，总次数是一个满二叉树的所有节点数量(所以极易栈溢出)                                                                                     """ def  get_order_diff_quot(xi  =  [], fi  =  []):      if  len (xi) >  2  and  len (fi) >  2 :          return  (get_order_diff_quot(xi[: len (xi)  -  1 ], fi[: len (fi)  -  1 ])  -  get_order_diff_quot(xi[ 1 : len (xi)], fi[ 1 : len (fi)]))  /  float (xi[ 0 ]  -  xi[ - 1 ])      return  (fi[ 0 ]  -  fi[ 1 ])  /  float (xi[ 0 ]  -  xi[ 1 ])

　　看上面的牛顿插值函数公式，有了差商，还差

　　这个就比较好实现了：

1 2 3 4 5 6 7 8 9 10 11 12 13 14

""" @brief:  获得Wi(x)函数;           Wi的含义举例 W1 = (x - x0); W2 = (x - x0)(x - x1); W3 = (x - x0)(x - x1)(x - x2) @param:  i  i阶(i次多项式) @param:  xi  所有插值节点的横坐标集合 @return: 返回Wi(x)函数 """ def  get_Wi(i  =  0 , xi  =  []):      def  Wi(x):          result  =  1.0          for  each  in  range (i):              result  * =  (x  -  xi[each])          return  result      return  Wi

　　　　

　　　　OK， 牛顿插值法最重要的两部分都有了，下面就是将这两部分组合成牛顿插值函数，如果是c之类的语言就需要保存一些中间数据了，我利用了Python的闭包直接返回一个牛顿插值函数，闭包可以利用到它所处的函数之中的上下文数据。

1 2 3 4 5 6 7 8 9 10 11

""" @brief: 获得牛顿插值函数 @ """ def  get_Newton_inter(xi  =  [], fi  =  []):      def  Newton_inter(x):          result  =  fi[ 0 ]          for  i  in  range ( 2 ,  len (xi)):              result  + =  (get_order_diff_quot(xi[:i], fi[:i])  *  get_Wi(i - 1 , xi)(x))          return  result      return  Newton_inter

　　　　上面这段代码就是对牛顿插值函数公式的翻译，注意get_Wi函数的参数是i-1，这个从函数的表达式可以找到原因。

 

　　构造一些插值节点

　

1 2 3

''' 插值节点, 这里用二次函数生成插值节点，每两个节点x轴距离位10 '''     sr_x  =  [i  for  i  in  range ( - 50 ,  51 ,  10 )]     sr_fx  =  [i * * 2  for  i  in  sr_x]

　　　

 

　获得牛顿插值函数

1 2 3 4

Nx  =  get_Newton_inter(sr_x, sr_fx)             # 获得插值函数   tmp_x  =  [i  for  i  in  range ( - 50 ,  51 )]           # 测试用例 tmp_y  =  [Nx(i)  for  i  in  tmp_x]                # 根据插值函数获得测试用例的纵坐标

　　

 

完整代码：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

# -*- coding: utf-8 -*- """ Created on Thu Nov 17 18:34:21 2016   @author: tete @brief: 牛顿插值法 """     import  matplotlib.pyplot as plt   """ @brief:   计算n阶差商 f[x0, x1, x2 ... xn] @param:   xi   所有插值节点的横坐标集合                                                        o @param:   fi   所有插值节点的纵坐标集合                                                      /   \ @return:  返回xi的i阶差商(i为xi长度减1)                                                     o     o @notice:  a. 必须确保xi与fi长度相等                                                        / \   / \            b. 由于用到了递归，所以留意不要爆栈了.                                           o   o o   o            c. 递归减递归(每层递归包含两个递归函数), 每层递归次数呈二次幂增长，总次数是一个满二叉树的所有节点数量(所以极易栈溢出)                                                                                     """ def  get_order_diff_quot(xi  =  [], fi  =  []):      if  len (xi) >  2  and  len (fi) >  2 :          return  (get_order_diff_quot(xi[: len (xi)  -  1 ], fi[: len (fi)  -  1 ])  -  get_order_diff_quot(xi[ 1 : len (xi)], fi[ 1 : len (fi)]))  /  float (xi[ 0 ]  -  xi[ - 1 ])      return  (fi[ 0 ]  -  fi[ 1 ])  /  float (xi[ 0 ]  -  xi[ 1 ])             """ @brief:  获得Wi(x)函数;           Wi的含义举例 W1 = (x - x0); W2 = (x - x0)(x - x1); W3 = (x - x0)(x - x1)(x - x2) @param:  i  i阶(i次多项式) @param:  xi  所有插值节点的横坐标集合 @return: 返回Wi(x)函数 """ def  get_Wi(i  =  0 , xi  =  []):      def  Wi(x):          result  =  1.0          for  each  in  range (i):              result  * =  (x  -  xi[each])          return  result      return  Wi                         """ @brief: 获得牛顿插值函数 @ """ def  get_Newton_inter(xi  =  [], fi  =  []):      def  Newton_inter(x):          result  =  fi[ 0 ]          for  i  in  range ( 2 ,  len (xi)):              result  + =  (get_order_diff_quot(xi[:i], fi[:i])  *  get_Wi(i - 1 , xi)(x))          return  result      return  Newton_inter                   """ demo: """ if  __name__  = =  '__main__' :          ''' 插值节点, 这里用二次函数生成插值节点，每两个节点x轴距离位10 '''      sr_x  =  [i  for  i  in  range ( - 50 ,  51 ,  10 )]      sr_fx  =  [i * * 2  for  i  in  sr_x]             Nx  =  get_Newton_inter(sr_x, sr_fx)             # 获得插值函数            tmp_x  =  [i  for  i  in  range ( - 50 ,  51 )]           # 测试用例      tmp_y  =  [Nx(i)  for  i  in  tmp_x]                # 根据插值函数获得测试用例的纵坐标               ''' 画图 '''      plt.figure( "I love china" )      ax1  =  plt.subplot( 111 )      plt.sca(ax1)      plt.plot(sr_x, sr_fx, linestyle  =  ' ', marker=' o ', color=' b')      plt.plot(tmp_x, tmp_y, linestyle  =  '--' , color = 'r' )      plt.show()