[Python] Computes the symbolic expression of the Lagrange interpolation polynomial and evaluates the function if needed

lagrange.py

from sympy import *
from sympy.matrices import *
import numpy
import pylab
import warnings

def Lagrange_interpolation(points, variable=None):
    """
    Compute the Lagrange interpolation polynomial.
    
    :var points: A numpy n×2 ndarray of the interpolations points
    :var variable: None, float or ndarray
    :returns:   * P the symbolic expression
                * Y the evaluation of the polynomial if `variable` is float or ndarray
                
    """
    x = Symbol("x")
    L = zeros(1, points.shape[0])
    i = 0

    for p in points:
        numerator = 1
        denominator = 1
        other_points = numpy.delete(points, i, 0)

        for other_p in other_points:
            numerator = numerator * (x - other_p[0])
            denominator = denominator * (p[0] - other_p[0])

        L[i] = numerator / denominator
        i = i+1
        
    # The Horner factorization will reduce chances of issues with floats approximations
    P = horner(L.multiply(points[..., 1])[0])
    Y = None
    
    try: 
        Y = lambdify(x, P, 'numpy')
        Y = Y(variable)
            
    except:
        warnings.warn("No input variable given - polynomial evaluation skipped")
            
    return P,Y

def test_Lagrange(sets):
    for points in sets:
        x = numpy.linspace(0, 100)
        P, Y = Lagrange_interpolation(points, x)
    
        print(P)
        pylab.plot(x, Y)
        

if __name__ == '__main__':

    sets = [numpy.array([ # Linear
                        [0,1],
                        [50, 50],
                        ]),
            numpy.array([ # Quadratic
                        [0,1],
                        [50, 50],
                        [100, 1]
                        ]),
            numpy.array([ # Cubic
                        [0,1],
                        [50, 50],
                        [75, 40],
                        [100, 1]
                        ])
           ]

    test_Lagrange(sets)

This is probably also good way

from sympy import Symbol, simplify, lambdify
import numpy as np
import matplotlib.pyplot as plt
from functools import reduce
import operator


def interpolate_lagrange(x, x_values, y_values):
    """
    x : value at which to evaluate y, should be between min and max x_values
    x_values: list or numpy array containing values of x
    y_values: list or numpy array contaning values of y
    """
    def _basis(j):
        p = [(x - x_values[m])/(x_values[j] - x_values[m]) for m in range(k) if m != j]
        return reduce(operator.mul, p)
    assert len(x_values) != 0 and (len(x_values) == len(y_values)), 'x and y cannot be empty and must have the same length'
    k = len(x_values)
    return sum(_basis(j)*y_values[j] for j in range(k)) 

x = Symbol('x')
poly = simplify(interpolate_lagrange(x,[-1, 0, 1, 2],[3,-4, 5, -6]))

x1 = np.linspace(-1, 2, 100)
y1 = lambdify(x, ploy)(x1)

fig, ax = plt.subplots()
ax.plot(x1, y1)
ax.scatter([-1, 0, 1, 2],[3,-4, 5, -6], c = 'r')
plt.show()

@Khalilsqu thank you for the code, you got a small typo in y1 definition. should be poly.