給定一些數據，生成函數的方式有兩種：插值，回歸。

插值而得到的函數通過數據點，回歸得到的函數不一定通過數據點。

下面給出拉格朗日插值，牛頓插值和Hermite插值的程序，

具體原理可參考課本，不再贅述。

拉格朗日插值法

線性插值 一次精度 需要2個節點

二次插值 二次精度 需要3個節點

n次插值 n次精度 需要n+1個節點

拉格朗日插值代碼段(根據傳入數據自動判斷次數)：

# 返回多項式

def p(x,a):

"""p(x,a)是x的函數，a是各冪次的系數"""

s = 0

for i in range(len(a)):

s += a[i]*x**i

return s

# n次拉格朗日插值

def lagrange_interpolate(x_list,y_list,x):

"""x_list 待插值的x元素列表y_list 待插值的y元素列表插值以后整個lagrange_interpolate是x的函數"""

if len(x_list) != len(y_list):

raise ValueError("list x and list y is not of equal length!")

# 系數矩陣

A = []

for i in range(len(x_list)):

A.append([])

for j in range(len(x_list)):

A[i].append(pow(x_list[i],j))

b = []

for i in range(len(x_list)):

b.append([y_list[i]])

# 求得各階次的系數

a = lu_solve(A, b) # 用LU分解法解線性方程組，可以使用numpy的類似函數

a = transpose(a)[0] # change col vec a into 1 dimension

val = p(x,a)

print(x,val)

return val

其中lu_solve(A,b)是自己寫的輪子，可以用numpy的numpy.linalg.sovle(A,b)來代替

到后面會有一期講矩陣方程直接法，會有講到如何寫lu_solve()

看一看插值的效果如何

import numpy as np

import matplotlib.pyplot as plt

# 待插值的元素值

x_points = [0,1,2,3,4,5]

y_points = [1,5,4,8,7,12]

# 拉格朗日插值

x = np.linspace(0,5)

y = list(map(lambda t: lagrange_interpolate(x_points,y_points,t),x))

# 畫圖

plt.scatter(x_points,y_points,color = "orange")

plt.plot(x,y)

plt.legend(["lagrange interpolation","scattered points"])

plt.show()拉格朗日插值

可以看到，插值之后的函數可以較好地反映原數據點的情況。

牛頓插值法

涉及到的兩個表，差分表和差商表：差分表差商表

遞歸打印差分表

def difference_list(dlist): # Newton

if len(dlist)>0:

print(dlist)

prev,curr = 0,0

n = []

for i in dlist:

curr = i

n.append(curr - prev)

prev = i

n.pop(0)

difference_list(n)

打印差商表，并以列表的形式返回圖中藍色圈起來的部分

def difference_quotient_list(y_list,x_list = []):

if x_list == []:

x_list = [i for i in range(len(y_list))]

print(y_list)

prev_list = y_list

dq_list = []

dq_list.append(prev_list[0])

for t in range(1,len(y_list)):

prev,curr = 0,0

m = []

k = -1

for i in prev_list:

curr = i

m.append((curr - prev)/(x_list[k+t]-x_list[k]))

prev = i

k+=1

m.pop(0)

prev_list = m

dq_list.append(prev_list[0])

print(m)

return dq_list

牛頓插值代碼段(調用了之前求差商表的代碼)

def newton_interpolate(x_list,y_list,x):

coef = difference_quotient_list(y_list,x_list)

p = coef[0]

for i in range(1,len(coef)):

product = 1

for j in range(i):

product *= (x - x_list[j] )

p += coef[i]*product

return p

if __name__ == '__main__':

import numpy as np

import matplotlib.pyplot as plt

# 待插值的元素值

x_points = [0,1,2,3,4,5]

y_points = [1,5,4,8,7,12]

# 牛頓插值

x = np.linspace(0,5)

y = list(map(lambda t: newton_interpolate(x_points,y_points,t),x))

# 畫圖

plt.scatter(x_points,y_points,color = "orange")

plt.plot(x,y)

plt.legend(["newton interpolation","scattered points"])

plt.show()

插值效果圖牛頓插值

可以看到，相同的插值節點，使用拉格朗日插值和牛頓插值的結果是相同的。

Hermite 插值法

三次Hermite插值不但要求在插值節點上插值多項式與函數值相等，還要求插值多項式的導數與函數導數相等。

進行Hermite插值需要六個信息：

這個比較好理解，通過

兩點的插值有無限種方式。比如：

但是，通過限制住兩個端點的導數值，就可以確定下來大致的插值形狀。(對于Hermite插值，六個條件能確定出唯一的插值結果)

例如，可以編程求出

def hermite(x0,x1,y0,y1,y0_prime,y1_prime,x):

alpha0 = lambda x: ((x-x1)/(x0-x1))**2 * (2*(x-x0)/(x1-x0)+1)

alpha1 = lambda x: ((x-x0)/(x1-x0))**2 * (2*(x-x1)/(x0-x1)+1)

beta0 = lambda x: ((x-x1)/(x0-x1))**2 * (x-x0)

beta1 = lambda x: ((x-x0)/(x1-x0))**2 * (x-x1)

H = alpha0(x)*y0 + alpha1(x)*y1 + beta0(x)*y0_prime + beta1(x)*y1_prime

return H

if __name__ == '__main__':

import numpy as np

import matplotlib.pyplot as plt

f = lambda x: hermite(0,1,0,1,-1,-4,x)

x = np.linspace(0,1)

y = list(map(f,x))

plt.scatter([0,1],[0,1],color = "orange")

plt.plot(x,y)

plt.show()

龍格現象(Runge phenomenon)

然而，并不是所有的函數都可以通過提高插值次數來提高插值準確度。

比如著名的龍格函數

，從[-1,1]取10個點，插值出來的函數是這樣的：

這就是龍格現象，主要原因是誤差項里面的高階導數導致的。什么是龍格現象(Rungephenomenon)？如何避免龍格現象？_MachineLearningwithTuring'sCat-CSDN博客_龍格現象?www.baidu.comhttps://en.wikipedia.org/wiki/Runge%27s_phenomenon?en.wikipedia.org

其實不單單是

，很多偶函數都有這樣的性質：

比如

：

用于生成上圖的代碼：

# 龍格現象的產生

if __name__ == "__main__":

import numpy as np

import matplotlib.pyplot as plt

f = lambda x: 1/(1+25*x**2)

# 待插值的元素值

x_points = np.linspace(-1,1,11)

print(x_points)

y_points = list(map(f,x_points))

# 牛頓插值

x = np.linspace(-1,1)

y = list(map(lambda t: lagrange_interpolate(x_points,y_points,t),x))

# 畫圖

plt.figure("lagrange interpolation")

plt.scatter(x_points,y_points,color = "orange")

plt.plot(x,y)

plt.legend(["lagrange interpolation","scattered points"])

plt.show()

如何避免龍格現象，就需要用到分段插值了。

下一篇將講解分段插值。