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:
It looks like there's a discontinuity, but it's actually smooth if you zoom in:
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:
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:
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}];
NIntegrate[
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...
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 over10^-12
, but I'm not sure whether the approach that produced it is general enough. Can you say a little more, if theNDSolveValue
trick doesn't work? $\endgroup$