Interpolate Lagrange

##########################################
#
#   ЧИСЛЕННЫЕ МЕТОДЫ. ПОСТРОЕНИЕ И АНАЛИЗ
#   ЗАДАНИЕ №1
#
##########################################

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# vim:fileencoding=utf-8

import matplotlib.pyplot as plt
from math import sqrt, sin, log, ceil

# вычисление значения функции f_k(x) в точке x
def f (k: int, x: float):
    if sin(k * x) > 0:
        return log(sin(k*x)) / k / k / k
    else:
        return 0


# интерполяция методом Лагранжа (вычисление PL(x0))
def PL (n0 : int, x : [float], y : [float], x0 : float):
    y0 = 0.0

    for i in range(n0):
        l = 1.0

        for j in range(n0):
            if i == j:
                continue
            else:
                l *= (x0 - x[j]) / (x[i] - x[j])

        y0 += l*y[i]

    return y0

# заданная точность EPS и N = N(EPS) - число итераций при суммировании ряда
EPS = 1e-5
N = ceil(sqrt(1 / 2 / EPS))

# отрезок [a, b] и число точек n0
n0 = 11
a = 1; b = 2

# шаг изменения точек x из [a, b]
h = (b - a) / (n0 - 1)

# узлы интерполирования и значения функции
x = [a + i*h for i in range(n0)]
y = [0.0 for i in range(n0)]

# подсчет значений
for i in range(n0):
    for j in range(1, N + 1):
        y[i] += f(j, x[i])

# создаём координатную плоскость
fig = plt.subplots()
# узлы интерполирования (точки (x, y)) изображаются на рисунке
plt.scatter(x, y)

# увеличиваем число точек чтобы графики были точными и красивыми
# ищем новое число точек таким образом, чтобы оно было больше заданного числа count
count = 101; n = n0

while n < count:
    n += n - 1
    h /= 2.0

# точки, по которым будет строиться графики

xnew = [a + i*h for i in range(n)]        # координаты x
ynew1 = [PL(n0, x, y, x0) for x0 in xnew] # значения PL(x), x <- xnew
ynew2 = [0.0 for i in range(n)]           # значения f(x), x <- xnew

for i in range(n):
    for j in range(1, N + 1):
        ynew2[i] += f(j, xnew[i])

# рисуем графики (соединяя точки прямыми)

plt.plot(xnew, ynew1) # график многочлена Лагранжа
plt.plot(xnew, ynew2) # график функции

# показываем график
plt.show()