Interpolation (pol, Lagrange, Newton)

import numpy as np
import matplotlib.pyplot as plt

def main():
    #максимальный индекс узловой точки
    '''
    n = int(input("n = "))

    #таблица узлов
    nodes_x = np.empty(n+1)
    nodes_f = np.empty(n+1)
    for i in range(n+1):
        nodes_x[i] = int(input("x[" + str(i) + "] ="))
        nodes_f[i] = int(input("f_" + str(i) + " = "))
    '''
    n = 4
    nodes_x = np.array([1, 2, 3, 4, 5])
    nodes_f = np.array([1, 16, 81, 256, 625])
    print("\nУзловые точки\n")
    print(nodes_x, end='\n')
    print(nodes_f, end='\n')

    #polynomial(nodes_x, nodes_f, n)

    #построение неузловых точек для формы Лагранжа
    x_to_plot = np.linspace(nodes_x[0], nodes_x[n], n*2+1)
    y_to_plot = [lagrange(nodes_x, nodes_f, i) for i in x_to_plot]
    print('\nЗначения функции в неузловых точках, согласно интерполянте в форме Лагранжа\n')
    print(x_to_plot)
    print(y_to_plot)
    plt.plot(nodes_x, nodes_f, 'o', x_to_plot, y_to_plot)
    plt.suptitle("Интерполирующий многочлен")
    plt.grid(True)
    plt.show()

    y_to_plot.clear()

    y_to_plot = [newton(nodes_x, nodes_f, i, n) for i in
                 x_to_plot]
    print("\nЗначение функции в неузловых точках, согласно интерполянте в форме Ньютона\n")
    print(x_to_plot)
    print(y_to_plot)

####################################################################
def polynomial(nodes_x, nodes_f, n):
    #матрица иксов для слау
    matrix = np.zeros((n+1, n+1))
    for i in range(n+1):
        for j in range(n+1):
            matrix[i][j] = np.power(nodes_x[i], n+1-j)
    print (matrix, "@")

    #вектор-столбец для слау
    vector = np.zeros((n+1, 1))
    for i in range(n+1):
        vector[i] = nodes_f[i]
    print (vector, "@")
####################################################################

####################################################################
def lagrange(x, f, t):
    L = 0
    for i in range(len(f)):
        const_1 = 1
        const_2 = 1
        for j in range(len(x)):
            if j != i:
                const_1 *= (t - x[j])
                const_2 *= (x[i] - x[j])
        L += (f[i] * const_1 / const_2)
    return L
####################################################################

####################################################################
def newton(x, y, r, n):
    #x — узловые точки
    #y — значения функции в узловых точках
    #r — точка x, для которой находится значение функции
    x.astype(float)
    y.astype(float)
    s = y[0]
    for k in range(1, n+1):
        temp = 1
        for i in range(k):
            temp *= (r - x[i])
        temp *= fk(x, y, k)
        s += temp
    return s
####################################################################

####################################################################
#разделенные разности
def fk(x, y, k):
    x.astype(float)
    y.astype(float)
    ans = 0
    for i in range(k+1):
        d = 1
        for j in range(k+1):
            if j != i:
                d *= (x[i] - x[j])
        ans += y[i] / d
    return ans
####################################################################



if __name__ == "__main__":
    main()