拉格朗日插值法的最大毛病就是每次引入一個新的插值節點，基函數都要發生變化，這在一些實際生產環境中是不合適的，有時候會不斷的有新的測量數據加入插值節點集，
因此，通過尋找n個插值節點構造的的插值函數與n+1個插值節點構造的插值函數之間的關系，形成了牛頓插值法。推演牛頓插值法的方式是歸納法，也就是計算Ln（x）- Ln+1（x），並且從n=1開始不斷的迭代來計算n+1時的插值函數。
牛頓插值法的公式是：
注意：在程序中我用W 代替 
計算牛頓插值函數關鍵是要計算差商，n階差商的表示方式如下：
關於差商我在這里並不討論
計算n階差商的公式是這樣：
很明顯計算n階差商需要利用到兩個n-1階差商，這樣在編程的時候很容易想到利用遞歸來實現計算n階差商，不過需要注意的是遞歸有棧溢出的潛在危險，在計算差商的時候
更是如此，每一層遞歸都會包含兩個遞歸，遞歸的總次數呈滿二叉樹分布：
這意味着遞歸次數會急劇增加:(。所以在具體的應用中還需要根據應用來改變思路或者優化代碼。
廢話少說，放碼過來。
首先寫最關鍵的一步，也就是計算n階差商：
"""
@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
OK， 牛頓插值法最重要的兩部分都有了，下面就是將這兩部分組合成牛頓插值函數，如果是c之類的語言就需要保存一些中間數據了，我利用了Python的閉包直接返回一個牛頓插值函數，閉包可以利用到它所處的函數之中的上下文數據。
"""
@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，這個從函數的表達式可以找到原因。
構造一些插值節點
''' 插值節點, 這里用二次函數生成插值節點，每兩個節點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] # 根據插值函數獲得測試用例的縱坐標
完整代碼：
# -*- 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()