Automatically Switch between Integrate and NIntegrate

Question

I have an algorithm for which I need to compute some multidimensional integrals. Sometimes, Integrate works but some others Integrate just fails to compute the integral.

A simple example:

$$ \int^{\infty}_{-\infty}\int^{\infty}_{-\infty} -\frac{2 x^2 \left(e^{2 y+1}+1\right) (\text{erf}(y)+\text{erf}(y+1)) e^{-x^2-(y+1)^2}}{\pi } {\rm d}x {\rm d}y$$

Calling integrate,

Integrate[-((2 E^(-x^2 - (1 + y)^2) (1 + E^(1 + 2 y)) x^2 (Erf[y] + Erf[1 + y]))/\[Pi]), {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]

yields

$$ \int_{-\infty }^{\infty } -\frac{e^{-(y+1)^2} \left(e^{2 y+1}+1\right) (\text{erf}(y)+\text{erf}(y+1))}{\sqrt{\pi }} \, {\rm d}y $$

which is 0 (Interestingly, WolframAlpha indeed handles this integral, while my Mathematica does not). Now, Calling NIntegrate,

NIntegrate[-((2 E^(-x^2 - (1 + y)^2) (1 + E^(1 + 2 y)) x^2 (Erf[y] + 
 Erf[1 + y]))/\[Pi]), {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, AccuracyGoal -> 14]

yields $-5.2\times 10^{-16}$ (After a bunch of warnings). I was wondering I there was any way of further calling NIntegrate whenever Integrate fails. I know that I can call

N[Integrate[...]]

but I have not figured out how to pass method options to this. I am looking for a generic solution that works for higher dimensional integrals. If it helps, my integrals all involve exponentials and error functions similar to that one.

Thanks in advance.

Answer 1

This can be achieved with quite a simple function. Let's define a few integrands to test on

integrand1 = -((2 E^(-x^2 - (1 + y)^2) (1 + 
         E^(1 + 2 y)) x^2 (Erf[y] + Erf[1 + y]))/π);
integrand2 = Exp[-x^2 - y^2];
integrand3 = Sin[Sin[Sin[x]]]/x;

Now define the integrator you requested. It checks whether integration has succeeded by seeing whether the result includes the symbol Integrate

int[u__] := Module[{ans},
  ans = Integrate[u];
  If[FreeQ[ans, Integrate], ans, N[ans]]]

Now test them

int[integrand1, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
(* 0 *)

[It looks like Mathematica V12 computes your integral is analytically]

Here's another case that is possible analytically

int[integrand2, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
(* π *)

and one that isn't

int[integrand3, {x, -2, 2}]
(* 2.75358 *)

------------------WITH OPTIONS-------------

Here is an alternative that (I hope) supports options for both Integrate and NIntegrate routing them appropriately.

int[u__, v___Rule] := Module[{ans},
  ans = Integrate[u, Sequence @@ FilterRules[{v}, Options[Integrate]]];
  If[FreeQ[ans, Integrate], ans, 
   NIntegrate[u, 
    Evaluate[Sequence @@ FilterRules[{v}, Options[NIntegrate]]]]]]

There are a lot of cases to test to prove that this works fully, so here are a few:

int[integrand1, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
(* 0 *)

int[integrand2, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
(* π *)

int[integrand3, {x, -2, 2}]
(* 2.75358 *)

int[integrand3, {x, -2, 2}, Method -> MonteCarlo]
(* 2.74324 *)

int[1/x, {x, -3, 3}, PrincipalValue -> True]
(* 0 *)