numerics - Numerical sum does not give consistent results

Consider the function

λ[n_] = Sqrt[2 n];
ϕ[n_, x_] = 1/Sqrt[2^n n! Sqrt[π]] HermiteH[n, x];
f[x_] = Sum[2 Re[(-1)^(n + m) (1/(λ[m] + λ[n + 1]) + 1/(λ[n] + λ[m + 1]))
            If[n != m, 1, 1/2]] ϕ[n, x] ϕ[m, x] E^-x^2, {n, 0, 100}, {m, n, 100}];

Why is

f[-4.0] == 4.9697

while

f[-4] // N == 1.07681

?

The second result is correct as you can check by using SetDelayed in the definition of f. This discrepancy has to do with approximations. It would be good to fix it, since computations are faster if the function is defined with Set.

I would also appreciate suggestions to evaluate this sum more efficiently.

Thank you!

Answer

First of all, as you seem to indicate, you don't have the problem with accuracy if you use := (SetDelayed) in your function definitions instead of = (Set), although, as you say, with this change it takes longer to evaluate f. The accuracy improves because Mathematica will calculate the Hermite polynomials accurately when x is Real. Using Set for the definition of f means that the sum is evaluated with x being an undefined variable. In this case the polynomial expression for the Hermite polynomials is substituted, and they can have large coefficients with alternating signs, which lead to a extreme loss of precision. But since you also ask how to do the calculation more efficiently, and since the answer also improves accuracy, I will address that next.

---

A simple, efficient fix

Basically your function f is a matrix multiplication, which is more efficiently done using Dot. Also, the functions λ and φ thread over lists, efficiently and automatically:

λ[Range[0, 3]]
φ[Range[0, 3], 0.1]

{0, Sqrt[2], 2, Sqrt[6]}

{1/\[Pi]^(1/4), 0.106225, -0.520503, -0.129231}

Just how to implement these ideas depends on how much control over the precision of the calculation you wish to have. If, as in the examples in the question, you wish to compute both exact and machine-precision reals, then the following will do:

λ[n_] := Sqrt[2 n];
φ[n_, x_] := 1/Sqrt[2^n n! Sqrt[π]] HermiteH[n, x];
With[{mat = Table[If[n > m, 0, 
       2 Re[(-1)^(n + m) (1/(λ[m] + λ[n + 1]) + 1/(λ[n] + λ[m + 1])) If[n < m, 1, 1/2]]],
       {n, 0, 100}, {m, 0, 100}]},
   With[{nmat = N[mat]}, 
    f2[x_Real] := With[{φ0 = φ[Range[0, 100], x]}, E^-x^2 nmat . φ0 . φ0]];

   f2[x_] := With[{φ0 = φ[Range[0, 100], x]}, 
     E^-x^2 mat. φ0 . φ0];
  ];

Timings and results:

f2[-4.] // Timing

{0.009106, 1.07681}

N[f2[-4]] // Timing

{0.115590, 1.07681}

Compared with your original f, the function f2 is many times faster, especially on machine precision reals, as well as accurate:

f[-4.] // Timing

{0.385012, 4.9697}

N[f[-4]] // Timing

{0.762620, 1.07681}

f[-4.`25] // Timing (* high precision *)

{1.041199, 1.07681}

---

Arbitrary Precision

Now if you wish to be able to specify the precision, then the matrix nmat needs to have a precision just as great as the input x; but with the above definition, it has only machine precision (the default of N) when x is Real. The following gives three definitions of f, one for machine precision x for speed; one for an arbitrary precision real x; and one for any other x (e.g., exact x):

Clear[f2];
With[{mat = Table[If[n > m, 0, 
     2 Re[(-1)^(n + m) (1/(λ[m] + λ[n + 1]) + 1/(λ[n] + λ[m + 1])) If[n < m, 1, 1/2]]],
     {n, 0, 100}, {m, 0, 100}]},
  With[{nmat = N@mat},
   f2[x_?MachineNumberQ] := With[{φ0 = φ[Range[0, 100], x]}, E^-x^2 nmat. φ0. φ0]];

  f2[x_Real] := With[{nmat = N[mat, Precision[x]], φ0 = φ[Range[0, 100], x]}, 
    E^-x^2 nmat. φ0. φ0];
  f2[x_] := With[{φ0 = φ[Range[0, 100], x]}, E^-x^2 mat. φ0. φ0];
 ]

The timings for f[4.] and f[4] are essentially the same. In the example below, we specify 25 digits of precision in the input. The calculation results in almost 22 digits of accuracy. If we calculate the exact result to that precision, we get the same result (up to the precision specified).

f2[-4.`25] // Timing
Precision@Last@%

{0.239016, 1.076809573991586482393}
21.8196

N[f2[-4], %] // Timing

{2.570414, 1.076809573991586482393}

As an additional note, I wish to point out that using a specific precision, even $MachinePrecision, costs time, since Mathematica takes extra care to keep track of the precision of its calculations.

f2[N[-4, $MachinePrecision]] // Timing

Precision@Last@%

{0.235853, 1.076809573992}
12.7742

Standard machine precision does not report any lost of precision, even though there is in fact a loss of about 3 digits:

f2[-4.] // Timing
N@Precision@Last@%

{0.007986, 1.07681}

15.9546

---

The precision of the original function

Now returning to the original problem in which f[4.] returned an incorrect value, we examine the precision.

f[N[-4, 25]] // Timing
Precision@Last@%

{1.068105, 1.07681}

5.2808

In other words, almost 20 digits of precision are lost. Since machine precision is about 16 digits, it is no wonder that the answer is inaccurate. (Note below that when Mathematica tracks precision, the answer 0*10^5 is different; in fact, there is an error indicated, "No significant digits are available to display.")

f[N[-4, $MachinePrecision]] // Timing
Precision@Last@%

{1.037480, 0.*10^5}
0.