3
$\begingroup$

I am working on a Mathematica project that involves many numeric integrals, each of which should be integrated over the domain {k, 0, Infinity}. I assumed that Mathematica would have an easier time integrating to a large finite number, so for all my integrals I used the domain {k, 0, kMax} so I could play with large values of kMax later.

The trouble is, each integral is a "step-like" function that drops off exponentially quickly after some cutoff value, though this cutoff is different for each integral. If I make kMax big enough to capture the important contribution from the integrals with the largest cutoffs, then the integrals with small cutoffs start to have issues.

As an example, observe that NIntegrate[E^(-x^2), {x, 0, 1000}] yields the correct value, whereas NIntegrate[E^(-x^2), {x, 0, 10000}] yields zero and throws an error.

What are some standard solutions to this issue? I've considered several approaches, but all have their setbacks:

  1. Use {k, 0, Infinity} for each integral. Perhaps Mathematica is smarter about finding cutoff values in this case. The downside is I don't really know what Mathematica is doing.
  2. Use "smart" cutoffs for each integral, that I have to figure out and program in. The downside is the messiness of this approach, and the assumptions I'd have to make about where the "right" cutoff location is for a good approximation.
  3. Use one large cutoff kMax value and also a very MinRecursion, MaxRecursion, and precision/accuracy settings. The downside is huge computation times, perhaps becoming unfeasible.
Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

Typically, one maps the infinite interval to a finite one by a rational function of the form $$x \mapsto \frac{a \, x+ b}{c \, x + d}$$ and employs the substitution rule. Then you do not need any cutoff (provided that the integrand is integrable in the first place).

This will probably be what Mathematica does for a symbolic input. But if you do that once, you can tell Mathematica to skip the symbolic processing step (Method -> {Automatic, "SymbolicProcessing" -> False}) and to just apply a numerical quadrature of your choice on a finite interval; typically Method->"GaussKronrod" works very well. For more details on numerical integration in Mathematica see, e.g., the Method section of NIntegrate's doc page or this documentation page: https://reference.wolfram.com/language/tutorial/NIntegrateIntroduction.html.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Interesting. Just to make sure I understand, if I have $\int_0^\infty f(x)\,dx$, you're saying to substitute $u=\frac{ax+b}{cx+d}$? Then we have $\int_{b/d}^{a/c}f\left(\frac{-ud+b}{cu-a}\right)\frac{ad-bc}{(cu-a)^2}\,du$. Presumably the singularity at the endpoint $a/c$ will not be an issue as long as the original integral was finite? $\endgroup$
    – WillG
    Commented Mar 3, 2020 at 18:40
  • $\begingroup$ Yes. At the infinite end, the function $f$ has to vanish (decay sufficiently quickly) so there should be no singularity. Btw., also the functions $\log$ and $\arctan$ might help to squeeze infinite intervals into a finite one. $\endgroup$
    – Henrik Schumacher
    Commented Mar 3, 2020 at 18:59
  • $\begingroup$ Unfortunately the new singularities at the endpoints, though integrable in theory, are causing new issues for Mathematica. It's throwing all sorts of errors about the integrals not converging to prescribed accuracy and WorkingPrecision being too small. Simple fiddling with these parameters doesn't seem to resolve it. $\endgroup$
    – WillG
    Commented Mar 3, 2020 at 21:21

Not the answer you're looking for? Browse other questions tagged or ask your own question.