Newton’s interpolating polynomial [python]

Question

I am trying to compute the finite divided differences of the following array using Newton's interpolating polynomial to determine y at x=8. The array is

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

The pseudo code that I have is at https://i.sstatic.net/X3WLr.png. Are there any available pseudocode, algorithms or libraries I could use to tell me the answer?

Also I believe this is how to do the program in matlab https://i.sstatic.net/1MOuN.jpg. I just can't figure out how to do it in python.

Answer 1

Try this

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

Answer 2

If you don't want to use any function definition, then here is the simple code:

## 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)

Answer 3

Here is the Python code. The function coef computes the finite divided difference coefficients, and the function Eval evaluates the interpolation at a given node.

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

Answer 4

The solution proposed by @Ledruid is optimal. It does not require 2 dimensional array. The tree of divided differences is in a way like a 2D array. A more primitive approach (from which @Leudruid's algorithm can be obtained) is thinking of a "matrix $F_{ij}$ corresponding to the nodes of the tree. In this way the algorithm would be.

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