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 [ ]: