# Problem with running this code in Wolfram Mathematica

I have a hard time running this:

Integrate[
Exp[-q^2 - (v^2 r^2/16/a^2/q^2) - (u^2 r^2/16/a^2/q^2)], {q,
r/(4 a t)^(1/2), Infinity},
Assumptions -> a >= 0 && u >= 0 && t >= 0 && r >= 0 && v >= 0]


What should I do?

• always try the basic integral first (indefinite) and see if that works. Integrate[Exp[-q^2 - (v^2 r^2/16/a^2/q^2) - (u^2 r^2/16/a^2/q^2)], q] does not integrate (with or without assumptions). Either Mathematica does not know how, or this is a non integrable. You can try numerical integration. Oct 15 at 23:27
• btw, your integrand has the form $e^{-\frac{r^2 \left(u^2+v^2\right)}{16 a^2 q^2}-q^2}$ Oct 15 at 23:36
• Re "problem" and "I have a hard time": Can you be more specific? What is happening (or not happening)? What are the symptoms? What is the problem? Please respond by editing (changing) your question, not here in comments (without "Edit:", "Update:", or similar - the question should appear as if it was written right now). Oct 16 at 8:50

Your integrand has the basic form

$$e^{-\frac{r^2 \left(u^2+v^2\right)}{16 a^2 q^2}-q^2}$$

Which Mathematica can't integrate. But it can integrate it when $$\left(u^2+v^2\right)$$ is a single symbol, like this

$$e^{-\frac{B r^2}{16 a^2 q^2}-q^2}$$

So by doing this replacement, (since $$u,v$$ do not show anywhere else), it can now do it.

Clear["Global*"]
rep = (v^2 + u^2) -> b;
term = Simplify[-q^2 - (v^2 r^2/16/a^2/q^2) - (u^2 r^2/16/a^2/q^2)]
term = term /. rep
res=Integrate[Exp[term], {q, r/(4 a t)^(1/2), Infinity},
Assumptions -> a >= 0 && t >= 0 && r >= 0 , GenerateConditions -> False] Now you can simply replace b back with u^2+v^2

res /. b -> (v^2 + u^2) To verify the above, we can do the indefinite integration and differentiate the antiderivative to see if we get the original integrand back

rep = (v^2 + u^2) -> b;
term = Simplify[-q^2 - (v^2 r^2/16/a^2/q^2) - (u^2 r^2/16/a^2/q^2)];
OriginalIntegrand = Exp[term];
term = term /. rep;
res = Integrate[Exp[term], q,
Assumptions -> a >= 0 && t >= 0 && r >= 0,
GenerateConditions -> False];
res = res /. b -> (v^2 + u^2);
backIntegrand = D[res, q] // Simplify;
backIntegrand == OriginalIntegrand
` 