Lagrange Interpolation

Concept

The basis function defined as

$$l_j(t_i)={\delta }_{ij}=\{\begin{array}{ll}1& ,i=j\\ 0& ,i\ne j\end{array}$$

where $t_i$ is given point then

$$p_{n-1}(t)=y_1{l}_1(t)+y_2{l}_2(t)+\cdot \cdot \cdot +y_n{l}_n(t)$$

Also we can define $l_j$ as

$$l_j(t)=\prod _{i\ne j}^n\frac{t-t_i}{t_j-t_i}$$

Codes

Lagrange interpolation

import numpy as np
 
def lagrange(t_vals, f_vals, t):
	"""
		`t_vals`: argument values used to construct the interpolants.
		`f_vals`: function values used to construct the interpolants.
		`t`: Values at which to evaluate the interpolant
	"""
	m = len(t_vals)
    res = np.empty_like(t)
    for i in range(len(t)):
        res[i] = 0
        for j in range(m):
            sum_up = 1
            sum_down = 1
            for k in range(m):
                if k != j:
                    sum_up *= t[i] - t_vals[k]
                    sum_down *= t_vals[j] - t_vals[k]
            res[i] += sum_up * f_vals[j] / sum_down
		
    return res
 

Fast Lagrange interpolation

import numpy as np
 
def update_w(w, old_m, m):
	for j in range(old_m):
        for k in range(old_m, m):
            if j != k:
                w[j] = w[j] / (t_vals[j] - t_vals[k])
    for j in range(old_m, m):
        w[j] = 1
        for k in range(m):
            if j != k:
                w[j] *= t_vals[j] - t_vals[k]
        w[j] = 1 / w[j]
    return w
 
def improved_lagrange(t_vals, f_vals, t):
	m = len(t_vals)
    res = np.empty_like(t)
    
    for i in range(len(t)):
        w = np.empty_like(f_vals)
		
        w = update_w(w, 0, m)
		
        right = 1
        for j in range(0, m):
            right *= t[i] - t_vals[j]
        s = 0
        for j in range(m):
            s += f_vals[j] * w[j] / (t[i] - t_vals[j])
        res[i] = right * s
    return res
 

barycentric Lagrange interpolation

import numpy as np
 
def update_w(w, old_m, m):
	for j in range(old_m):
        for k in range(old_m, m):
            if j != k:
                w[j] = w[j] / (t_vals[j] - t_vals[k])
    for j in range(old_m, m):
        w[j] = 1
        for k in range(m):
            if j != k:
                w[j] *= t_vals[j] - t_vals[k]
        w[j] = 1 / w[j]
    return w
 
def barycentric(t_vals, f_vals, t):
    res = np.empty_like(t)
    m = len(t_vals)
    
    for i in range(len(t)):
        w = np.empty_like(f_vals)
        w = update_w(w, 0, m)
        
        s_up = 0
        s_down = 0
        for j in range(m):
            s_up += w[j] * f_vals[j] / (t[i] - t_vals[j])
            s_down += w[j] / (t[i] - t_vals[j])
        res[i] = s_up / s_down
		
    return res
 

Lagrange with Newton table

import numpy as np
 
def lagrange(t_vals, f_vals, t):
    def newton_table(x, y, m):
        x = np.copy(x[:m])
        a = np.copy(y[:m])
        for k in range(1, m):
            a[k:m] = (a[k:m] - a[k - 1]) / (x[k:m] - x[k - 1])
        return a
	
    res = np.empty_like(t)
    m = len(t_vals)
    for i in range(len(t)):
        tab = newton_table(t_vals, f_vals, m)
		    
        right = np.empty_like(t_vals)
        right[0] = 1
        for iii in range(1, m):
	        right[iii] = right[iii - 1] * (t[i] - t_vals[iii - 1])
        res[i] = tab @ right[:m].T
		
    return res