本篇说说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处的值。虽然比较两种拟合方式比较接近,但是由于我们的数据不太符合多项式的特性,因此拟合还是出现了较大的偏差,因此,再次说明,多项式拟合并不一定能够构造出客观反映数据规律的多项式。