# Wrong divergence with numerical value

this is my first question in this forum.

I'm trying to evaluate some complicated function of, say, $$x$$ near $$x=0$$ (in order to integrate it later). The problem is that the numerical value of this function when $$x \rightarrow 0$$ goes to infinity while I know for sure that the limit of this function at the origin is $$0$$. This is happening because the function has the following structure (but much more complicated):

$$f(x) \sim (A+\frac{B}{x}) x^2$$

and for some reason mathematica is sending the second term to infinity faster than $$x^2$$ goes to zero.

Here I post a screenshot of my output, where it can be seen that the limit in the origin is indeed $$0$$: Do you know any way of dealing with this problem?

Thank you.

Since I think that I wasn't accurate enough i will explain my problem better:

I'm trying to solve the following integral

$$f(\theta) \propto \int_0^\infty dq \frac{q}{(\mu^2+q^2)^2} \int_0^{2 \pi} dx \frac{q^2 \cos(x)^2}{[q^2+c_1-2q c_2 cos(x)]^2[q^2+c_1'-2q c_2' cos(x+\theta)]^2[q^2+c_1''+c_3'' cos(x+\theta)-2q c_2'' cos(x)]^2}$$

The code that I made solves the $$x$$ integral using the residue theorem. So, defining $$z=e^{ix}$$, $$a_1=\frac{q^2+c_1}{2 q c_2}$$, $$a_2=\frac{q^2+c_1'}{2 q c_2'}$$, $$a_3=\frac{q^2+c_1''}{2 q c_2''}$$, $$b=\frac{c_3''}{2 q c_2''}$$ and $$c=\theta$$ (where all the constants are real) the integral transforms to

$$\int_0^\infty dq \frac{q}{(\mu^2+q^2)^2q^2} \int_0^{2 \pi} dx \frac{C}{q^4} \frac{q^2 \cos(x)^2}{[a_1 - cos(x)]^2[a_2 - cos(x+c)]^2[a_3+b cos(x+c)-cos(x)]^2}$$,

where $$C$$ is some constant. The code that I have made to solve this integral is

z1 = (a3 + Sqrt[(a3)^2 - 1 - b^2 + 2 b Cos[c]])/((1 - b Exp[I c]));
z2 = (a3 - Sqrt[(a3)^2 - 1 - b^2 + 2 b Cos[c]])/((1 - b Exp[I c]));
z3 = a1 + Sqrt[a1^2 - 1];
z4 = a1 - Sqrt[a1^2 - 1];
z5 = (a2 + Sqrt[a2^2 - 1]) Exp[-I c];
z6 = (a2 - Sqrt[a2^2 - 1]) Exp[-I c];
f[z_] := ((16 E^(-4 I c)
z^2 A6 (1 + E^(2 I c) z^2)^2))/((z - z3)^2 (z - z4)^2 (z -
z5)^2 (z - z6)^2 ((1 - b E^(I c) ) (z - z1) (z - z2))^2) z
peq1 = 2 Pi (Residue[f[z], {z, z2}] + Residue[f[z], {z, z4}] +
Residue[f[z], {z, z6}]);


Check if the integral is fine:

In:= N[peq1 /. {A6 -> 6, a1 -> 2, a2 -> 2, a3 -> 3, b -> 1, c -> 2}]
NIntegrate[( 6 Cos[2 + x]^2)/((2 - Cos[x])^2 (2 -
Cos[2 + x])^2 (3 + 1 Cos[x + 2] - Cos[x])^2), {x, 0, 2 Pi}]

Out= 0.338923 + 8.99042*10^-13 I

Out= 0.338923

In:= prue = peq1 /. {A6 -> c3^2};

In:=
prue2[k1_, k2_, q_, \[Mu]_, \[Mu]p_, \[CapitalDelta]\[Phi]_] :=prue /. {c3 -> 1/(4 q k1), a1 -> (k1^2 + q^2 + \[Mu]p^2)/(2 k1 q),a2 -> (k2^2 + q^2 + \[Mu]p^2)/(2 k2 q), a3 -> (q^2 + k1^2 + k2^2 + \[Mu]^2 +2 k1 k2 Cos[\[CapitalDelta]\[Phi]])/(2 q k1), b -> -k2/k1,c -> -\[CapitalDelta]\[Phi]}

In:= example[q_] := prue2[1, 1, q, 1/10, 2/10, 1/10]

In:= N[example[1/10000]]

Out= -2.82812*10^10 - 7.3361*10^11 I

In:= Limit[example[q], q -> 0]

Out= 0


I think that the problem with the divergence have to be in the definition of the variables $$a_1$$, $$a_2$$, $$a_3$$ and $$a_4$$ but I don't how to solve the $$x$$ integral in a way that there is no $$1/q$$ factors.

Thank you all.

• Perhaps there is a precision error. Try converting all the numbers in both your function and in the argument to exact numbers (no decimals) and see what happens. For instance, you'd want to enter example[1/10000] instead of example[0.0001]. Alternately, you could use non-exact numbers, but specify a sufficiently high precision, e.g. example[0.000150]. Dec 15 '18 at 17:11
• Please provide the complete code that allows to fully reproduce the issue. Dec 15 '18 at 17:11
• A major part of numerics is about rearranging terms such that are evaluate in an order that is less susceptible to catastrophic cancellation. In the current case, f[x] := A x^2 + B x should work better. Dec 15 '18 at 17:25
• You may be able to use Collect to implement suggestion by @HenrikSchumacher. Also, since you are going to be integrating, if you have version 11.3, look at the EXPERIMENTAL function AsymptoticIntegrate Dec 15 '18 at 17:59
• Machine precision neither tracks/controls precision. You get what you get--but rapidly. That is often good enough. When you get unexpected results or just to check results, use arbitrary precision. When you calculate with arbitrary-precision numbers, Mathematica tracks precision at all points and tries to give you results which have the highest possible precision, given the precision of the input. Using your code, for N[example[1/10000], 15] I get 6.71359093775999*10^-17 + 0.*10^-32 I Dec 15 '18 at 18:32

If you make a plot you can see that Mathematica has trouble evaluating example[q] near q=0:

LogLinearPlot[Re[example[q]], {q, 0, 2}, PlotRange -> {-100000, 300000}]
LogLinearPlot[Im[example[q]], {q, 0, 2}, PlotRange -> All]  If the purpose is to establish a value to begin integration, it seems that you could pick a number around q = 0.1.

• Try LogLinearPlot[Re[example[q]], {q, 0, 2}, PlotRange -> {-100000, 300000}, WorkingPrecision -> 15] LogLinearPlot[Im[example[q]], {q, 0, 2}, PlotRange -> All, WorkingPrecision -> 15] Dec 15 '18 at 18:36
• @Bob Hanlon Yes, this removes the oscillations from the plot. But is there a way to achieve the same in evaluating the function? Dec 15 '18 at 20:34
• The plot is evaluating the function. Note that it is better to use exact input and then convert to numerical rather than use an approximate input. Compare Precision[example[1/10000]//N[#, 15]&] with Precision[example[N[1/10000, 50]]]` Dec 15 '18 at 20:51