2
$\begingroup$

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
$\endgroup$
2
$\begingroup$

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
$\endgroup$
  • $\begingroup$ thanks again!!! $\endgroup$ – fr14 Dec 4 '18 at 1:59
  • $\begingroup$ I posted my final question, this one has to do with simpsons rule but I am running into an error again with my section function $\endgroup$ – fr14 Dec 4 '18 at 2:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.