Mathematica Asked by John Taylor on January 12, 2021
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.
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}]
Correct answer by Wen Chern on January 12, 2021
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 by David G. Stork on January 12, 2021
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 by bbgodfrey on January 12, 2021
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}
Answered by Carl Woll on January 12, 2021
A few comments to the posted comments and answers follow.
@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} *)
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.33227*10^-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 *)
Answered by Anton Antonov on January 12, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP