Gauss Quadrature in numerical methods

Question

In the following code I have implemented gauss quadrature. It is working correctly for my first function but for my second function I am getting an error.

So,

Is it possible to do this code without sing p_roots? and if not how do I fix it so my second function works?

from scipy.special.orthogonal import p_roots
def gaussquad(f,a,b,n):
    [x,w] = p_roots(n+1)
    G = 0.5*(b-a)*sum(w*f(0.5*(b-a)*x+0.5*(b+a)))
    return G
print(gaussquad(lambda x: x**2,5,10,100))
print(gaussquad(lambda x: sin(x),0,pi/2,100))

giving the output

291.6666666666665

---------------------------------------------------------------------------
TypeError                                 Traceback (most   recent call last)
<ipython-input-22-b4248b816207> in <module>()
  5     return G
  6 print(gaussquad(lambda x: x**2,5,10,100))
----> 7 print(gaussquad(lambda x: sin(x),0,pi/2,100))

<ipython-input-22-b4248b816207> in gaussquad(f, a, b, n)
  2 def gaussquad(f,a,b,n):
  3     [x,w] = p_roots(n+1)
----> 4     G = 0.5*(b-a)*sum(w*f(0.5*(b-a)*x+0.5*(b+a)))
  5     return G
  6 print(gaussquad(lambda x: x**2,5,10,100))

<ipython-input-22-b4248b816207> in <lambda>(x)
  5     return G
  6 print(gaussquad(lambda x: x**2,5,10,100))
----> 7 print(gaussquad(lambda x: sin(x),0,pi/2,100))

TypeError: only size-1 arrays can be converted to Python scalars

Answer 1

Vectorization issues. Use numpy instead

import numpy 
from scipy.special.orthogonal import p_roots

def gaussquad(f,a,b,n):
    [x,w] = p_roots(n+1)
    G = 0.5*(b-a)*sum(w*f(0.5*(b-a)*x+0.5*(b+a)))
    return G
print(gaussquad(lambda x: x**2,5,10,100))
print(gaussquad(lambda x: numpy.sin(x),0,numpy.pi/2,100))

Result is

1.000000000000001