Overflow, indeterminate, or infinity errore in numerics

I need your help once again. I am trying to calculate some numerical stuff with Mathematica and always run into the following error:

NIntegrate::inumri: "The integrand [...] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{7.36315*10^13, 1.10447*10^14}}"

My integral is of the following form (details are negligible, has to be convergent in any case):

$$\int_0^\infty{\rm d}x \sum_{n=1}^{\infty} f_n(x)$$

and I need to cut off the integration and sum at some high value of x and n, respectively. The functions $$f_n(x)\sim \frac{a_n-b_n}{c_n-d_n}$$ look like this and diminish for higher values of n.

WithNIntegrate I can get around the error message by increasing the WorkingPrecision to a certain to a certain point. This was my first choice solution to the problem but then the whole program takes a lot of computational time, as I need to compute a huge amount of points.

In the following, I tried to replace the integral over x by a Gaussian quadrature to improve the computational time. With this solution I am running again into the problem of overflow. This time though, I have no control parameter like WorkingPrecision in the case of NIntegrate to get around the problem.

I suspect that at certain values of x and n the denominator consists of a substraction of very small numbers. Since Mathematica resorts to machine number precision if not specified otherwise, it is like dividing over zero.

I would like to stick to the idea of replacing the integral by a gaussian quadrature, but is there a way to control the precision of numbers in a sum?

I know about N[x, Precision] but it is not working in my sums. The program always sticks to machine precision instead, thus leading to the overflow error.

Edit

I think the code is too long to post it here, as it is a very longy and complicated function. f is the function under study, n the summation index. R is a parameter, om is the integration variable.

In: f[2, n, 1, 1, om] /. {n -> 1, R -> 3.16*10^-6, om -> 3*10^13}
Out: 0.000484331 + 0.0216804 I

For higher sum index n, the number gets very small or rather zero for the specified parameters R and om.

In: f[2, n, 1, 1, om] /. {n -> 100, R -> 3.16*10^-6, om -> 3*10^13}
Out: -5.057392809172020*10^-482 + 3.553535804695703*10^-477 I

And for even higher sum index, I get the error

In: f[2, n, 1, 1, om] /. {n -> 200, R -> 3.16*10^-6, om -> 3*10^13}
Out: Indeterminate

along with the following messages:

Power::infy: Infinite expression 1/(0. + 0. I) encountered. >>
Infinity::indet: Indeterminate expression (0. + 0. I) ComplexInfinity encountered. >>

The problem is that in the final program I need the sum to go to this high values of n, as there is a non negligible contribution for other values of parameters R and om.

Edit 2

The code is:

f[P_, l_, muM_, epsM_, om_] :=
-KroneckerDelta[P, 1] *
(mus[om]/muM SphericalBesselJ[l, Sqrt[epss[om] mus[om]] R om/c0] *
D[Rs SphericalBesselJ[l, Rs], Rs] -
SphericalBesselJ[l, Rs] D[Rstilde SphericalBesselJ[l, Rstilde], Rstilde]) /
(mus[om]/muM SphericalBesselJ[l, Sqrt[epss[om] mus[om]] R om/c0] *
D[Rs SphericalHankelH1[l, Rs], Rs] -
SphericalHankelH1[l,Rs]D[Rstilde SphericalBesselJ[l, Rstilde],Rstilde]) -
KroneckerDelta[P,2]*
(epss[om]/epsM SphericalBesselJ[l, Sqrt[epss[om] mus[om]] R om/c0] *
D[Rs SphericalBesselJ[l, Rs], Rs] -
SphericalBesselJ[l, Rs] D[Rstilde SphericalBesselJ[l, Rstilde], Rstilde]) /
(epss[om]/epsM SphericalBesselJ[l, Sqrt[epss[om] mus[om]] *
R om/c0]D[Rs SphericalHankelH1[l, Rs], Rs] -
SphericalHankelH1[l, Rs] D[Rstilde SphericalBesselJ[l, Rstilde], Rstilde])
/. {Rs-> Sqrt[epsM muM] R om/c0, Rstilde -> Sqrt[epss[om] mus[om]] R om/c0}

c0 = 299792458;
omp := (903./100.*eV)/hbar;
omt:=(267./100.*10.^-2.*eV)/hbar;
eV = 1602/1000*10^-19;
hbar = 1055/1000*10^-34;

epss[om_] := 1 - omp^2/(om(om + I omt))
mus[om_] := 1

Leading to the inderminate result posted above for the given parameters.

• This really needs full code in order to have any chance of a proper diagnosis. – Daniel Lichtblau Aug 5 '16 at 13:55
• First I would plug in numbers in {7.36315*10^13,1.10447*10^14} to see which is happening, Overflow[], Indeterminate, or Infinity. From there I could probably figure out why. What to do about it depends on how one could feasibly transform the function, at which point heed Daniel's comment. – Michael E2 Aug 5 '16 at 14:06
• Clearly you have a 0/0 situation: Fix it. (Sorry, can't tell you how without the code.) It might be like (Exp[x] - 1)/x /. x -> 0, which won't be handled by Simplify. But something like Simplify[f[2, n, 1, 1, om] /. n -> 200] /. {R -> 3.16*10^-6, om -> 3*10^13} might work. – Michael E2 Aug 6 '16 at 15:35
• It looks to me you need to become acquainted with Mathematica's arbitrary precision arithmetic features before proceeding with your work. Suggest you google on "Mathematica arbitrary precision arithmetic". – m_goldberg Aug 6 '16 at 19:35
• Can you come up with a simpler function which demonstrates the same problem? – Simon Woods Aug 6 '16 at 21:12