import numpy as np
import matplotlib.pyplot as plt

def runge(x):
    return 1 / (1 + 25 * x**2)

def divided_differences(x, y):
    n = len(y)
    a = np.zeros((n,n))
    
    a[:,0] = y
    
    for i in range(1,n):
        for j in range(1,i+1):
            a[i,j] = (a[i,j-1] - a[i-1,j-1]) / (x[i] - x[i-j]) # (links daneben - links darüber) / (x[zeile] - x[zeile - spalte])
    
    return a
    
def newton_interpolation(a, data, x):
    n = len(a)
    
    p = np.zeros(len(x))
    
    for i in range(1,n+1):
        p = a[n-i] + (x - data[n-i]) * p # Horner Schema; '-' and '*' are overloaded for numpy arrays: so at least one argument 'data' or 'x' has to be numpy array 
    
    return p

####################################################################################################################################################################
n = 12

x = np.linspace(-1, 1, 200)
x_e = np.linspace(-1, 1, n) # equidistant grid points

# Chebyshev grid points
x_c = np.zeros(n)
for i in range(0,n):
    x_c[i] = np.cos((2 * i + 1) * np.pi / (2 * n))
    
f = runge(x)
y_e = runge(x_e) # values for grid points for interpolation with equidistant grid points
y_c = runge(x_c) # values for grid points for interpolation with Chebyshev grid points

# Interpolation with equidistant grid points and evaluation of interpolated values at x
a_e = np.diag(divided_differences(x_e, y_e))
p_e = newton_interpolation(a_e, x_e, x)

# Interpolation with Chebyshev grid points and evaluation of interpolated values at x
a_c = np.diag(divided_differences(x_c, y_c))
p_c = newton_interpolation(a_c, x_c, x)
       
       
# Plotting of Runge and the two interpolated polynomials
plt.rcParams['text.usetex'] = True

plt.plot(x, f, label=r"$f(x) = \frac{1}{1 + 25 x^2}$")
plt.plot(x, p_e, label=r"$p_n(x)$ equidistant")
plt.plot(x, p_c, label=r"$p_n(x)$ Chebyshev")
plt.plot([], [], ' ', label="$n = {}$".format(n))

plt.title("A3)")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()

plt.show()

# braendel@math.tu-freiberg.de