minpoly.py

import numpy as np

def find_min_poly(x, d, n):
    """Find the polynomial that minimizes max|p(x)| over the given
       range of x values and among polynomials of the form

       x^n + a_d x^d + a_{d-1} x^{d-1} + ... + a_0

       Inputs:

        x: a vector of x values; conveniently generated with numpy.linspace
        d, n: integers

        I wrote this in response to a Mathoverflow question:
        http://mathoverflow.net/questions/215411/a-generalization-of-chebyshev-polynomials
    """
    from cvxopt import matrix
    from cvxopt.modeling import variable, op, max

    X = matrix(x)
    y = np.empty( (len(x), d+2) )
    y[:,0] = x**n
    for i in range(d, -1, -1):
        y[:, d+1-i] = x**i

    Y = matrix(y)
    c = variable(d+1)

    try:
        op(max(abs(Y[:,0] + Y[:,1:]*c))).solve()
    except(ValueError):
        return 0, 0

    return c.value, max(abs((Y[:,0] + Y[:,1:]*c.value)))

def shanks(x):
    """Apply the Shanks transformation to accelerate sequence convergence."""
    S = (x[2:]*x[:-2] - x[1:-1]**2) / (x[2:] - 2*x[1:-1] + x[:-2])
    return S

if __name__ == "__main__":
    x = np.linspace(-1,1,10000)
    N = range(1,51) # Actually fails for n bigger than about 25
    maxabs = np.zeros( (len(N),len(N)) )

    for n in N:
        for d in range(n-1,-1,-1):
            print n, d
            coeffs, M = find_min_poly(x, d, n)
            maxabs[d, n-1] = M
            print M

    print maxabs

    # Investigate conjecture that diagonal ratios converge to 2:
    diag = 2
    z = np.zeros(50)
    for i in range(diag+1,49):
        z[i] = maxabs[i-diag,i]/maxabs[i-diag+1,i+1]
    plt.plot(z)
    plt.plot(shanks(z))
    plt.plot(shanks(shanks(z)))
    # The evidence is inconclusive.