Calculate the lebesgue constant of a polynomial interpolation operator.

lebesgue.py

import math

n = 20

#x_k = [-1 + 2 * k / float(n) for k in range(0, n + 1)]

x_k = [math.cos((2.0 * k + 1) / (2.0 * n +2.0) * math.pi)
       for k in range(0, n + 1)]

def l_j(j, x):
    product = 1.0
    for i in range(0, n + 1):
       if i != j:
           product *= (x - x_k[i]) / (x_k[j] - x_k[i])
    return product

def y(x):
    csum = 0
    for j in range(0, n + 1):
        csum += abs(l_j(j, x))
    return csum

def drange(start, stop, step):
    r = start
    while r < stop:
        yield r
        r += step

def d():
    greatest = None
    for x in drange(-1.0, 1.0, 0.000001):
        current = y(x)
        if greatest is None or current > greatest:
            greatest = current
    return greatest

print d()