Python数值计算（14）——Neville方法

C囧囧

已于 2024-08-03 17:02:50 修改

阅读量665

点赞数 4

分类专栏： Python 文章标签：算法 numpy python

于 2024-07-28 18:17:21 首次发布

本文链接：https://blog.youkuaiyun.com/cdinten/article/details/140751714

版权

Python 专栏收录该内容

45 篇文章

订阅专栏

本篇说说Neville方法。Neville方法的基础是，插值多项式可以递归的生成，有时进行插值的目的是为了计算某个点的值，这个时候并不需要将拟合曲线完全求出，而是可以通过递归的方式进行计算，具体操作如下：

例如有多个点，为了表达方便，这里取五个点 $x_0,x_2,...x_4$ ，与之对应的值为 $P_0,P_1,...P_4$ ，则可以构造列表：

其中

$P_{0,1}=\frac{(x-x_0)P_1+(x-x_1)P_0}{x_1-x_0}$

$P_{1,2}=\frac{(x-x_1)P_2+(x-x_2)P_1}{x_2-x_1}$

等等，随后 $P_{0,1,2}$ 计算方式为：

$P_{0,1,2}=\frac{(x-x_0)P_{1,2}+(x-x_2)P_{0,1}}{x_2-x_0}$

最后递归地计算到 $P_{0,1,2,3,4}$ ，该值就是在对应值为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对象。

如何验证其拟合值呢？

第一组，用原来的函数 $f(x)=1-2x+3x^3$ ，数据为：

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，两者几乎具有完全相同的结果。

事实上，这组数据是通过如下函数生成的（当然在拟合时，我们并不知道）：

$f(x)=\frac{e^x}{x^2-x+1}$

该函数在该点处的结果为： $f(0.6)=2.3975247373559325$

上图汇总，红色是原函数的曲线，而绿色是通过fit得到的多项式函数，而蓝色点则是使用neville方法计算的0.6处的值。虽然比较两种拟合方式比较接近，但是由于我们的数据不太符合多项式的特性，因此拟合还是出现了较大的偏差，因此，再次说明，多项式拟合并不一定能够构造出客观反映数据规律的多项式。