Lagrange interpolating polynomial for both arbitrary and equidistant points

lagrange_interpolation.py

from typing import Callable

import numpy as np
import matplotlib.pyplot as plt


def lagrangian_interpolation(xs: np.ndarray, ys: np.ndarray) -> Callable:
    """Make a Lagrangian interpolating polynomial function."""
    def _interpolate(x: float) -> float:
        result = 0
        for j in range(xs.shape[0]):
            prod = 1
            for k in range(xs.shape[0]):
                if k != j:
                    prod *= (x - xs[k]) / (xs[j] - xs[k])
            result += ys[j] * prod
        return result

    return _interpolate

xs = np.array([1, 1.5, 2, 4])
ys = np.array([4, 4.1, 2.5, 0.9])
f = np.vectorize(lagrangian_interpolation(xs, ys))
plt.plot(xs, ys, 'o')
xs1 = np.arange(-1, 5, 0.01)
plt.plot(xs1, f(xs1), '-')
plt.show()


def equidistant_lagrangian_interpolation(xs: np.ndarray, ys: np.ndarray) -> Callable:
    """Make a Lagrangian interpolating polynomial function given equidistant
    points."""
    x_min, x_max = np.min(xs), np.max(xs)
    h = (x_max - x_min) / xs.shape[0]
    
    def _interpolate(x: float) -> float:
        result = 0
        s = (x - x_min) / h
        for j in range(xs.shape[0]):
            prod = 1
            for k in range(xs.shape[0]):
                if k != j:
                    prod *= (s - k) / (j - k)
            result += ys[j] * prod
        return result + x_min * h

    return _interpolate

xs = np.array([1, 2, 3, 4])
ys = np.array([4, 4.1, 2.5, 0.9])
f = np.vectorize(lagrangian_interpolation(xs, ys))
plt.plot(xs, ys, 'o')
xs1 = np.arange(-1, 5, 0.01)
plt.plot(xs1, f(xs1), '-')
plt.show()