7
$\begingroup$

I am trying to understand the numerical integration routine by using as a benchmark the function

f[x_] := -x^2 + x^4

used to define the following function

g[z_] := NIntegrate[1/Sqrt[-f[x]], {x, z, 1}, WorkingPrecision -> 16]

As $\lim_{x \to 0} g(x) = \infty$

I was checking how the computation is performed as the argument of $g$ approaches $0$.

Trying

g[10*^-32]

I get the warning

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive        bisections in x near {x} = {1.64072643374912601558232961177989984273356228783177323546417714519*10^-29}. NIntegrate obtained 71.8065640502810565832993434204331739193329356186861158845012728424`66. and 2.36220417967229179569077470511748314210153999062657561896751085287`66. for the integral and error estimates.

and a result of

71.80656405028106

For smaller arguments I get

g[10*^-100]
217.5365551442017

and an error estimate of

44.1017

I checked this against Julia and its standard integration package QuadGK

julia> j = quadgk(h,10^-100,1)
(230.9516545085585, 3.0963683972298146e-6)

no sweat, much smaller error.

In Python with Scipy.integrate I obtain similar results, but the maximum number of subdivision has to be increased, and I set it to 500.

Setting

MaxRecursion -> 500 

makes the calculation run for x10 more than Julia, and I still get some error (probably fully acceptable)

228.6490713869369

What should I improve in my handling of similar divergent integrals in Mathematica (assuming that Julia's claim are realistic..)?

Thanks

$\endgroup$
  • 1
    $\begingroup$ Mentioned in a comment to the answer, but 10*^-100 != 10^-100 (you want the latter) $\endgroup$ – b3m2a1 Apr 9 at 22:19
  • $\begingroup$ @b3m2a1, noted thanks $\endgroup$ – Smerdjakov Apr 10 at 7:49
7
$\begingroup$

Based on the warning, I would try increasing either the MaxRecursion or WorkingPrecision. Increasing the WorkingPrecision only seems to help marginally, so I went with MaxRecursion -> 20.

f[x_] := -x^2 + x^4
g[z_] := NIntegrate[
           1/Sqrt[-f[x]], 
           {x, z, 1}, 
           MaxRecursion -> 20, 
           WorkingPrecision -> 16
         ]
g[10^-100]

230.9516564798752

It does still give a warning, but the result seems to be correct to quite a few decimal places. We can verify this by performing the integration with infinite precision first:

Integrate[1/Sqrt[-f[x]], {x, 10^-100, 1}]

Log[10000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000] + 3 Sqrt[11111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111]

and

N[%, 16]

230.9516564799645

Also, in case you didn't already know, if you ever put a number like 10.^-100 that has a decimal, unless you explicitly tell Mathematica otherwise, it will assume that it is a machine precision number. If there is no decimal in the number (like in your numbers), it will assume it is exact and has infinite precision.

| improve this answer | |
$\endgroup$
  • $\begingroup$ well I certainly did not know. I was inputting g[10*^-100] and was getting $\approx 228.6$. Using g[10*^-100] gets better, thanks a lot $\endgroup$ – Smerdjakov Apr 9 at 21:53
  • 5
    $\begingroup$ @Smerdjakov Oh, that's another thing: 10*^-100 means 10*10^-100 which is actually just 10^-99. If you want to use that notation, I would use 1*^-100. $\endgroup$ – MassDefect Apr 9 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.