Skip to content

Instantly share code, notes, and snippets.

@aurelienpierre
Last active December 8, 2022 10:44
Show Gist options
  • Save aurelienpierre/1d9826e7db078e048bf437e516a7a4b2 to your computer and use it in GitHub Desktop.
Save aurelienpierre/1d9826e7db078e048bf437e516a7a4b2 to your computer and use it in GitHub Desktop.
[Python] Computes the symbolic expression of the Lagrange interpolation polynomial and evaluates the function if needed
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)
@sjodiel
Copy link

sjodiel commented Jan 12, 2018

Hello, is it a version for python 2 or python3?

@sjodiel
Copy link

sjodiel commented Jan 12, 2018

Hello, The following error occurs: raise TypeError ("can not convert expression to float") in L[i] = numerator / denominator.

@Khalilsqu
Copy link

Khalilsqu commented Mar 26, 2018

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

@shelddin
Copy link

shelddin commented Dec 8, 2022

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment