本篇介绍一下多项式插值中,拉格朗日法的原理及其实现。
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类型。