本篇介绍一下多项式插值中，拉格朗日法的原理及其实现。
1. 一点数学知识
先引用数学背景。如果给定N个点，然后要求一个多项式通过这N个点，最简单直接的方式是列出线性方程求解，N个点可以确定N个未知量，则所求的拟合多项式，其最高次幂就是（N-1）。
在很多教材上，关于该插值多项式的算法如下：
这个公式具有很强烈的轮换性，和线性方程组的克莱姆法则很相似（其实本质上就是一样的），所以实现起来也是非常简单。
2. 算法实现
def lagrangeIntp(x:ndarray,y:ndarray):
'''
返回基于x,y点对的拉格朗日插值多项式
x,y必须有相同的长度
'''
n=len(x)
ret=P([0])
for k in range(n):
roots=list(x)
s=roots[k]
del roots[k]
p=P.fromroots(roots)
p=p*y[k]/p(s)
ret+=p
return ret当然，上述函数并未对输入进行检测，例如x，y的尺寸应该是一致的等。
然后测试一下该拟合函数，仍旧以之前的函数为例：
a=P([1,-2,0,3])
print(a) # 1.0 - 2.0·x + 0.0·x² + 3.0·x³
x=np.arange(-2,2)
print(type(x))
y0=a(x)
b=lagrangeIntp(x,y0)
print(b) # 1.0 - 2.0·x + 0.0·x² + 3.0·x³
c=P.fit(x,y0,deg=3,domain=[-1,1])
print(c) # 1.0 - 2.0·x - (3.87289473e-15)·x² + 3.0·x³
X=np.linspace(-2,1,100)
plt.plot(X,a(X),'r')
plt.plot(X,b(X),'g--')
plt.plot(X,c(X),'b--')
plt.grid()
plt.show()上面的例子中，先用原多项式a生成了四个点(-2,-19)，(-1,0)，(0,1)，(1,2)，然后将其作为拟合的输入，计算得到多项式b，并加入了通过fit拟合出来的多项式c，最后绘制三者的图形，如下：
可以看到三个图形是几乎完全重叠的，特别是使用拉格朗日法的插值多项式b，与原多项式a完全相同。
3. 现成工具
目前在numpy中暂未找到拉格朗日插值的模块，但是在Scipy中有对应的模块，那就是scipy.interpolate，用起来很简单，只要几行代码就可以搞定：
from scipy.interpolate import lagrange
import numpy as np
from numpy.polynomial import Polynomial as P
np.polynomial.set_default_printstyle('unicode')
x=np.array([-2,-1,0,1,2])
y=np.array([-19,0,1,2,21])
g=P(lagrange(x,y).coef[::-1])
print(g)
# 1.0 - 2.0·x¹ + 1.1102230246251565e-16·x² + 3.0·x³
上面这个代码中，x,y是根据前面的函数生成的，注意，这里一共有五个点，所以最高次数有可能会达到4次，但是算法运行的结果基本上和前面的函数相同，只有x的2次幂有一个非常小的项，可以认为是浮点运算误差造成的。另外就是scipy.interpolate.lagrange返回的是一个numpy.poly1d的对象，这个在numpy已经是一个过时的类型了，语句P(lagrange(x,y).coef[::-1])是将其转化为了常用的Polynomial类型。
