# Numerical Integration convergence

u0[p_, al_] :=
(1 - E^((-al)*p))/1;
K1[y_, al_] :=
E^(al*(-1 + y))*
(Gamma[0, -al] -
Gamma[0, al*(-1 +
y)]) - Log[1 - y];
Plot[
{NIntegrate[
((1 - E^((-al)*(1 -
y)))*(1 -
E^((-al)*y)))/
(y*(1 - y)*(1 -
E^(-al))),
{y, 0, 1}]*NIntegrate[
(u0[y, al]/y)*
(1 - E^((-al)*(1 -
y))), {y, 0, 1},
MaxRecursion -> 30,
Method ->
al*NIntegrate[
((((1 - E^((-al)*
(y - y1)))*(1 -
E^((-al)*y1)))/
(al*y1*(1 - y1)))*
(1 - E^((-al)*
(1 - y)))*(2 - y))/
(y*(1 - y)),
{y, 1.*^-7,
0.99999999}, {y1, 0,
y}, MaxRecursion ->
30, Method ->
NIntegrate[
(1/(y*(1 - y)))*
(((E^(al*(1 - y)) - 1)/
(1 - E^al))*
K1[y, al] +
((E^((-al)*(1 - y)) -
1)/(1 - E^(-al)))*
K1[y, -al]),
{y, 0, 1},
MaxRecursion -> 30,
Method ->
{al, 100, 1000},
PlotRange -> All]


I get the warning that integrand is highly oscillating and because of it, the integral is not the smooth function it should be. Can you please suggest can I trust the result obtained from this integral?

• Use the option Method -> "LocalAdaptive" of NIntegrate Commented Apr 2, 2021 at 8:00
• The method specification does not solve the problem, it just don't show the warning. I edited the question, the output I get is not smooth which it should be, I suspect the possible convergence issue of the NIntegerate. Commented Apr 2, 2021 at 9:43

All of the integrals can be integrated symbolically in one dimension. The one-dimensional integrals can be integrated completely and the two-dimensional integrals can be reduce to one-dimensional numerical integrals.

u0[p_, al_] := (1 - E^((-al)*p))/1;
K1[y_, al_] :=
E^(al*(-1 + y))*(Gamma[0, -al] - Gamma[0, al*(-1 + y)]) - Log[1 - y];

nIntegrate[args___] := Block[{NIntegrate},
Integrate[args] /.
HoldPattern@Integrate[f_, y__List, opts___?OptionQ] :>
NIntegrate[f, y]
];

Block[{NIntegrate},
i1[al_] =
nIntegrate[((1 - E^((-al)*(1 - y)))*(1 - E^((-al)*y)))/(y*(1 -
y)*(1 - E^(-al))), {y, 0, 1}, Assumptions -> al > 1];
i2[al_] =
nIntegrate[(u0[y, al]/y)*(1 - E^((-al)*(1 - y))), {y, 0, 1},
Assumptions -> al > 1];
i3[al_] =
al*nIntegrate[((((1 - E^((-al)*(y - y1)))*(1 - E^((-al)*y1)))/(al*
y1*(1 - y1)))*(1 - E^((-al)*(1 - y)))*(2 - y))/(y*(1 -
y)), {y, 0, 1}, {y1, 0, y}, Assumptions -> al > 1];
i4[al_] =
nIntegrate[(1/(y*(1 - y)))*(((E^(al*(1 - y)) - 1)/(1 - E^al))*
K1[y, al] + ((E^((-al)*(1 - y)) - 1)/(1 - E^(-al)))*
K1[y, -al]), {y, 0, 1}, Assumptions -> al > 1];
]

opts = Options@NIntegrate;
SetOptions[
NIntegrate, {MaxRecursion -> 20, WorkingPrecision -> 20,
PrecisionGoal -> 6, AccuracyGoal -> 20}];
Block[{\$MaxExtraPrecision = 500},
Plot[{i1[al]*i2[al] - i3[al] - i4[al]}, {al, 100, 1000},
PlotRange -> All, WorkingPrecision -> 25]
]
SetOptions[NIntegrate, opts];
(* warnings are emitted as Plot tests the integrand at machine precision *)

• Many thanks to you, Michael E2. Commented Apr 3, 2021 at 2:57
• I get the warnings that the precision of the argument function is less than the Working precision (20). Can I ignore the warning and trust the data? Commented Apr 3, 2021 at 3:08
• @SachinKaushik Yes. I mentioned them in a comment at the end of the code. It's something that Plot does, plug in MachinePrecision values to test the integrand, even when you specify WorkingPrecision -> 25. You could plot it this way instead: ListLinePlot[Table[{al, i1[al]*i2[al] - i3[al] - i4[al]}, {al, 100, 1000, 20}]] Commented Apr 3, 2021 at 3:29
• Sir if i use log log plot whether with list plot or plot, I get the previous problem back. Commented Apr 7, 2021 at 11:27

Update 1: I made a mistake in my function sum of original post: I had coded (f1f2)-f3-f3 when it should have been (f1f2)-f3-f4. Made the correction below and still plotting a discontinuity. Please feel free to remove upvote or I'll remove it in a bit.

Update 2: Since f3[x] is causing the problem I tried separating this double integral into two integral functions and this at least initially, seems to remove the discontinuity around a=649. Try to incorporate this into the original problem and see if it resolves the problem.

g[y_?NumericQ, a_?NumericQ] :=
NIntegrate[((((1 - E^((-a)*(y - x)))*(1 - E^((-a)*x)))/(a*
x*(1 - x)))*(1 - E^((-a)*(1 - y)))*(2 - y))/(y*(1 - y)), {x,
0, y}];
h[a_?NumericQ] := NIntegrate[g[y, a], {y, 1/10000, 9999/10000}]
Plot[{h[a]}, {a, 649, 650}]


I tried separating each integration as integral functions f1, f2, f3, and f4 and then plotting (f1*f2)-f3-f4:

f1[a_?NumericQ] :=
NIntegrate[((1 - E^((-a)*(1 - y)))*(1 - E^((-a)*y)))/(y*(1 - y)*(1 -
E^(-a))), {y, 0, 1}, Method -> "LocalAdaptive"];

f2[a_?NumericQ] :=
NIntegrate[(u0[y, a]/y)*(1 - E^((-a)*(1 - y))), {y, 0, 1},
MaxRecursion -> 30, Method -> "LocalAdaptive"];

f3[a_?NumericQ] :=
a*NIntegrate[((((1 - E^((-a)*(y - y1)))*(1 - E^((-a)*y1)))/(a*
y1*(1 - y1)))*(1 - E^((-a)*(1 - y)))*(2 - y))/(y*(1 -
y)), {y, 1.*^-7, 0.99999999}, {y1, 0, y}, MaxRecursion -> 30,

f4[a_?NumericQ] :=
NIntegrate[(1/(y*(1 - y)))*(((E^(a*(1 - y)) - 1)/(1 - E^a))*
K1[y, a] + ((E^((-a)*(1 - y)) - 1)/(1 - E^(-a)))*K1[y, -a]), {y,
0, 1}, MaxRecursion -> 30, Method -> "LocalAdaptive"]

Plot[(f1[x] f2[x]) - f3[x] - f4[x], {x, 100, 1000}]


• Could you please suggest how to resolve this because if I plot these terms separately, they are smooth functions so this discontinuity does not make sense. Is there any way I can correct this? Commented Apr 2, 2021 at 10:48
• @MichaelE2 Could you look into this problem. Commented Apr 2, 2021 at 10:58
• The functions are not smooth: f3[x] above has a small discontinuity around t=650 which you can see if you plot it separately. Perhaps focus on this function. Commented Apr 2, 2021 at 11:07
• Could you please check and suggest whether the nested integral definition is written properly or not in f3 function? and also can working precision or precision goal is helpful in this particular case? Commented Apr 2, 2021 at 11:21
• see update 2 above. I did this quick. You will need to check if the integration order is ok. Commented Apr 2, 2021 at 12:33