import numpy as np
import numpy.linalg as la
a = 0
b = 1
n = 5
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
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)
get_quad_rule(0, 1, 5)
get_quad_rule(17, 18, 5)
get_quad_rule(17, 19, 5)
get_quad_rule_cheap(17, 19, 5)