Follow up: how to evaluate this double sum quickly

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?

• You'd better add the definitions of μ and Ω and nmax into the code rather than state in the text of your question. Also, you forgot to add the definition of Ngrid (though it's not hard to guess that it's 160.) Finally, your code cotains typo, the With in pt1 loses its right bracket, and the position of WorkingPrecision option is wrong. – xzczd Oct 14 '15 at 10:34

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][];
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.