Integral of the integral using NIntegrate

Question

Consider the integral

NIntegrate[Exp[-NIntegrate[x, {x, 0, y}]], {y, 0, 2.8}]

It displays an error x = y is not a valid limit of integration, however, gives some number. What is a reason for the error and how to interpret the numeric result?

Of course, I can easily replace the integration symbol in the exponent by the value of the integral, however, this is only a toy example demonstrating the problem.

Of course you cannot perform the inner numerical integral with an un-specified (free) variable. — David G. Stork, Jul 26 '19 at 22:39
@DavidG.Stork : however, it turned out that Mathematica gives the correct answer. Why this is possible? — John Taylor, Jul 26 '19 at 22:40
Hmmm.... a tricky question about the internal methods of Mathematica. Frankly, I don't know why the correct answer is given (despite the error messages). — David G. Stork, Jul 26 '19 at 22:43
Mathematica does not have to evaluate the integrand symbolically, in order to apply numerical quadrature. However, it does so usually in order to analyze the integrand. If it fails, it throws an error and tries purely numerical methods. This why it works in the end. — Henrik Schumacher, Jul 26 '19 at 22:46
@HenrikSchumacher I'm wondering if I turn off SymbolicProcessing, would this error still present? — Turgon, Jul 27 '19 at 08:34
@Turgon How would you do that? SymbolicProcessing is not an option of NIntegrate. — Henrik Schumacher, Jul 27 '19 at 12:59
@HenrikSchumacher actually it is an option of NIntegrate. Check the help page, it's under Options->Method->SymbolicProcessing. — Turgon, Jul 27 '19 at 13:41
@Turgo Ah, you mean Method -> {Automatic, "SymbolicProcessing" -> False}. Yes. Unfortunately, that does not remove the error messages. — Henrik Schumacher, Jul 27 '19 at 13:57
Using Method -> {Automatic, "SymbolicProcessing" -> False} indeed does not eliminate the error messages, but does reduce the number of them to two when applied to the outer integration. It seems to have no effect on the inner integration. — bbgodfrey, Jul 27 '19 at 18:32
Possible duplicate: https://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users/26037#26037, esp. ref. 6) under NIntegrate. — Michael E2, Jul 27 '19 at 20:13

score 8 · Accepted Answer · answered Jul 26 '19 at 22:44

8

You may try to restrict the definition of the integrand to avoid symbolic calculation. Like:


f[y_?NumericQ] := Exp[-NIntegrate[x, {x, 0, y}]];
NIntegrate[f[y], {y, 0, 2.8}]

answered Jul 26 '19 at 22:44

Wen Chern

716
4
8

score 5 · Answer 2 · answered Jul 27 '19 at 19:43

Here's a refinement of bbgodfrey's approach, where I combine both integrals into a single NDSolveValue call:

sol = NDSolveValue[
    {
    z'[x] == x, z[0]==0,
    int'[x] == Exp[-z[x]], int[0]==0
    },
    int,
    {x, 0, 2.8}
];
sol[2.8] //AbsoluteTiming

{7.*10^-6, 1.24691}

The nice thing about this approach is that computing the integral for various y limits is quick. Compare this timing to that of the other numerically based answers:

f[y_?NumericQ] := Exp[-NIntegrate[x, {x, 0, y}]];
NIntegrate[f[y], {y, 0, 2.8}] //AbsoluteTiming

{0.072347, 1.24691}

NDSolveValue[{z'[x] == x, z[0] == 0}, z[x], {x, 0, 2.8}];
NIntegrate[Exp[-%], {x, 0, 2.8}] //AbsoluteTiming

{0.026102, 1.24691}

Your approach indeed is much faster than mine (about a factor of eight, +1), but not by as much as your answer suggests. AbsoluteTiming as applied in your answer determines only the time required to evaluate sol[2.8], but omits the time required to compute the interpolation function sol itself. — bbgodfrey, Jul 27 '19 at 20:18

score 4 · Answer 3 · answered Jul 27 '19 at 18:25

4

A typically more efficient approach for complicated integrands is

NDSolveValue[{z'[x] == x, z[0] == 0}, z[x], {x, 0, 2.8}];
NIntegrate[Exp[-%], {x, 0, 2.8}]
(* 1.24691 *)

The inner integration is performed only once here.

answered Jul 27 '19 at 18:25

bbgodfrey

61,439
17
89
156

score 3 · Answer 4 · answered Jul 26 '19 at 22:41

3

Of course you cannot perform the inner numerical integral with an un-specified (free) variable, $y$.

Why not simply perform the analytic integrals and be done with it?

Integrate[Exp[-Integrate[x, {x, 0, y}]], {y, 0, 2.8}]

$1.24691$

answered Jul 26 '19 at 22:41

David G. Stork

41,180
3
34
96

2

I have used this simple integral as an example. Realistic integrals I want to use are much more complicated. – John Taylor Jul 26 '19 at 22:46
1

@JohnTaylor: How do you know your "more complicated" integrals cannot be computed analytically? – David G. Stork Jul 28 '19 at 06:05

score 3 · Answer 5 · edited Jun 16 '20 at 09:23

A few comments to the posted comments and answers follow.

Faster computations without messages

@Turgo Ah, you mean Method -> {Automatic, "SymbolicProcessing" -> False}. Yes. Unfortunately, that does not remove the error messages. – Henrik Schumacher Jul 27 at 13:57

Using Method -> {Automatic, "SymbolicProcessing" -> False} indeed does not eliminate the error messages, but does reduce the number of them to two when applied to the outer integration. It seems to have no effect on the inner integration. – bbgodfrey Jul 27 at 18:32

Here is a definition that produces faster results without messages. The setting Method -> {Automatic, "SymbolicProcessing" -> 0} is used for both integrals.

Clear[f];
f[y_?NumericQ] := 
  Exp[-NIntegrate[x, {x, 0, y}, Method -> {Automatic, "SymbolicProcessing" -> 0}]];
AbsoluteTiming[
 res = NIntegrate[f[y], {y, 0, 2.8}, Method -> {Automatic, "SymbolicProcessing" -> 0}]
 ]
(* {0.00553, 1.24691} *)

Precision of the results

OP says the computation in his post is a simplified example of a computation with more complicated integrands. Depending on the integrands and precision goal requirements the NDSolve approach might be not precise enough.

Illustrating computations follow. (I show WorkingPrecision->30 and PrecisionGoal->20, but similar results are obtained with machine precision and PrecisionGoal->12.)

`Integrate`

F[y0_] := Integrate[Exp[-Integrate[x, {x, 0, y}]], {y, 0, y0}]
F[y0] 
(* Sqrt[\[Pi]/2] Erf[y0/Sqrt[2]] *)
F[2.8]
(* 1.24691 *)

`NDSolve`:

sol2 = NDSolveValue[{z'[x] == x, z[0] == 0, int'[x] == Exp[-z[x]], 
    int[0] == 0}, int, {x, 0, 28/10}, WorkingPrecision -> 30, 
   PrecisionGoal -> 20];
sol2[28/10]
(* 1.24690937538389563405549789088 *)
Abs[sol2[2.8] - F[28/10]]
(* 1.3322710^-14 )

`NIntegrate`

f[y_?NumericQ] := 
  Exp[-NIntegrate[x, {x, 0, y}, 
     Method -> {Automatic, "SymbolicProcessing" -> 0}, 
     WorkingPrecision -> 30, PrecisionGoal -> 20]];
res = NIntegrate[f[y], {y, 0, 28/10}, 
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}, 
  MaxRecursion -> 20, WorkingPrecision -> 30, PrecisionGoal -> 20]
(* 1.24690937538388234380595431698 *)
Abs[res - F[28/10]]
(* 0.10^-30 )

Integral of the integral using NIntegrate

5 Answers5

Faster computations without messages

Precision of the results

Integrate

NDSolve:

NIntegrate

`Integrate`

`NDSolve`:

`NIntegrate`