0
$\begingroup$

I have this integrand given by tinteg[z,zm,zh] where I want to integrate it in the range [0,zh]. However, the integrand is divergent at zh so I added a small parameter eps away from zh so that the integral range is now [0,zh+eps]. Doing this, I can then take the principal value of the integral for some eps.

d = 3;
ag = 20;
pg = 20;
wp = 40;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
tinteg[z_, zm_, zh_] := -1/(f[z, zh] Sqrt[1 - (zm^(2 d) f[z, zh])/(z^(2 d) f[zm, zh])])

th[eps_?NumericQ, zh_?NumericQ] := Module[{epsr, zhr}, {epsr, zhr} = Rationalize[{eps, zh}, 0]; NIntegrate[tinteg[z, zhr + epsr, zhr], {z, 0, zhr, zhr + epsr}, Method -> PrincipalValue, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp]]

th[10^-9, 100]
0.0003162277660196573967272197281906499

th[10^-10, 100]
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {99.9999999998213805278900613760597337931183240142657527684191416030996934771392010899316650}. NIntegrate obtained -57.3106917390753700387564367567637817259282112119845402535003746786552197747577990043682551 and 1.85452812514993194537323109011792300865115544959913889998941281161057638475694312252915879`90.*^-10 for the integral and error estimates.
0.0001000000000000889405834231780766118

th[10^-15, 100]
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {99.9999999999997866007398081727804784995174546220786282909845407978056563093393638396629520}. NIntegrate obtained -57.3107914228724344179090132179535225292646529296309563458710154094184035242065488421646708 and 0.000096232917814206821120246506287335370896029249848796293573094556559070111933328788951276428`90. for the integral and error estimates.
3.161570082363756029776335213343*10^-7

The problem is, as I take eps to be smaller and smaller it will reach a point where error occurs which is related to convergence issues. For eps=10^-9, I can still get a result without error. However, for eps=10^-10 error already occurs.

I have set the corresponding AccuracyGoal,PrecisionGoal,WorkingPrecision to try and go over this problem however it seems not to be working. I think the smallest eps that I might need is around eps=10^-15 in order to justify the accuracy of what I'm doing, so if I can just push the convergence to this limit then everything will be good.

Improve this question
$\endgroup$
5
  • $\begingroup$ If you need the integral in [0,zh] limits why do you expand the integration domain and to [0,zh+eps], where zh>0 and eps>0 and take the principal value? These two are not equivalent. It does make sense mathematically. $\endgroup$
    – yarchik
    Commented Aug 8, 2021 at 5:53
  • $\begingroup$ @yarchik The reason is because the integrand is divergent at zh. So I added an extremely small parameter eps, say, eps=10^-10. Then I would say it is approximately the same as the original integral to some tolerance eps. The integrand diverges from below for less than zh and diverges from above for greater than zh so they cancel near zh so I can take the principal value. $\endgroup$
    – mathemania
    Commented Aug 8, 2021 at 7:32
  • $\begingroup$ Let us say you have an improper integral $\int_0^{a} f(x)\, \mathrm{d}x$, where $f(x)\stackrel{x\rightarrow a}{\rightarrow a}$ and $a>0$. Then, the integral can be computed as $\lim_{b\rightarrow a^-}\int_0^{b} f(x)\, \mathrm{d}x$. See en.wikipedia.org/wiki/Improper_integral. This is different from what you are trying to compute $\endgroup$
    – yarchik
    Commented Aug 8, 2021 at 9:27
  • $\begingroup$ @yarchik I think you can just take that my problem is integrated in the range [0,zh+eps] originally then I just want to have convergence as I make eps smaller, say, eps=10^-15. $\endgroup$
    – mathemania
    Commented Aug 8, 2021 at 10:24
  • $\begingroup$ Ok, I plotted your integrand function. There seems to be 2 singularities. $\endgroup$
    – yarchik
    Commented Aug 8, 2021 at 10:31

1 Answer 1

Reset to default
2
$\begingroup$

A common approach, though it should perhaps be better known, is to subtract off the singular parts and integrate them separately. In this case, the integrand is complicated and the subtractive cancellation is significant. With work, that the loss of significance could probably be reduced. But Mathematica lets us be lazy and increase WorkingPrecision.

sc1 = SeriesCoefficient[tinteg[z, zhr + epsr, zhr], {z, zhr, -1}];
sint1 = sc1/(z - zhr);
sc2 = SeriesCoefficient[
   tinteg[z, zhr + epsr, zhr] - sint1, {z, zhr + epsr, -1/2}];
sint2 = sc2/Sqrt[z - (zhr + epsr)];
singpart = 
  Integrate[sint1, {z, 0, zhr + epsr}, PrincipalValue -> True, 
    Assumptions -> zhr > 0 && epsr > 0 && z > 0] + 
   Integrate[sint2, {z, 0, zhr + epsr}, 
    Assumptions -> zhr > 0 && epsr > 0 && z > 0];
th[eps_?NumericQ, zh_?NumericQ] := 
 Block[{epsr, zhr}, {epsr, zhr} = Rationalize[{eps, zh}, 0];
  singpart +
   NIntegrate[
    tinteg[z, zhr + epsr, zhr] - sint1 - sint2,
    {z, 0, zhr, zhr + epsr},
    AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, 
    MaxRecursion -> 9]]

OP's examples:

wp = 80;
th[10^-10, 100] // N[#, 20] &
(*  0.00010000000000008915927  *)

wp = 200;
Block[{$MaxExtraPrecision = 500},
  th[10^-15, 100]] // N[#, 20] &
(*  3.1622776601683869256*10^-7  *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ Similar strategy used here and here. $\endgroup$
    – Michael E2
    Commented Aug 9, 2021 at 4:12

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.