本篇说说Neville方法。Neville方法的基础是，插值多项式可以递归的生成，有时进行插值的目的是为了计算某个点的值，这个时候并不需要将拟合曲线完全求出，而是可以通过递归的方式进行计算，具体操作如下：
例如有多个点，为了表达方便，这里取五个点，与之对应的值为
，则可以构造列表：
其中
等等，随后计算方式为：
最后递归地计算到，该值就是在对应值为x时的拟合值。
实现代码为：
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial as P
np.polynomial.set_default_printstyle("unicode")
def neville(x:np.ndarray,y:np.ndarray,p:float):
n=len(x)
Q=np.array([0.0]*n*n).reshape(n,n)
Q[:,0]=y.copy()
for c in range(1,n):
for r in range(c,n):
Q[r,c]=((p-x[r-c])*Q[r,c-1]-(p-x[r])*Q[r-1,c-1])/(x[r]-x[r-c])
return Q注意返回的是一个ndarray对象。
如何验证其拟合值呢？
第一组，用原来的函数，数据为：
| X | Y |
| -2 | -19 |
| -1 | 0 |
| 0 | 1 |
| 1 | 2 |
| 2 | 21 |
现在要计算在0.7处的拟合值。测试代码如下：
x=[-2,-1,0,1,2]
y=[-19,0,1,2,21]
print(x)
print(y)
ret=neville(x,y,0.7)
print(ret)运行结果为：
[[-19. 0. 0. 0. 0. ]
[ 0. 32.3 0. 0. 0. ]
[ 1. 1.7 -9.01 0. 0. ]
[ 2. 1.7 1.7 0.629 0. ]
[ 21. -3.7 -0.19 0.629 0.629]]函数值应该为：
print(ret[4,4]) # 0.6289999999999998
a=P([1,-2,0,3])
print(a(0.7)) # 0.6289999999999998可见结果堪称完美。
测试另外一组数据：
| X | Y |
| 0 | 1 |
| 0.25 | 1.58033897 |
| 0.5 | 2.19829503 |
| 0.75 | 2.60553848 |
| 1 | 2.71828183 |
需要估算其在0.6处的值：函数调用如下：
x=np.array([0,0.25,0.5,0.75,1])
y=np.array([1,1.58033897,2.19829503,2.19829503,2.71828183] )
a=neville(x,y,0.6)
print(a)
print(a[4,4]) # 2.2488997156640003运行结果是：
[[1. 0. 0. 0. 0. ]
[1.58033897 2.39281353 0. 0. 0. ]
[2.19829503 2.44547745 2.45601024 0. 0. ]
[2.19829503 2.19829503 2.27244976 2.30916185 0. ]
[2.71828183 1.88630295 2.13589661 2.20872496 2.24889972]]因此，可知在0.6处的拟合值，其结果为2.2488997156640003。
然后再看看使用fit函数的结果：
F=P.fit(x,y,deg=4)
print(F)
# 2.19829503 + 0.53756111·x - 1.53483146·x² + 0.32157981·x³ + 1.19567734·x⁴
print(F(0.6)) # 2.2488997156639994这种方式得到的在0.6处的拟合值结果为2.2488997156639994，两者几乎具有完全相同的结果。
事实上，这组数据是通过如下函数生成的（当然在拟合时，我们并不知道）：
该函数在该点处的结果为：
上图汇总，红色是原函数的曲线，而绿色是通过fit得到的多项式函数，而蓝色点则是使用neville方法计算的0.6处的值。虽然比较两种拟合方式比较接近，但是由于我们的数据不太符合多项式的特性，因此拟合还是出现了较大的偏差，因此，再次说明，多项式拟合并不一定能够构造出客观反映数据规律的多项式。
