from __future__ import division

import numpy as np
import numpy.linalg as la
import scipy.special as sps
import matplotlib.pyplot as pt

# Set parameters

n = 5
left_end_nodes = 0
right_end_nodes = 0

# Ignore (don't even evaluate) this cell at the outset.

# If enabled, it computes weights for variants of Gauss 
# quadratures that include interval ends.

# Enable (only) this for Gauss-Radau:
if 1:
    left_end_nodes = 1
    
# Enable (only) this for Gauss-Lobatto:
if 0:
    left_end_nodes = 1
    right_end_nodes = 1

deg_reduction = left_end_nodes+right_end_nodes

nodes = np.empty(n)
nodes[left_end_nodes:n-right_end_nodes] = \
    sps.legendre(n - deg_reduction).weights[:, 0] 

nodes[:left_end_nodes] = -1
nodes[n-right_end_nodes:] = 1

mesh = np.linspace(-1, 1, 300)

pt.plot(mesh, sps.eval_legendre(n - deg_reduction, mesh))
pt.plot(nodes, 0*nodes, "o")

# Use method of undetermined coefficients to find
# Newton-Cotes rule for Gauss nodes.

max_degree = len(nodes) - 1
powers = np.arange(max_degree+1)

Vt = nodes ** powers.reshape(-1, 1)

a, b = -1, 1
rhs = 1/(powers+1) * (b**(powers+1) - a**(powers+1))

weights = la.solve(Vt, rhs)

for i in range(2*len(nodes) + 1):
    approx = np.dot(weights, nodes**i)
    true = 1/(i+1)*(1. - (-1)**(i+1))
    
    print "Error at degree %d: %g" % (i, approx-true)