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.
