Problem with setting working precision in NIntegrate

Question

I want to obtain a good numerical approximation (up to 10 decimal place would be ok for me) to an integral:

$$ \int^{\infty}_{0} f(r)r^2dr $$

I am using the function $f(r)$, which is related to the function

$$g(r)=-\frac{\sqrt[3]{3} \sqrt[3]{e^{-2 r}}}{\pi ^{2/3}}-\frac{\sqrt[3]{2 \pi }}{5 \sqrt[3]{e^{-2 r}} \left(\frac{3 \sqrt[3]{\pi } \sinh ^{-1}\left(\frac{2 \sqrt[3]{2 \pi }} {\sqrt[3]{e^{-2 r}}}\right)}{5\ 2^{2/3} \sqrt[3]{e^{-2 r}}}+1\right)}$$

by the relation

$$ f(r)=-\frac{1}{4\pi}\nabla^2_{r,\theta,\phi} g(r) $$

Obviously, explicit integration is impossible. The product $f(r)r^2$ is well-behaved and integrable for sure. The function $f(r)$ decays faster than $\frac{1}{r^2}$.

When I try to increase WorkingPrecision in NIntegrate, Mathematica says the expression I am integrating itself is not specified so precisely. How can I overcome this? Any tips/ hints?

I am asking for a general strategy to obtain a precise value of the integral:

NIntegrate[f[r]*4*π*r^2, {r, 0, y}, WorkingPrecision -> x]   

where y and x are some numbers.

P.S., I've been using Mathematica for only two days.

Answer 1

Next time it would be really nice of you to give actual Mathematica code for your function definitions. Well, for now, thanks to @ssch and ToExpression we have the following definition of g[r]:

g[r_]=-((3^(1/3) (E^(-2 r))^(1/3))/\[Pi]^(2/3))-
(2 \[Pi])^(1/3)/(5 (E^(-2 r))^(1/3) (1+(3 \[Pi]^(1/3) 
ArcSinh[(2 (2 \[Pi])^(1/3))/(E^(-2 r))^(1/3)])/(5 2^(2/3) (E^(-2 r))^(1/3))));

g[r]//TraditionalForm

Find function f[r]:

f[r_] = Simplify[PowerExpand[-(1/(4 Pi)) D[r^2 g'[r], r]/r^2]];

Note, you need the simplifications above so NIntegrate won't choke on too many terms. Now your function has a cool asymptotic behavior:

Limit[r^4 f[r], r -> Infinity]

So indeed r^2 f[r] will decay nicely. Let's plot the monster you're trying to integrate:

Plot[4 Pi f[r] r^2, {r, 0, 25}, PlotRange -> All, Filling -> Axis]

Define the integrating function:

h[s_] := NIntegrate[f[r] 4 Pi r^2, {r, 0, s}, MaxRecursion -> 12]

It works fine:

h /@ {.1, 1, 10, 100, Infinity}

{-0.00396818, -0.2323, -0.647235, -0.942626, -0.999967}

You can even plot it:

ListPlot[Table[{s, h[s]}, {s, 0, 10, .1}], 
 PlotStyle -> Directive[Red, PointSize[.01]], Filling -> 0, AspectRatio -> 1/3]

Answer 2

This is not an answer, just too long for a comment.

To start of with your NIntegrate call looks peculiar, the f(r) part is not a function call but is interpreted as f*r a function call looks like f[r]

(* Obtained by doing ToExpression[...I pasted latex stuff here..., TeXForm]
 and replaced e with E ( /.e->E) *)
g[r_] := -((3^(1/3) (E^(-2 r))^(1/3))/\[Pi]^(2/3)) - (2 \[Pi])^(1/3)/(
   5 (E^(-2 r))^(
    1/3) (1 + (
      3 \[Pi]^(1/3) ArcSinh[(2 (2 \[Pi])^(1/3))/(E^(-2 r))^(1/3)])/(
      5 2^(2/3) (E^(-2 r))^(1/3))));
(* I presume the laplacian simplifies to 1/r^2 D[r^2 D[g[r], r], r] in this case *)
f[r_] := Evaluate[Simplify[ -1/(4 Pi) 1/r^2 D[r^2 D[g[r], r], r], r \[Element] Reals]]
(* Sanity check, gives 0, so that's good *)
Limit[f[r] r^2, r -> Infinity]

This is basically how far I got, this spits out errors about under/overflow and $MaxExtraPrecision being reached

Block[{$MaxExtraPrecision = 1000000.}, 
  N[f[10^20] (10^20)^2, {Infinity, 10}]]