The doit function for integrals is expected to be used to evaluate integrals symbolically:
In [47]: f = Integral(exp(-x**2), (x, 0, oo))
In [48]: f
Out[48]:
∞
⌠
⎮ 2
⎮ -x
⎮ ℯ dx
⌡
0
In [49]: f.doit()
Out[49]:
√π
──
2
In [50]: _.evalf()
Out[50]: 0.886226925452758
When that works it's nice because an explicit expression without involving integral signs that can then be evaluated more efficiently numerically. You can do that evaluation with evalf for multiprecision calculations or you can use lambdify and have numpy/scipy etc do the calculation much faster in machine precision floating point.
However not all integrals can be evaluated symbolically so doit will sometimes just return the same integral:
In [53]: h2 = Integral(2*r*exp(I*r**2*u/2)*besselj(0, r*v), (r, 0, 1))
In [54]: h2
Out[54]:
1
⌠
⎮ 2
⎮ ⅈ⋅r ⋅u
⎮ ──────
⎮ 2
⎮ 2⋅r⋅ℯ ⋅besselj(0, r⋅v) dr
⌡
0
In [55]: h2.doit() # slow
Out[55]:
1
⌠
⎮ 2
⎮ ⅈ⋅r ⋅u
⎮ ──────
⎮ 2
2⋅⎮ r⋅ℯ ⋅besselj(0, r⋅v) dr
⌡
0
Here we see that the integral did not evaluate. Sometimes where an integral has symbols other than the integration variable we can get a symbolic result by substituting values in. That's what happens here if you set either u or v to zero:
In [59]: h2.subs({u:0, v:5}).doit()
Out[59]:
2⋅besselj(1, 5)
───────────────
5
In [60]: h2.subs({u:5, v:0}).doit()
Out[60]:
5⋅ⅈ
───
2
2⋅ⅈ⋅ℯ 2⋅ⅈ
- ──────── + ───
5 5
Because the integrand became simpler it was possible to evaluate the integrals symbolically. However in this case if neither u or v is zero then the integral cannot be evaluated symbolically:
In [61]: h2.subs({uu:2,vv:3}).doit()
Out[61]:
1
⌠
⎮ 2
⎮ ⅈ⋅r
2⋅⎮ r⋅ℯ ⋅besselj(0, 3⋅r) dr
⌡
0
I expect that this particular integral does not have a straight-forward analytic formula. I checked Wolfram Alpha which also does not find a symbolic result (only a numeric one):
https://www.wolframalpha.com/input/?i=Integral%282*r*exp%28I*r**2%29*besselj%280%2C+3*r%29%2C+%28r%2C+0%2C+1%29%29
The impossibility of finding analytic results from integrals is a basic fact of symbolics: there do not exist recognisable mathematical functions to express the results of most possible integrals. Even in many cases where the functions do exist there do not exist algorithms that can find them.
So we accept that symbolic integration is impossible which means that there is no point in calling doit which is by far the slowest operation in the code that you are using. The purpose of doit would be to get a symbolic result which we could then use more efficiently for many different values of u and v but instead you are calling doit for every possible value of u and v and every time it does not return something useful.
If symbolic integration is not possible then numeric integration is needed. For numeric integration we first need to substitute values for u and v and then we can use a numeric integration algorithm. Since we don't have a symbolic result for the integral the numeric integration algorithm will need to be run once for each possible value of the parameters. There are two straight-forward ways to do this from SymPy: .evalf or lambdify (note that N(obj) is the same as obj.evalf() and also obj.n()).
Using .evalf is slower because it uses multiprecision arithmetic but that also makes it more accurate. Every digit in the returned result can be reasonably be trusted as correct (unless there'a a bug...) and you can also compute more digits if desired:
In [64]: %time h2.subs({uu:2,vv:3}).evalf()
CPU times: user 1.35 s, sys: 0 ns, total: 1.35 s
Wall time: 1.35 s
Out[64]: 0.236305280388988 + 0.0146270775320418⋅ⅈ
In [65]: %time h2.subs({uu:2,vv:3}).evalf(50)
CPU times: user 3.64 s, sys: 8 ms, total: 3.65 s
Wall time: 3.65 s
Out[65]:
0.23630528038898808262356997039156552180149071451802 + 0.014627077532041813385103750863443465733130
09416169⋅ⅈ
This is faster than calling doit but is still slow compared top standard machine precision numerical code. It is very accurate though.
For faster evaluation you can use lambdify:
In [66]: f2 = lambdify((uu, vv), h2)
In [67]: %time f2(2, 3)
/home/oscar/current/sympy/sympy.git/38venv/lib/python3.8/site-packages/scipy/integrate/quadpack.py:463: ComplexWarning: Casting complex values to real discards the imaginary part
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 438 µs
Out[67]: 0.2363052803889881
This now returns very quickly and gives an accurate result except that it only returns the real part. I guess that the scipy routine only handles real integrals. We can instead do the real and imaginary parts separately e.g. the imaginary part is:
In [91]: %time f2(2, 3)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 387 µs
Out[91]: 0.01462707753204179
