# 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!

-
As a matter of course: computations that are stable when done in an orthogonal polynomial basis often become unstable if you re-express the orthogonal polynomials as powers of the argument. This would seem to be the case here. – J. M. is back. May 12 '13 at 4:04

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.

-
This is very good and helpful. Thanks a lot! It is not clear to me why in this case the matrix multiplication is more efficient than the sum. – lagoa May 14 '13 at 20:35
First, in the sum the coefficients are evaluated every time, compared to being evaluated once to construct the matrix. Second, some functions have internally optimized versions to use in commonly occurring situations. There are probably other reasons I'm not recalling right now. You might be interested in this and this. – Michael E2 May 14 '13 at 21:34

The problem has to do with precision in numerical approximations. During the evaluation of your function Mathematica has to deal with really huge numbers, so you have to make sure that during this computation everything is evaluated up to a very high precision. One way to achieve that is by specifying the initial argument with high precision. For example, if you evaluate

f[-4.0000000000000000000000000]


you will get the same answer as in the exact case.

-