I'm doing a numerical integration that involves an integrand with a few exponentials whose values over the integration region range enormously. I'm really only interested in the part of the range that actually contributes to the integral, but the other parts are giving me errors relating to precision and extremely small values. I've read this and this but I'm still having trouble putting the concepts into action.

Here is an example that's contrived but similar to my actual application:

fn1 = 1/(1 + Exp[500*(x - 2)]);
Plot[fn1, {x, 0, 5}]
LogPlot[fn1, {x, 0, 5}]

Here it is, in the range I'm looking at:

enter image description here

It looks like there's a discontinuity, but it's actually smooth if you zoom in:

enter image description here

This varies enormously over the integration range of {0,5}. At x = 1, it's about 1, but at x = 3, it's about 10^-218. That's so small, I don't really need it in my answer. But obviously, where you decide to cut off the number is your choice; depending on the application, I could decide that 0.01 is 'small' compared to the values of the integrand over most of the range.

So let's say that what I want to do is calculate my integral to a factor of 10^9 below the maximum value in the integration range. How can I adjust the Precision of my integrand, the WorkingPrecision, and the AccuracyGoal so that it can be done automatically without throwing errors or making me manually adjust it every time?

Here's what I tried. I thought that if I look at the maximum value in my integration range, I can use that to 'tune' the accuracy I'd like, and thus the WP.

mv = NMaxValue[{fn1, x <= 5, x >= 0}, x]
ag = 9;
wp = ag*2;
NIntegrate[SetPrecision[fn1, wp], {x, 0, 5}, WorkingPrecision -> wp, 
 AccuracyGoal -> ag]

I thought that setting my AccuracyGoal to 9 would give me the desired accuracy I wanted, and setting the WorkingPrecision to twice that would give it enough workingprecision to do that, but it gives me this error:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

We know it doesn't have a singularity, the integration isn't equal to 0, and it's not oscillating. So, the WP must be too small. If I instead set wp=ag*3, it doesn't throw the error. So how can I figure out the WP necessary, in advance, to get the Accuracy/Precision I want, without it throwing an error?

EDIT: Someone suggested I show my actual application in case a given solution that applies to my simple example isn't general enough for mine.

Here is a pastebin of the awful massive function that's my actual integrand.

For some clarity, the integrand is actually the product of two functions, plotted here:

enter image description here

The functions quickly go to absurdly small numbers away from the small section they barely overlap at. Here is a plot of their product, the integrand:

enter image description here

The peak of the region that actually contributes most to the integral is about 6*10^-4, which isn't really that small of a number. However, if you look at a region away from that peak, for example the value at Ex = Er = 4, it's about 10^-8400... which we probably don't need to bother with.

So how can I make my function do the following: Find the max value of my integrand in the integration region, divide that number by something conservative like 10^9, and then don't bother integrating over any values less than that? Essentially, if a value is less than that, take it to be 0.

edit2: It seems like Chop[] could potentially work for what I want to do. I could do something like

mv = NMaxValue[{integrand, Ex <= upperlim, Ex >= 0, Er <= upperlim, 
    Er >= 0}, {Ex, Er}];
  Chop[integrand, mv/10^9], {Ex, 0, upperlim}, {Er, 0, upperlim}];

However, it seems to be throwing a lot more errors than before, so I suspect I'm using it wrong...

  • $\begingroup$ NDSolveValue[{y'[x] == fn1, y[0] == 0}, y[5], {x, 0, 5}, WorkingPrecision -> 20]? -- Also, I feel your MWE might be too simple. Some solvers will work quickly on it, but might not on the actual use case. E.g. NIntegrate[fn1, {x, 0.`, 0.000021073377463201208`, 0.3161217353254813`, 1.8161217353254813`, 2.1854546039338802`, 5.`}]] has an error just over 10^-12, but I'm not sure whether the approach that produced it is general enough. Can you say a little more, if the NDSolveValue trick doesn't work? $\endgroup$
    – Michael E2
    Commented Mar 15, 2016 at 3:36
  • $\begingroup$ @MichaelE2, thanks for the advice. It's actually a 2 variable function, so I'm not sure how to apply your method (I'm sure it can be done, I just don't know how). I'll edit my OP to show the actual example, thanks! $\endgroup$ Commented Mar 15, 2016 at 15:27
  • $\begingroup$ I think i made the same comment on another question of yours. If you require high precision then all the numerical values in your integrand need to be suitably high precision. If you wrote your integrand in terms of (looks like) 4 or 5 symbols for the numeric values i think it would be small enough to post here as well. $\endgroup$
    – george2079
    Commented Mar 18, 2016 at 16:55

3 Answers 3

fn1[x_] = 1/(1 + Exp[500*(x - 2)]);

The integral can be evaluated exactly

int = Integrate[fn1[x], {x, 0, 5}]

(*  (1/500)*Log[1 + E^1000]  *)

int // N

(*  2.  *)

To see how close to 2 the integral is

Block[{$MaxExtraPrecision = 1000},
 N[Log10[N[int - 2, 500]], 6]]

(*  -436.993  *)

So it differs by approximately 10^-437

  • 2
    $\begingroup$ Thanks for the response, but I don't think this answers the question. It asked how to adjust the WP, and AG or PG so that it achieves the desired accuracy, in a general case. The example I gave was a very simplified version; the actual one definitely can't be solved analytically. $\endgroup$ Commented Mar 15, 2016 at 15:05

You probably need a LocalAdaptive Method.

Basically, the LocalAdaptive strategy is to recursivelly calculate the integral over smaller disjoint regions using an integration rule - the recursion stopping criteria is user-specified.

NIntegrate[fn1[x], {x, 0, 5}, AccuracyGoal -> ag, PrecisionGoal -> pg,
  WorkingPrecision -> wp, Method -> "LocalAdaptive"]

NIntegrate computes a initialIntegral over an initial region and then partitions it to disjoint regions. NIntegrate compares the regionError of the disjoint regions to the initialIntegral. If the error is significant, the recursion is called. The error will be insignificant if

initialIntegral + regionError = initialIntegral

The region's integral estimate is the sum of the integral estimates returned from these recursive calls. Since you want to compute the integral estimate to user-specified precision and accuracy goals, the following stopping criteria is used:

integralEst = Min[initialIntegral 10^-pg / eps, 10^-ag / eps]
integralEst + regionError == integralEst

A slightly different strategy is to integrate only over the region above some threshold value, which can be fraction of the maximum: mv/a.

r = Reduce[0 <= x <= 5 && fn1 > mv/a, x] (* 0 <= x < 2.05526 *)
NIntegrate[fn1, {x, r[[1]], r[[5]]}] (* 2. *)

I lowered the fraction, because with a=10^9 the result didn't seem reasonable (2.041446531671895 compared to 2.00000000000000090701555082 that results from your initial code)

I tried using Implicit Region, which would be a one-liner in theory, but couldn't get it to work.

NIntegrate[fn1, {x} \[Element] ImplicitRegion[fn1 > mv/a, {{x, 0, 5}}]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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