0

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.

mathemania
  • 713
  • 5
  • 12
  • 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. – yarchik Aug 08 '21 at 05:53
  • @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. – mathemania Aug 08 '21 at 07:32
  • 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 https://en.wikipedia.org/wiki/Improper_integral. This is different from what you are trying to compute – yarchik Aug 08 '21 at 09:27
  • @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. – mathemania Aug 08 '21 at 10:24
  • Ok, I plotted your integrand function. There seems to be 2 singularities. – yarchik Aug 08 '21 at 10:31

1 Answers1

2

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.162277660168386925610^-7  )
Michael E2
  • 235,386
  • 17
  • 334
  • 747