#!/usr/bin/env python3

# Runge's Phenomenon
#----------------------------------

import numpy as np

def func_to_interpolate(x):
    return 1/( 1 + 25 * x**2 )

def equidistant_nodes(n):
    return 2*np.arange(0,n+1)/n -1

def chebychev_nodes(n):
    return np.cos( (2*np.arange(1,n+1) -1) * np.pi / (2*n) )

def random_nodes(n):
    return 2*np.random.rand(n) -1

def main():
    import matplotlib.pyplot as plt
    from scipy.interpolate import lagrange as scipy_lagrange
    
    np.random.seed(0)


    n_special = 20
    n_cheby = n_special
    n_equi = n_special
    n_random = n_special

    n_func = 1e4
    x_func = np.linspace(-1, 1, n_func, endpoint = True )

    fig, ax = plt.subplots()

    ax.set_title("Lagrange Interpolations with different types of nodes")
    ax.grid()
    ax.set_xlabel("x")
    ax.set_ylabel("y")

    func = func_to_interpolate
    # Plot the formula
    ax.plot(x_func, func(x_func), '--', label="Function $1/(1+25x^2)$")

    # Using equidistant nodes
    x_equi = equidistant_nodes(n_equi)
    ax.plot(x_equi, scipy_lagrange(x_equi, func(x_equi))(x_equi), '^', label="Equidistant (n={})".format(n_equi))

    # Using chebychev nodes
    x_cheby = chebychev_nodes(n_cheby)
    ax.plot(x_cheby, scipy_lagrange(x_cheby, func(x_cheby))(x_cheby), '>',label="Chebyshev (n={})".format(n_cheby))

    # Using random nodes
    x_random = random_nodes(n_random)
    ax.plot(x_random, scipy_lagrange(x_random, func(x_random))(x_random), '.', label="Random (n={})".format(n_random))


    ax.legend()


    print("""
    Runge's Phenomenon
    ------------------

    According to the docs [1], the scipy Lagrange
    interpolation is unstable and there should not be more
    than 20 points.


    The region where the equidistant nodes break is near
    x = -1, where the interpolation is fluctuating between
    large negative and large positive values. From n > 25,
    this starts to get visible.

    The chebyshev nodes break near x = 1 when n > 30, where
    it interpolates the function to be negative.

    Funnily enough, the equidistant nodes vary much more
    than the chebyshev nodes.

    For the random nodes, you just have be lucky that not
    too many points are close together at the extremes of
    the interval.

    [1] = https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.lagrange.html
""")

    plt.show()


if __name__ == "__main__":
    main()