This Python script compares Chebyshev polynomials against polynomials pp = p / ||p||_2

integral_twonorm_comparisons.py

import hashlib

import matplotlib.pyplot as plt
import numpy as np

from scipy import integrate

from chebyshev import ChebyshevPolynomial

def pid(p):
    sha1 = hashlib.sha1()
    sha1.update(p.r.tobytes())
    return sha1.hexdigest()[:7]

    
def poly_repr(coeffs):
    if isinstance(coeffs, np.poly1d):
        coeffs = coeffs.coef
    terms = []
    n = len(coeffs) - 1 

    for i, coeff in enumerate(coeffs):
        # Skip zero coefficients
        if np.isclose(coeff, 0):
            continue
        
        # Determine the exponent
        exp = n - i
        
        # Format the coefficient to 2 decimal places or as an integer
        if coeff % 1 == 0:
            # Integer coefficient, convert to int for display
            formatted_coeff = str(int(coeff))
        else:
            # Non-integer, format to 2 decimal places
            formatted_coeff = f"{coeff:.2f}"
        
        # Construct the monomial term
        if exp == 0:
            term = formatted_coeff
        elif exp == 1:
            term = f"{formatted_coeff}z"
        else:
            term = f"{formatted_coeff}z^{exp}"
        
        # Handle positive coefficients (except the first)
        if coeff > 0 and i != 0:
            term = "+ " + term
        
        terms.append(term)

    return " ".join(terms)

def twonorm(p):
    """ Compute the two norm  \int_{-1}^1 p(x)^2 dx/sqrt{1-x^2}
        
        This uses scipy.quad, which uses the Fortran library QUADPACK
        to perform adaptive quadrature numerical integration
    """
    return np.sqrt(integrate.quad(lambda x: p(x)**2 / np.sqrt(1 - x**2), -1, 1))


def twonorm(p):
    """ Compute the two norm  \int_{-1}^1 p(x)^2 dx/sqrt{1-x^2}
        
        This uses scipy.quad, which uses the Fortran library QUADPACK
        to perform adaptive quadrature numerical integration
    """
    return np.sqrt(integrate.quad(lambda x: (p(x)*p(x) / np.sqrt(1 - x**2)), -1, 1))


# Set up T, points on the inner circles
T = ChebyshevPolynomial.classical(2)
circle_points = T.Ek_circle_points(n=100_000)
# Compute a grid for plotting over in case we find counterexamples
xv, yv, zv = T.grid()
# and normalize the monic Chebyshev polynomial in the integral 2norm
T2norm, error = twonorm(T.T) # T.T is monic
T2 = T.T/T2norm

# Evaluate T on the inner circles
Td = np.abs(T2(circle_points))

# Evaluate T over the grid
Tzv = np.abs(T2(zv))

for i in range(10):
    # Every 10 comparisons, print a message so that we can track status
    if not i % 10:
        print(f"Checked {i} polynomials")
        
    # Generate random quadratic polynomials by picking 2 roots uniformly on [-1, 1]
    roots = np.random.uniform(-1, 1, 2)
    p = np.poly1d(roots, r=True)
    
    # Compute the twonorm of p
    norm, error = twonorm(p)
    # Normalize p with the 2 norm we just computed
    pp = (p/norm)
    print("Here is the norm/error of the original polynomial along with the norm of the normalized polynomial as a sanity check", (norm, error), twonorm(pp))
    # Evaluate p on the circle points (inner disks)
    ppd = np.abs(pp(circle_points))
    if any(ppd > Td):
        # If we found a counterexample, print a message and generate a plot
        print(f"Found a possible counterexample: {pp=} with roots at {pp.r=}")
        fig, ax = plt.subplots()
        ppzv = np.abs(pp(zv))
        print(twonorm(p)[0], twonorm(pp)[0])
        ax.contourf(xv, yv, Tzv > ppzv)
        T.plot_disks(ax=ax)
        fig.suptitle(f"Locations where $|T_2(z)| < |p(z)|$ (blue), $p(z) \\approx " + poly_repr(pp) + "$")
        fig.savefig(f"counterexample-2norm-{pid(pp)}.png", facecolor='white', transparent=False)