I have the following integral (inspired by this question) involving two real variables a and b (and the imaginary unit I):
f = integrate(exp((-a*t**2+I*b*t)/(3*t**2+1)+I*t*x)*(x/(3*t**2+1)), (t, -oo, oo), (x, 0, oo))
I would like to iterate both variables a and b, let us say in the range from -10 to 10, e.g. for a in range(-10,11): followed by an inner loop for b in range(-10,11): and then evaluate the integral, which allows me to draw a surface depending from both variables a and b.
The approach looks as follows:
from sympy.abc import a, b, x, t
from sympy import *
from sympy import var, Integral, init_printing
init_printing()
a, b, x, t = symbols('a b x t')
integrand = exp((-a*t**2+I*b*t)/(3*t**2+1)+I*t*x)*(x/(3*t**2+1))
for i in range(-1,1,0.1):
for j in range(-1,1,0.1):
print(Integral(integrand.subs(a,i).subs(b,j), (t, -oo, oo), (x, 0, oo)).evalf())
I tried to evaluate it using:
Integral(exp((-a*t**2+I*b*t)/(3*t**2+1)+I*t*x)*(x/(3*t**2+1)), (t, -oo, oo), (x, 0, oo)).evalf()
but without success (syntactically it seems to be correct). It would be great to know how we can approach such evaluations. I am open to any alternative approach. I am really curious to know what the surface F(a,b) looks like.
Under the hood when you use evalf for an integral SymPy will use the corresponding routines from the mpmath library:
https://mpmath.org/doc/current/calculus/integration.html
Although mpmath supports numerically evaluating double integrals SymPy will only attempt this for single integrals. You can use mpmath directly though:
In [46]: import mpmath
In [47]: a = 1
In [48]: b = 1
In [49]: I = mpmath.sqrt(-1)
In [50]: f = lambda x, t: mpmath.exp((-a*t**2+I*b*t)/(3*t**2+1)+I*t*x)*(x/(3*t**2+1))
In [51]: mpmath.quad(f, [0, mpmath.inf], [-mpmath.inf, mpmath.inf])
Out[51]: mpc(real='-1.1728243698636993e+37', imag='0.0')
I think that -1e37 result probably just means that the integral doesn't converge (or at least that the numerical algorithm trying to compute it doesn't converge). Picking finite values for the upper limit of the outer integral suggests that this does not converge:
In [54]: for c in [1, 10, 100, 1000, 10000]:
...: print(mpmath.quad(f, [0, c], [-mpmath.inf, mpmath.inf]))
...:
(0.496329567173036 + 0.0j)
(3.47061267913388 + 0.0j)
(6.20670960862542 + 0.0j)
(-245.035824572364 + 0.0j)
(149439.86145905 + 0.0j)
Maybe it converges for some values of a and b but not others.