Wednesday, April 12, 2017

Newton Polynomial in Python

I wrote code for this in Delphi. This time I want to rewrite it in Python based on this wiki.

I use this set of data point
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)

and I use xc=3 for the test data.

It's obvious that these sets of data points have quadratic form and f(xc) must have value of 9.




The heart of code lay on this




def n(j,xc,x):
    n = 1
    for i in arange(j):
        n *= (xc-x[i])

    return n
def a(j,l,x,y):
    if j==0:
        return y[0]
    elif j-l==1 :
        return (y[j]-y[l])/(x[j]-x[l])
    else:
        return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])
        

def N(xc,x,y):
    N = 0
    for j in range(len(x)):
        N += a(j,0,x,y)*n(j,xc,x)
    
    return N


Look at the function a(j,l,x,y), that's recursive function to obtain Newton divided difference value.

the whole code is this
from pylab import *

def n(j,xc,x):
    n = 1
    for i in arange(j):
        n *= (xc-x[i])

    return n
def a(j,l,x,y):
    if j==0:
        return y[0]
    elif j-l==1 :
        return (y[j]-y[l])/(x[j]-x[l])
    else:
        return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])
        

def N(xc,x,y):
    N = 0
    for j in range(len(x)):
        N += a(j,0,x,y)*n(j,xc,x)
    
    return N

x = []
y = []
#initial value
x.append(0)
x.append(1)
x.append(2)
x.append(4)
x.append(5)

y.append(0)
y.append(1)
y.append(4)
y.append(16)
y.append(25)

#for testing
xc = 3

yc = N(xc,x,y)

print ''
print xc, yc
#plot
t = linspace(-7,7,100)
u = N(t,x,y)
plot(t,u)
grid(True)
show()



.


here's the graphics


with another sets of data points, I have another result
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)