插值法的偉大作用我就不說了。。。。
那么貼代碼？
首先說一下下面幾點：
1. 已有的數據樣本被稱之為 “插值節點”
2. 對於特定插值節點，它所對應的插值函數是必定存在且唯一的（關於這個的證明我暫時不說了，如果哪天我回頭看看我的blog有點寒磣，我再再補上）
也就是說對於同樣的插值樣本來說，用不同方法求得的插值函數本質上其實是一樣的。
3. 拉格朗日插值法依賴於每個插值節點對應的插值基函數，也就是說每個插值節點都有對應的插值基函數。
4. 拉格朗日插值函數最終由所有插值節點中每個插值節點的縱坐標值與它對應的插值函數的積的和構成。
5. 也就是說拉格朗日插值法關鍵在於求基函數
6. 拉格朗日插值法並不好，當每一次加入新的插值節點的時候，所有的系數都要重算一遍
計算插值函數參數的最本質的方法是解下面的矩陣方程：
我先構造一組插值樣本：
sr_x = [i for i in range(-50, 50, 10)] sr_fx = [i**2 for i in sr_x]
sr_x 為樣本的橫坐標，sr_fx為樣本的縱坐標，樣本的定義域是[-50, 50]，樣本之間距離為10
縱坐標為橫坐標的平方，也就是說這是一個二次函數 sr_fx = sr_x ^2
這樣一組樣本的圖像是這樣的：
然后我用這一組插值節點來獲得每個節點對應的插值基函數：
基函數的計算公式是 li = (x - x0)(x - x1) ... (x - x(i-1))(x - x(i+1)) ... (x - xn)/(xi - x0)(xi - x1) ... (xi - x(i-1))(xi - x(i+1)) ... (xi - xn)
其中W = (x - x0)(x - x1) ... (x - x(i-1))(x - x(i+1)) ... (x - xn)
c = (xi - x0)(xi - x1) ... (xi - x(i-1))(xi - x(i+1)) ... (xi - xn)
因此 li = W / c
注意：在計算i的插值基函數的時候，公式c與W里面是沒有被減數為Xi項的。（要是有的話W與c就為0了）
這一段代碼是這樣的：
def get_li(xi, x_set = []):
def li(Lx):
W = 1; c = 1
for each_x in x_set:
if each_x == xi:
continue
W = W * (Lx - each_x)
for each_x in x_set:
if each_x == xi:
continue
c = c * (xi - each_x)
# 這里一定要轉成float類型，否則極易出現嚴重錯誤. 原因就不說了
return W / float(c)
return li
這段代碼用到了閉包，這樣可以返回一個基函數，並且這個函數以get_li的內部為上下文（上下文這個翻譯總讓人感覺怪怪的，似乎與寫作文有某種聯系）。
當獲得基函數之后就是累加基函數與插值節點縱坐標的乘積構成拉格朗日插值函數了
這一段的代碼是這樣：
"""
@brief: 獲得拉格朗日插值函數
@param: x 插值節點的橫坐標集合
@param: fx 插值節點的縱坐標集合
@return: 參數所指定的插值節點集合對應的插值函數
"""
def get_Lxfunc(x = [], fx = []):
set_of_lifunc = []
for each in x: # 獲得每個插值點的基函數
lifunc = get_li(each, x)
set_of_lifunc.append(lifunc) # 將集合x中的每個元素對應的插值基函數保存
def Lxfunc(Lx):
result = 0
for index in range(len(x)):
result = result + fx[index]*set_of_lifunc[index](Lx) #根據根據拉格朗日插值法計算Lx的值
print fx[index]
return result
return Lxfunc
到這里就大功告成了，我用上面給的插值節點計算出的插值函數是這樣的：
完整代碼如下：
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 16 13:26:58 2016
@author: tete
@brief: 拉格朗日插值法
"""
import matplotlib.pyplot as plt
"""
@brief: 獲得拉格朗日插值基函數
@param: xi xi為第i個插值節點的橫坐標
@param: x_set 整個插值節點集合
@return: 返回值為參數xi對應的插值基函數
"""
def get_li(xi, x_set = []):
def li(Lx):
W = 1; c = 1
for each_x in x_set:
if each_x == xi:
continue
W = W * (Lx - each_x)
for each_x in x_set:
if each_x == xi:
continue
c = c * (xi - each_x)
# 這里一定要轉成float類型，否則極易出現嚴重錯誤. 原因就不說了
return W / float(c)
return li
"""
@brief: 獲得拉格朗日插值函數
@param: x 插值節點的橫坐標集合
@param: fx 插值節點的縱坐標集合
@return: 參數所指定的插值節點集合對應的插值函數
"""
def get_Lxfunc(x = [], fx = []):
set_of_lifunc = []
for each in x: # 獲得每個插值點的基函數
lifunc = get_li(each, x)
set_of_lifunc.append(lifunc) # 將集合x中的每個元素對應的插值基函數保存
def Lxfunc(Lx):
result = 0
for index in range(len(x)):
result = result + fx[index]*set_of_lifunc[index](Lx) #根據根據拉格朗日插值法計算Lx的值
print fx[index]
return result
return Lxfunc
"""
demo:
"""
if __name__ == '__main__':
''' 插值節點, 這里用二次函數生成插值節點，每兩個節點x軸距離位10 '''
sr_x = [i for i in range(-50, 50, 10)]
sr_fx = [i**2 for i in sr_x]
Lx = get_Lxfunc(sr_x, sr_fx) # 獲得插值函數
tmp_x = [i for i in range(-45, 45)] # 測試用例
tmp_y = [Lx(i) for i in tmp_x] # 根據插值函數獲得測試用例的縱坐標
''' 畫圖 '''
plt.figure("play")
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()