python插值程序_计算方法（2）——插值法（附Python程序）

给定一些数据，生成函数的方式有两种：插值，回归。

插值而得到的函数通过数据点，回归得到的函数不一定通过数据点。

下面给出拉格朗日插值，牛顿插值和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_phenomenonen.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)