Python数值计算（11）

本篇介绍一下多项式插值中，拉格朗日法的原理及其实现。

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完全相同。