In [2]:
import numpy as np
import numpy.linalg as la
In [16]:
a = 0
b = 1
n = 5
In [17]:
def get_quad_rule(a, b, n):
    x = np.linspace(a, b, n)
    V = np.zeros((n,n))

    for i in range(n):
        V[:, i] = x**i

    integrals = np.array([
            1/(i+1) * (b**(i+1)-a**(i+1))
            for i in range(n)])
    w = la.solve(V.T, integrals)
    
    return x, w
In [22]:
def get_quad_rule_cheap(a, b, n):
    x01, w01 = get_quad_rule(0, 1, n) # This could just be stored once and for all--for each value of n
    return a + x01 * (b-a), w01 * (b-a)
In [18]:
get_quad_rule(0, 1, 5)
Out[18]:
(array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ]),
 array([ 0.07777778,  0.35555556,  0.13333333,  0.35555556,  0.07777778]))
In [19]:
get_quad_rule(17, 18, 5)
Out[19]:
(array([ 17.  ,  17.25,  17.5 ,  17.75,  18.  ]),
 array([ 0.07777778,  0.35555556,  0.13333333,  0.35555556,  0.07777778]))
In [23]:
get_quad_rule(17, 19, 5)
Out[23]:
(array([ 17. ,  17.5,  18. ,  18.5,  19. ]),
 array([ 0.15555556,  0.71111111,  0.26666667,  0.71111111,  0.15555556]))
In [24]:
get_quad_rule_cheap(17, 19, 5)
Out[24]:
(array([ 17. ,  17.5,  18. ,  18.5,  19. ]),
 array([ 0.15555556,  0.71111111,  0.26666667,  0.71111111,  0.15555556]))
In [ ]: