# Contradictory results of numerical evaluations of a function

I have the following function:

f[a_,b_,A_,m1_,x_] := Sum[2^(3 + 4 m) (1/a + 1/b)^(-2 (1 + m))Gamma[2 + 2 m] Hypergeometric1F1[2 (1 + m), 1/2, -((a b x^2)/(A^2 (a + b)))], {m, 0, m1}];


This is the result of the following integral:

Sum[Integrate[y^(4 m + 3) Cos[x/A y] E^(-y^2/4 (1/a + 1/b)), {y, 0, \[Infinity]}, Assumptions -> Re[m] > -1 && Re[1/a + 1/b] > 0], {m,0,m1}]


The limit of the function when $$x \rightarrow \infty$$ should be 0, as the terms of the sum go to $$0$$ as $$m \to \infty$$:

In: Limit[((2 a)/(a + b) - 1)^(2 m)   ((2^(3 - 2 m)) Gamma[2 + 2 m] )/Gamma[1 + m]^2, m -> \[Infinity], Assumptions -> Vxx > 0 && Vyy > 0]
Out: 0


and, using the asymptotic representation of the confluent hypergeometric function $$_1F_1(a,b,z)$$: $$\Gamma(b)\left(\frac{e^x x^{a - b}}{\Gamma(a)} + \frac{-x^{-a}}{\Gamma(b - a)}\right)$$

We get:

In: Limit[Gamma[b] ((E^x x^(a - b))/Gamma[a] + (-x)^-a/Gamma[b - a]) /. a -> 2 (1 + m) /. b -> 1/2 /. x -> -((x^2 a b)/(a + b)), x -> \[Infinity], Assumptions -> m >= 0 && Vxx > 0 && Vyy > 0 && A > 0]
Out: 0


So, for large $$x$$, I would expect the function to decrease to $$0$$.

Now, I have two different problems that, I suspect, are related.

1. When I evaluate the function at large values of $$x$$, I get contradictory results:

       In[1]: q = f[2, 3, 2, 3, x] // N;
q /. x -> 100
Out[1]: -5.76 + 0. I

In[2]: f[2, 3, 2, 3, 100] // N
Out[2]: 3.03645*10^26 + 0. I

In[3]: f[2., 3., 2., 3., 100.] // N
Out[3]: 9.61733*10^-7

In[4]: FullSimplify[f[2, 3, 2, 3, 100]] // N
Out[4]: -3.02833*10^29


And I tried to use the integral form to check, but the same thing happened:

       In[1]: Integrate[y^(4 m + 3) Cos[x/A y] E^(-y^2/4 (1/a + 1/b)) /. a -> 2. /. b -> 3. /. A -> 2. /. m -> 3. /. x -> 100., {y, 0, \[Infinity]}, Assumptions -> Re[m] > -1 && Re[1/a + 1/b] > 0]
Out[1]: 8.76703*10^-16

In[2]: Integrate[y^(4 m + 3) Cos[x/A y] E^(-y^2/4 (1/a + 1/b)) /. a -> 2 /. b -> 3 /. A -> 2 /. m -> 3 /. x -> 100, {y, 0, \[Infinity]}, Assumptions -> Re[m] > -1 && Re[1/a + 1/b] > 0]
Out[2]: (15479341056 (-1528415185117230896460916 + 30563207285404851812924775 Sqrt[30] DawsonF[10 Sqrt[30]]))/78125
In[3]: (15479341056 (-1528415185117230896460916 + 30563207285404851812924775 Sqrt[30] DawsonF[10 Sqrt[30]]))/78125 // N
Out[3]: -3.02833*10^29

In[4]: NIntegrate[FullSimplify[y^(4 m + 3) Cos[x/A y] E^(-y^2/4 (1/a + 1/b))] /. a -> 2 /. b -> 3 /. A -> 2 /. m -> 3 /. x -> 100, {y, 0, \[Infinity]}]
Out[4]: 104.409


For In[3], I got the following error:

And for In[4]:

2. When I plot the function(s), I get unexpected graphs:

First graph: f[2, 3, 2, x, 3] = s (blue); N[f[2, 3, 2, x, 3]] = N[s] (red); FullSimplify[N[f[2, 3, 2, x, 3]]] = FullSimplify[N[s]] (black). There are some oscillations at larger $$x$$, and different graphs for s, N[s] and FullSimplify[N[s]].

Second graph: f[2, 3, 2, x, 3] = s (blue); N[f[2, 3, 2, x, 3]] = N[s] (red). Very large oscillations at larger x, different graphs for s, N[s] and FullSimplify[N[s]], and we can see that the function stabilizes at a constant value, that's different from $$0$$.

Third graph: f[2, 3, 2, x, 3] = s (blue); N[f[2, 3, 2, x, 3]] = N[s] (red). Close-up of s and N[s], with different values.

Fourth graph: f[2, 3, 2, x, 3] = s (blue); N[f[2, 3, 2, x, 3]] = N[s] (red); FullSimplify[N[f[2, 3, 2, x, 3]]] = FullSimplify[N[s]] (black). Although s and N[s] stabilize at values different from $$0$$, FullSimplify[N[s]] seems to decrease (or grow) exponentially as $$x$$ gets larger, which contradicts the limit computed earlier.

As the function is an infinite sum in terms of $$m$$, as $$m$$ gets larger, the terms should be closer and closer together. But I plotted three variants of f[2,3,2,110,m] for values of $$m$$ going from $$0$$ to $$20$$, and this is the result:

Fifth graph: f[2,3,2,110,m]//FullSimplify//N; it grows exponentially as m gets larger.

Sixth (and final) graph: f[2.,3.,2.,110.,m] (blue); f[2,3,2,110,m] (red).

Here is the code for the graphs:

 (*1, 2, 3, 4*) Plot[Evaluate[{f[2, 3, 2, 3, t], f[2, 3, 2, 3, t]//N, f[2, 3, 2, 3, t]//N//FullSimplify}], {t, 0, 60}, ImageSize -> Large, PlotLegends -> {"s", "N[s]","FullSimplify[N[s]]"}, PlotStyle -> {{Red, Opacity[1]}, {Blue, Opacity[0.5]}, {Black, Opacity[1]}}]
(*5*) ListLogLogPlot[Table[f[2, 3, 2, n, 110] // N // FullSimplify, {n, 0, 20}]]
(*6*) Show[ListPlot[Table[f[2., 3., 2., n, 110.], {n, 0., 20.}]], ListPlot[Table[f[2, 3, 2, n, 110.], {n, 0, 20}], PlotStyle -> Red], PlotRange -> Automatic]


Is there an avoidable reason behind the strange behavior of Mathematica?

• Machine precision is insufficient for such small numbers. Don't use just N to obtain the numerical value; instead use arbitrary-precision by specifying a desired precision with e.g. N[..., 20] or however many digits of precision you think you need. Mar 17, 2021 at 15:19

With differences of numbers of the order of E^3000 you are in danger to overshoot machine precision.

You can see this e.g.:

f[a_, b_, A_, m1_, x_] :=
Sum[2^(3 + 4 m) (1/a + 1/b)^(-2 (1 + m)) Gamma[
2 + 2 m] Hypergeometric1F1[2 (1 + m),
1/2, -((a b x^2)/(A^2 (a + b)))], {m, 0, m1}];

q = f[2, 3, 2, 3, x];

f[2, 3, 2, 3, 100] == q /. x -> 100

(*True*)


Therefore, your calculation is correct up to where you use N. If subsequent calculations are needed, they must be done using arbitrary precision and not machine precision.