Problems integrating using quadpy

Question

I am trying to use quadpy as I want to do 2D numerical vector integration for 2D integrals. To see how fast quadpy works, I wanted to test it and compare it with scipy 1D vector integration. Thus I wrote a simple code:

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import quadpy

q = np.linspace(0.03, 1.0, 500)


def f(t):
    return t * 0.5 * (erf((t - 40) / 3) - 1) * j0(q * t)


y, _ = integrate.quad_vec(f, 0, 50)
y1, _ = quadpy.quad(f, 0, 50)

print(y)
print(y1)

Having no experience with quadpy, I get the following errors:

Traceback (most recent call last):
  File "test.py", line 8, in <module>
    y1 = quadpy.quad(lambda t: t*0.5*(erf((t-40)/3)-1)*j0(q*t),0.0,50)
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/_scipy_compat.py", line 44, in quad
    g, [a, b], eps_abs=epsabs, eps_rel=epsrel, max_num_subintervals=limit
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/line_segment/_tools.py", line 42, in integrate_adaptive
    range_shape=range_shape,
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/line_segment/_gauss_kronrod.py", line 124, in _gauss_kronrod_integrate
    fx_gk = numpy.asarray(f(sp))
  File "/Users/shankardutt/opt/anaconda3/envs/Cylinder_Fitting/lib/python3.7/site-packages/quadpy/_scipy_compat.py", line 41, in g
    return f(x, *args)
  File "test.py", line 8, in <lambda>
    y1 = quadpy.quad(lambda t: t*0.5*(erf((t-40)/3)-1)*j0(q*t),0.0,50)
ValueError: operands could not be broadcast together with shapes (500,) (15,)

Change your lambda to a regular def. Then check the t that quad passes it. I suspect it will be a (15,) vector, which won't compute with your (500,) shaped q. As written your function only work with a scalar t. Also study the quadpy docs about how to properly use vector valued functions. With the scipy package you had to use quad_vec not the regular quad. — hpaulj
– hpaulj, Commented Mar 28, 2020 at 20:21
@hpaulj Thanks a lot for your comment. I definitely tried using def and got the same error. After your comment I checked and you are right t is a (15,) vector. I cannot find much in quadpy docs. Again, regarding the comment that my function can work with scalar t, I am not sure how to implement that. Can you help with the same? And regarding scipy, I used quad_vec instead of quad inside a for loop because the former is quite fast. I was hoping to get some good speed from quadpy but I failed. — Shankar_Dutt
– Shankar_Dutt, Commented Mar 28, 2020 at 21:11
Looks like it's trying to gain speed by asking for values at 15 t at once. t*q[:,None] would give an outer product (15,500). But it's up to you to figure out how to use that in a meaningful way. — hpaulj
– hpaulj, Commented Mar 28, 2020 at 21:20

Nico Schlömer · Accepted Answer · 2020-03-29 12:33:42Z

2

You're missing one multiply.outer:

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import quadpy

q = np.linspace(0.03, 1.0, 500)


def f(t):
    return t * 0.5 * (erf((t - 40) / 3) - 1) * j0(np.multiply.outer(q, t))


y, _ = integrate.quad_vec(f, 0, 50)
y1, _ = quadpy.quad(f, 0, 50)

print(y - y1)

answered Mar 29, 2020 at 12:33

Nico Schlömer

59.7k36 gold badges216 silver badges291 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Shankar_Dutt Over a year ago

Thanks a lot for your reply. I have an additional question, how would one implement 2D integral using quadpy. Please see this stackoverflow.com/questions/61026523/…

Collectives™ on Stack Overflow

Problems integrating using quadpy

1 Answer 1

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related