我正在尝试使用牛顿插值多项式来计算以下数组的有限除法，以确定 x=8 处的 y。数组是
x = 0 1 2 5.5 11 13 16 18
y= 0.5 3.134 5.9 9.9 10.2 9.35 7.2 6.2
我拥有的伪代码位于http://imgur.com/gallery/Lm2KXxA/new。是否有任何可用的伪代码、算法或库可以用来告诉我答案？
我也相信这是如何在 matlab http://imgur.com/gallery/L9wJaEH/new中执行程序。我只是不知道如何在python中做到这一点。
尝试这个
def _poly_newton_coefficient(x, y):
"""
x: list or np array contanining x data points
y: list or np array contanining y data points
"""
m = len(x)
x = np.copy(x)
a = np.copy(y)
for k in range(1, m):
a[k:m] = (a[k:m] - a[k - 1])/(x[k:m] - x[k - 1])
return a
def newton_polynomial(x_data, y_data, x):
"""
x_data: data points at x
y_data: data points at y
x: evaluation point(s)
"""
a = _poly_newton_coefficient(x_data, y_data)
n = len(x_data) - 1 # Degree of polynomial
p = a[n]
for k in range(1, n + 1):
p = a[n - k] + (x - x_data[n - k])*p
return p
这是Python代码。该函数coef计算有限除差系数，并且该函数Eval评估给定节点处的插值。
import numpy as np
import matplotlib.pyplot as plt
def coef(x, y):
'''x : array of data points
y : array of f(x) '''
x.astype(float)
y.astype(float)
n = len(x)
a = []
for i in range(n):
a.append(y[i])
for j in range(1, n):
for i in range(n-1, j-1, -1):
a[i] = float(a[i]-a[i-1])/float(x[i]-x[i-j])
return np.array(a) # return an array of coefficient
def Eval(a, x, r):
''' a : array returned by function coef()
x : array of data points
r : the node to interpolate at '''
x.astype(float)
n = len( a ) - 1
temp = a[n] + (r - x[n])
for i in range( n - 1, -1, -1 ):
temp = temp * ( r - x[i] ) + a[i]
return temp # return the y_value interpolation
如果您不想使用任何函数定义，那么这里是简单的代码：
## Newton Divided Difference Polynomial Interpolation Method
import numpy as np
x=np.array([0,1,2,5.5,11,13,16,18],float)
y=np.array([0.5, 3.134, 5.9, 9.9, 10.2, 9.35, 7.2, 6.2],float)
n=len(x)
p=np.zeros([n,n+1])#creating a Tree table (n x n+1 array)
value =float(input("Enter the point at which you want to calculate the value of the polynomial: "))
# first two columns of the table are filled with x and y data points
for i in range(n):
p[i,0]=x[i]
p[i,1]=y[i]
## algorithm for tree table from column 2 two n+1
for i in range(2,n+1): #column
for j in range(n+1-i):# defines row
p[j,i]=(p[j+1,i-1]-p[j,i-1])/(x[j+i-1]-x[j])#Tree Table
np.set_printoptions(suppress=True)## this suppress the scientific symbol(e) and returns values in normal digits
# print(p) ## can check the complete Tree table here for NDDP
b=p[0][1:]#This vector contains the unknown coefficients in the polynomial which are the top elements of each column.
print("b= ",b)
print("x= ",x)
lst=[] # list where we will append the values of prouct terms
t=1
for i in range(len(x)):
t*=(value-x[i]) ##(x-x0), (x-x0)(x-x1), (x-x0)(x-x1)(x-x2) etc..
lst.append(t)
print("The list of product elements ",lst,end = ' ')
## creating a general function
f=b[0]
for k in range(1,len(b)):
f+=b[k]*lst[k-1] ## important**[k-1]** not k because in list we use one step earlier element. For example for b1 we have to use (x-x0), for b2, we use (x-x0)(x-x1) for b3 we use (x-x0)(x-x1)(x2)
print("The value of polynomial: ","%.3f"%f)
@Ledruid 提出的解决方案是最优的。它不需要二维数组。划分差异树在某种程度上类似于 2D 数组。一种更原始的方法（可以从中获得@Leudruid 的算法）是考虑一个“矩阵 $F_{ij}$ 对应于树的节点。这样算法就可以了。
import numpy as np
from numpy import *
def coeficiente(x,y) :
''' x: absisas x_i
y : ordenadas f(x_i)'''
x.astype(float)
y.astype(float)
n = len(x)
F = zeros((n,n), dtype=float)
b = zeros(n)
for i in range(0,n):
F[i,0]=y[i]
for j in range(1, n):
for i in range(j,n):
F[i,j] = float(F[i,j-1]-F[i-1,j-1])/float(x[i]-x[i-j])
for i in range(0,n):
b[i] = F[i,i]
return np.array(b) # return an array of coefficient
def Eval(a, x, r):
''' a : retorno de la funcion coeficiente()
x : abcisas x_i
r : abcisa a interpolar
'''
x.astype(float)
n = len( a ) - 1
temp = a[n]
for i in range( n - 1, -1, -1 ):
temp = temp * ( r - x[i] ) + a[i]
return temp # return the y_value interpolation