0
$\begingroup$

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:enter image description here

And for In[4]: enter image description here

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$. enter image description here

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. enter image description here

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. enter image description here

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. enter image description here

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

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?

Improve this question
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Mar 17, 2021 at 15:19

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.