12
$\begingroup$

Sometimes when I do integrals in Mathematica (M), it keeps thinking and thinking and I have no idea what is going on inside M. For how long should one wait or how does one know whether M has not got stuck in some loop or something? These integral may be algebraic or numerical ones. Is there a rule of thumb for how much time one should give M to evaluate an integral before aborting?

$\endgroup$
  • 1
    $\begingroup$ Have you seen TimeConstrained? $\endgroup$ – J. M. will be back soon Jun 8 '15 at 13:57
  • $\begingroup$ No I have not, but what I meant was that if there is a way to know whether one should wait some more or simply abort? I.e. if M is really working or just simply stuck. $\endgroup$ – Your Majesty Jun 8 '15 at 14:11
  • $\begingroup$ It can be hard to tell, admittedly; hence my suggestion. $\endgroup$ – J. M. will be back soon Jun 8 '15 at 14:12
  • 2
    $\begingroup$ The numeric and analytic cases are very different. NIntegrate might reasonably take a very long time just because the integrand evaluation takes a long time. As to analytic Integrate, I cant say I've seen it take more than an hour or so and come back with a useful result. If its important maybe let it chug away over night, but no longer than that.. $\endgroup$ – george2079 Jun 8 '15 at 15:16
17
$\begingroup$

For what it's worth, you can monitor the progress of Integrate and NIntegrate. I'm not sure how helpful the tools below are, but I feel it is probably worth mentioning them.

Integrate

Internal`Integrate`debugSwitch

If you set Internal`Integrate`debugSwitch to the magic number 10, it will print its (major) steps in searching for the answer. For instance:

Block[{Internal`Integrate`debugSwitch = 10},
 Integrate[Exp[-x]/Sqrt[1 + x^3], {x, 0, Infinity}]
 ]

Be prepared for lots of output in terms of internally named functions. It's rather hard to understand what's going on, though not completely impossible. For a computation that goes on for a long time, it would be a lot of stuff to wade through.

Other debug printers

Also, each submethod of Integrate seems to have its own debug-printing functions. Here are some (I don't know them all). Each has to be loaded, since most of the integration code is loaded as needed. This can be done by calling Integrate.

Integrate[Exp[-x]/Sqrt[1 + x^3], x];
Integrate[Exp[-x]/Sqrt[1 + x^3], {x, 0, Infinity}];
Names["*Int*`*db*"]
(*
 {"Integrate`EllipticTrig`db", "Integrate`ImproperDump`dbgPrint", \
  "Integrate`AlgebraicFunctionDump`dbgPrintalg", \
  "Integrate`ConvolutionGGDump`dbgPrintcongg", \
  "Integrate`NLtheoremDump`dbgPrintDI", \
  "Integrate`Elliptic`dbgPrintEE", \
  "Integrate`EllipticTrig`dbgPrintEET", \
  "Integrate`ExponentialsDump`dbgPrintexp", \
  "Integrate`ImproperDump`dbgPrintII", \
  "Integrate`LogarithmDump`dbgPrintlog", \
  "Integrate`QuickLookUpDump`dbgPrintQT", \
  "Integrate`FindIntegrandDump`dbgPrintslatr", \
  "Integrate`TableDump`dbgPrintT", "Integrate`TableDumpExp`dbgPrintTE", \
  "Integrate`TableDumpSpec`dbgPrintTS", \
  "Integrate`TableDumpTrig`dbgPrintTT", "Integrate`EllipticTrig`db$"}
*)

They can be used like this

Block[{Integrate`QuickLookUpDump`dbgPrintQT = Print},
 Integrate[x Sin[x], x]
 ]

or this

Reap@Block[{Integrate`QuickLookUpDump`dbgPrintQT = Sow[{##}] &},
  Integrate[x Sin[x], x]
  ]

Again, probably of limited usefulness.

NIntegrate

I don't know of any debug-Print-like controls within NIntegrate, but you do have EvaluationMonitor and the undocumented IntegrationMonitor.

EvaluationMonitor

For EvaluationMonitor, perhaps the best use is to count the number of evaluations, which can be displayed dynamically in real time while the computation is ongoing.

f[x_?NumericQ] :=  (* N.B. designed to be particularly SLOW *)
  Interpolation[Table[{t, Sin[100 t]}, {t, 0, 1, 1./2^10}], InterpolationOrder -> 1][x];
k = 0;
Dynamic[k]
NIntegrate[f[x], {x, 0, 1}, EvaluationMonitor :> k++, MaxRecursion -> 20]
(* k is 34279 at the end *)

One could also store the evaluation points. I thought something like this would work, since dynamic updating is designed to have a low priority, but when the number of evaluation points gets large, thing slow down. The kernel would actually crash at low settings of UpdateInterval. Bummer...and beware.

f[x_?NumericQ] := 
  Interpolation[Table[{t, Sin[100 t]}, {t, 0, 1, 1./2^10}], InterpolationOrder -> 1][x];

Clear[g];
k = 0;
pts = Sequence[];
g := Graphics[         (* Graphical monitor *)
   Dynamic@Refresh[
     Point[#[[;; ;; Ceiling[Length[#]/800]]] &@Level[pts, {-2}]], 
     TrackedSymbols :> {}, UpdateInterval -> 1],
   PlotRange -> All, AspectRatio -> 1/2, Frame -> True];
Dynamic@g

NIntegrate[f[x], {x, 0, 1}, 
 EvaluationMonitor :> (pts = {pts, {{++k, x}}}), MaxRecursion -> 20]

g = Graphics[Point[Level[pts, {-2}]],  (* set g to be static *)
   PlotRange -> All, AspectRatio -> 1/2, Frame -> True];

enter image description here

IntegrationMonitor

For IntegrationMonitor option (at least in the "GlobalAdaptive" strategy), its setting should be a function, and the function seems to be passed a data structure like this:

{NIntegrate`IntegrationRegion[{{0}, {1}}, 
  Experimental`NumericalFunction[{x}, 1/Sqrt[x], "-NumericalFunctionData-"], {}, 
  NIntegrate`GeneralRule[{
    {0.007957319952578756`, 0.046910077030668074`, 0.12291663671457531`,
     0.2307653449471585`, 0.3601847934191085`, 0.5`, 0.6398152065808915`, 
     0.7692346550528415`, 0.8770833632854247`, 0.9530899229693319`, 0.9920426800474212`},
    {0.0212910183755409`, 0.057616658311237225`, 0.0934003982782463`, 0.12052016961432388`,
     0.13642490095627946`, 0.1414937089287456`, 0.13642490095627946`, 0.12052016961432388`,
     0.0934003982782463`, 0.057616658311237225`, 0.0212910183755409`},
    {0.0212910183755409`, -0.0608467842168579`, 0.0934003982782463`, -0.11879416563535941`, 
     0.13642490095627946`, -0.14295073551569884`, 0.13642490095627946`,
    -0.11879416563535941`, 0.0934003982782463`, -0.0608467842168579`, 0.0212910183755409`}}
  ]]}

The first argument is the subinterval being integrated. One could monitor the progress with

NIntegrate[1/Sqrt[x], {x, 0, 1}, MaxRecursion -> 3, 
 IntegrationMonitor :> (Print[#[[1, 1]]] &)]
(*
  {{0},{1}}
  {{0},{0.5}}
  {{0},{0.25}}
  {{0},{0.125}}
*)

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 3 recursive bisections in x near {x} = {0.124005}. NIntegrate obtained 1.9779036090196027and 0.033906515833193304 for the integral and error estimates. >>

(*  1.9779  *)

This may vary with the method strategy and/or rule.

Notes

I'm not sure where I heard of debugSwitch. A search of the site turned up nothing (neither did Google). I'm pretty sure it was mentioned in a comment on this site, or maybe by spekunking. But you can see it in a Trace, e.g. Trace[Integrate[x Sin[x], {x, 0, 1}], TraceInternal -> True]. Lots of solvers have *Print methods built into them. They can be found with Names["*`*Print"], although, like Integrate the names might not show up until the corresponding module has been loaded;. There are several kinds of *Print methods, but the debug ones often have db, dbg, debug or Debug preceding Print and start with a lower-case letter. The output is not particularly user-friendly and I suspect they were not meant for end-users.

I discovered IntegrationMonitor via Trace: For example,

Trace[
 NIntegrate[1/Sqrt[x], {x, 0, 1}],
 _NIntegrate`InitializeIntegrationStrategy,
 TraceInternal -> True]

Searching the site, I found it has been used once so far, by Simon Woods in his answer to Is there a way to see the result of NIntegrate's symbolic preprocessing?

$\endgroup$
  • $\begingroup$ @LoveLearning You're welcome. $\endgroup$ – Michael E2 Jun 9 '15 at 4:03
  • 1
    $\begingroup$ @LoveLearning perhaps you could consider editing the question title to something like "How to monitor Integrate and NIntegrate progress?". That way the answerable part is emphasized and the other one ("How long does it take?") is kept only in the inner text. The main purpose of the proposed change would be to insert useful keywords (monitor, integrate, progress) in the question title for future searchers $\endgroup$ – Dr. belisarius Jun 10 '15 at 15:38
8
$\begingroup$

There probably cannot be general answer for such a general question because the bigger integrand (other conditions being equal) the more time Integrate needs to handle it. But for basic tabular integrands there are some benchmarks made by Albert Rich, the developer of Rubi - rule-based integrator. On the linked page a table is given where the "Timeout" column means "the number of problems the system fails to integrate within a 25 second timelimit". From this table it is clear that Mathematica 9 on his computer (characteristics are absent but I'm quite sure it is an ordinary cheap computer, not a super-computer) cannot finish in 25 seconds evaluation of only 2% of the integrals from his test suite (which can be downloaded from the same page). So you can expect that Mathematica is able to take almost any standard indefinite integral within 25 seconds. If it fails then probably it is just unable to take it at all.

Also you can download his test suite and produce your own benchmarks.

$\endgroup$
  • $\begingroup$ Cool answer. So in short, one should use other software to numerically integrate? I used to use Matlab for that and liked it, but thought that I should learn that in Mathematica too. RUBI seems awesome btw. $\endgroup$ – Your Majesty Jun 8 '15 at 19:06
  • 3
    $\begingroup$ @Love, not necessarily; note that NIntegrate[] supports a number of methods that can be tried if the default setting is not up to snuff. But more importantly, the nature of your integrand (and its accompanying singularities) will often demand the use of a non-default setting. $\endgroup$ – J. M. will be back soon Jun 8 '15 at 19:37
  • 2
    $\begingroup$ this page has CPU time of problems used, for V 10.1. No time out is used in this one. $\endgroup$ – Nasser Jun 8 '15 at 22:17
  • $\begingroup$ Cool interesting @Nasser $\endgroup$ – Your Majesty Jun 8 '15 at 22:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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