# NIntegrate vs quadGK in Julia

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.806564050281056583299343420433173919332935618686115884501272842466. and 2.3622041796722917956907747051174831421015399906265756189675108528766. 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

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

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.

• 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 – Smerdjakov Apr 9 at 21:53
• @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. – MassDefect Apr 9 at 21:57