1

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,)
Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
Shankar_Dutt
  • 115
  • 9
  • 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 Mar 28 '20 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 Mar 28 '20 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 Mar 28 '20 at 21:20

1 Answers1

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)
Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
  • Thanks a lot for your reply. I have an additional question, how would one implement 2D integral using quadpy. Please see this https://stackoverflow.com/questions/61026523/what-would-be-the-computationally-faster-way-to-implement-this-2d-numerical-inte?noredirect=1#comment107965006_61026523 – Shankar_Dutt Apr 04 '20 at 20:41