1
$\begingroup$

I'm evaluating an integral using NIntegrate, I checked the expression and there shouldn't be any singularity at least on paper. The parameter zs can take on values in [0,9.9] and zh=10. I took zs=9 as an example, but when I evaluate it gave me an error about singularity/WorkingPrecision. I have no idea why this is happening since I already tried to adjust WorkingPrecison and other ways. Any help?

d = 3;
ag = 8;
pg = 8;
wp = 30;
f[z_, zh_] := 1 - (z/zh)^(d + 1);
Sinteg[z_?NumericQ, zs_?NumericQ, zh_?NumericQ] := Module[{zr, zsr, zhr}, {zr, zsr, zhr} = Rationalize[{z, zs, zh}, 0]; SetPrecision[(Sqrt[1/(f[zr, zhr] (1 - (zr/zsr)^(2 d) (1 + ((2 d f[zsr, zhr])/(d + 1)) (zhr/zsr)^(d + 1) (1 - (f[zsr, zhr]/f[zr, zhr])))))]), wp]]
SFull[zs_?NumericQ, zh_?NumericQ] := Module[{zsr, zhr}, {zsr, zhr} = Rationalize[{zs, zh}, 0]; NIntegrate[SetPrecision[Sinteg[z, zsr, zhr], wp], {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 300]]

SFull[9,10]
(*During evaluation of In[214]:= 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.
During evaluation of In[214]:= NIntegrate::zeroregion: Integration region {{9.00000000000000000000000000000,9.0000000000000000000000000000000000000000000000000000000000000000000000000000000}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions.
330.387118746134454732010635945*)
```
$\endgroup$
4
  • 1
    $\begingroup$ Your integrand diverges as $1/(z_s-z)$ as $z\to z_s$, and so the integral is unbounded. More precisely, for $z\to z_s$ the integrand is $\frac{z_h^2 z_s}{\sqrt{6} (z_s-z) \sqrt{5 z_h^4-z_s^4}} + O(1)$. $\endgroup$
    – Roman
    May 30, 2022 at 10:48
  • $\begingroup$ @Roman The original expression is under a square root sign so it should be finite right under an integral? As $z \rightarrow z_s$ the denominator approaches $\sqrt{0}$. $\endgroup$
    – mathemania
    May 30, 2022 at 10:59
  • $\begingroup$ @Roman Did you use Mathematica to expand out the integrated as $z \rightarrow zs$? If yes, how did you do it? $\endgroup$
    – mathemania
    May 30, 2022 at 15:29
  • $\begingroup$ See answer below, too long for a comment. $\endgroup$
    – Roman
    May 30, 2022 at 15:46

1 Answer 1

4
$\begingroup$

As per comment, the integrand diverges for $z\to z_s$:

d = 3;
f[z_, zh_] = 1 - (z/zh)^(d + 1);
Sinteg[z_, zs_, zh_] = Assuming[0 < z < zs < zh, 
  Sqrt[1/(f[z, zh] (1 - (z/zs)^(2 d) (1 + ((2 d f[zs, zh])/(d + 1))
    (zh/zs)^(d + 1) (1 - (f[zs, zh]/f[z, zh])))))] // FullSimplify]

$$ \frac{\sqrt{2} z_h^2 z_s^5}{\sqrt{3 z^{10} z_h^4+3 z^6 z_s^8+2 z_s^{10} \left(z_h^4-z^4\right)-z^6 z_s^4 \left(z^4+5 z_h^4\right)}} $$

Assuming[0 < z < zs < zh, 
  Series[Sinteg[z, zs, zh], {z, zs, 0}] // FullSimplify]

$$ -\frac{z_h^2 z_s}{(z-z_s) \sqrt{30 z_h^4-6 z_s^4}}+\frac{65 z_h^6-17 z_h^2 z_s^4}{6 \sqrt{6} \left(5 z_h^4-z_s^4\right)^{3/2}}+O(z-z_s) $$

$\endgroup$

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.