Follow up: how to evaluate this double sum quickly

Question

This question is a follow-up to my previous question.

The code I use at the moment is the following:

Bose[k_, μ_, nx_, ny_] := 1/(Exp[( k (nx + ny + 1) - μ)] - 1)
wavefunc[x_, n_] := 1/Sqrt[2^n Factorial[n] π] Exp[-x^2/2] HermiteH[n, x]
Intensity[k_, x_, y_, μ_, nmax_] := Sum[Bose[k, μ, nx, ny] (wavefunc[x, nx]*wavefunc[y, ny])^2, {nx, 0, nmax}, {ny, 0, nmax}]

Using this, I can now easily make plots with the following line (which was a suggestion in the answer to the previous question)

(pt1 = With[{expr = Intensity[Ω, x, 0, μ1, nmax]}, 
      Plot[expr, WorkingPrecision -> 50, {x, -40, 40}];)
Show[pt1]

This sums up $500^2$ terms in about 4 minutes if I take a large working precision. (Typical values are $\Omega = 1/150$ and $\mu = 1/2 * \Omega$).

However, I would like to export the graph as data points in a .dat file, such that I can use the data elsewhere. This again suffers from extremely long running times. The code I use is

linearmesh[a_, b_, n_Integer] := Array[# &, n, {a, b}]
xarr = linearmesh[-40, 40, 160];
soly1 = Monitor[Table[Intensity[Ω, xarr[[i]], 0, μ, 150], {i, 1,Ngrid}], i];
Export["fit.dat", soly1]

but even for 150 terms this takes forever. Is there an easy way to speed this up significantly?

Answer 1

Since you've already got pt1 in 4 minutes (well, your computer really owns a huge memory, my 8G laptop exceeds the memory limitation while ploting pt1!), you can just dig out those data points in pt1:

data = Cases[pt1, Line[a_] :> a, Infinity][[1]];
Export["fit.dat", data]
SystemOpen@Directory[]

My memory exhausted, so I can't test the code, but it should work.

Then I'd like to say something about your new trial. I think its slowness is mainly because of the use of exact number rather than arbitrary precision number. Indeed, though still a little faster than symbolic calculation, arbitrary precision calculation is already slow:

Ω = 1/150; μ = Ω/2;
nmax = 150;
eg1 = Intensity[Ω, 1`50, 0, μ, nmax]; // AbsoluteTiming 
eg2 = Intensity[Ω, 1, 0, μ, nmax]; // AbsoluteTiming 
(* {10.150728, Null} *)
(* {13.089546, Null} *)

But don't forget eg2 is a symbolic expression! When exported as a .dat, it becomes something like:

Export["a.dat", eg2] // SystemOpen

Pictured by Simon Wood's shadow.

Such expression is undoubtedly slow in exporting, and most of all, undesired for your subsequent processing. You can add something like N[…, 50] to soly1 to fix your new code.