Python code for Lagrange interpolation - determining the equation of the polynomial

Question

The following code takes in a single value, x, and a list of points, X, and determines the value of the Lagrange polynomial through the list of points at the given x value.

def chunkIt(seq, num):
    avg = len(seq) / float(num)
    out = []
    last = 0.0
    while last < len(seq):
        out.append(seq[int(last):int(last + avg)])
        last += avg
    return out

def product(list):
    p = 1
    for i in list:
        p *= i
    return p

def Lagrange(x,X):
    T = np.zeros((2,len(X)))
    list = []
    for i in range(len(X)):
        for j in range(len(X)):
            if i != j:
                list.append((x-X[j][0])/(X[i][0]-X[j][0]))
    p = []
    for i in chunkIt(list,len(X)):
        p.append(product(i))
    for i in range(len(X)):
        T[0][i] = p[i]
        T[1][i] = X[i][1]

    list2 = []
    for i in range(len(X)):
        list2.append(T[0][i]*T[1][i])
    return sum(list2)

For example:

x, X = 3, [[0,0],[1,1],[2,0.5]]

gives a value of -1.5.

How do I modify this code to determine the equation of the polynomial through the list of points? i.e. if I put x = 'x' as the input, I want it to return -0.75x**2 + 1.75x [for the given example]

Solution

import numpy as np
from pypoly import Polynomial

x, X = 3, [[0, 0], [1, 1], [2, 0.5]]

order = len(X)

This is the order of the resulting Lagrange polynomial. For your example, order is 3.

equations = np.array([[point[0] ** i for i in range(order)] for point in X])
values = np.array([point[1] for point in X])
coefficients = np.linalg.solve(equations, values)

This sets up simultaneous equations by substituting the points into a general polynomial. For order 3, the general polynomial is:

a * x ** 2 + b * x ** 1 + c * x ** 0 = y

It solves the system of simultaneous equations to find coefficients. For order 3, we get the values of a, b, c.

print 'coefficients', list(coefficients)

coefficients [0.0, 1.75, -0.75]

p = Polynomial(*coefficients)

Here, the * operator splits the elements of the array-like into individual values to be passed as arguments to Polynomial().

print p

1.75 * X - 0.75 * X**2

print p(x)

-1.5

To install PyPolynomial with pip, use:

for Python 2:

pip install PyPolynomial

for Python 3:

pip3 install PyPolynomial

Comments

[hiro protagonist]

maybe this is a place to start: sympy.org