2

The following program succeeds in getting matrix CC, but it takes time badly, especially in the case varNumber becomes larger just as the following varNumber = 35. Who can speed up the process of calculation? Thanks!

ClearAll["Global`*"];
varNumber = 35; end = Infinity;
s1 = 112*^9; s2 = 25*^9; s3 = 15.1; s4 = 5.5*10^-9;

a[m_] := Exp[-x/2]*LaguerreL[m, x];
b[m_, i_, j_, l_] := Integrate[a[m]*x^i*D[a[l], {x, j}], {x, 0, end}];
 d[m_, i_, j_, l_] := 
 Integrate[
  a[m]*x^i*D[
    a[l], {x, j}]*(DiracDelta[x] - 
     DiracDelta[x - end]), {x, -Infinity, Infinity}];

c[1, 1][m_, l_] := s2*d[m, 0, 1, l] + s2*b[m, 0, 2, l];
c[1, 2][m_, l_] := 0;
c[1, 3][m_, l_] := 0;
c[2, 1][m_, l_] := 0;
c[2, 2][m_, l_] := s1*d[m, 0, 1, l] + s1*b[m, 0, 2, l];
c[2, 3][m_, l_] := s3*d[m, 0, 1, l] + s3*b[m, 0, 2, l];
c[3, 1][m_, l_] := 0;
c[3, 2][m_, l_] := s3*d[m, 0, 1, l] + s3*b[m, 0, 2, l];
c[3, 3][m_, l_] := -s4*d[m, 0, 1, l] - s4*b[m, 0, 2, l];

CC = ArrayFlatten@
    Table[c[m, n][i, j], {m, 3}, {n, 3}, {i, 0, varNumber - 1}, {j, 0,
       varNumber - 1}]; // AbsoluteTiming
{2283.69, Null}
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
likehust
  • 693
  • 3
  • 8

1 Answers1

7

Try this instead of your definitions of b and d:

b[m_, 0, 2, l_] /; l == m = 1/4;
b[m_, 0, 2, l_] /; l > m = l - m;
b[m_, 0, 2, l_] /; l < m = 0;
d[m_, 0, 1, l_] = -l - 1/2;

With these I can assemble CC in 0.047589 seconds (varNumber = 35).

For different values of $i$ and $j$, I find the fast definition

d[m_, i_, j_, l_] := If[i == 0, a[m]*D[a[l], {x, j}] /. x -> 0, 0];

following directly from the integral over Dirac $\delta$-functions. As for b[m,i,j,l] I would recommend asking at math.SE to see if anyone knows a closed formula for these integrals.

More generally, in the absence of such formulas you can use classical memoization, which gains you a big factor because you don't keep re-calculating the same values. Alternatively, you can use persistent memoization that will remember values forever, even when the kernel is restarted:

cacheloc = PersistenceLocation["Local", 
  FileNameJoin[{$UserBaseDirectory, "caches", "bintegrals"}]];
end = Infinity;
a[m_] = Exp[-x/2]*LaguerreL[m, x];
b[m_Integer, i_Integer, j_Integer, l_Integer] := b[m, i, j, l] =
  Once[Integrate[a[m]*x^i*D[a[l], {x, j}], {x, 0, end}], cacheloc];

Here I've used a combination of classical memoization and persistent storage, which complement each other: the former is very fast but impermanent, while the latter is a bit sluggish but permanent. Together, we get both advantages: the first lookup from permanent storage is still sluggish, but afterwards we get very fast lookup.

Roman
  • 47,322
  • 2
  • 55
  • 121
  • The speed of your codes surprised me in the extreme with correct answer, but I don't know why b[m_, 0, 2, l_] gives 1/4, 1-m, 0 in different conditions and d[m_, 0, 1, l_] gives -1-1/2, Can you explain the process of getting their body of b[ ] and d[ ] in different conditions, thank you. – likehust May 24 '19 at 16:19
  • Careful, it gives l-m, not 1-m. It's a bad idea to use l as a variable for exactly this reason. – Roman May 24 '19 at 17:23
  • Just look at the matrices Table[b[m, 0, 2, l], {m, 0, 10}, {l, 0, 10}] and Table[d[m, 0, 1, l], {m, 0, 10}, {l, 0, 10}], it's pretty obvious what the patterns are. I didn't do more than that. For higher values of i and j the formulas may be a bit more complex. – Roman May 24 '19 at 17:25
  • Thanks, I will avoid using l as argument as you suggested in the future, and have known how to get the body of functions b[ ] and d[ ]. However, I have to say, In fact, the program mentioned above is only a mini part of a long complex one, there exists lots of functions similar to b[ ] and d[ ], and their body are more complex than the listed above, so it's very difficult to find the patterns of them through eyes. Can you propose a more effective method? Thank you! – likehust May 26 '19 at 05:58
  • @likehust that's difficult to answer in general. Maybe you can ask this as a separate question. – Roman May 26 '19 at 06:02
  • Thanks ,your solution is great. I also want to know: why it doesn't support parallel calculation? After my rewriting CC = ArrayFlatten@Table to CC = ArrayFlatten@ParallelTable and run it, which gives lots of warning such as: (kernel 4) Get::noopen : "Cannot open "C:\\ Users "..., LocalObject::garbled : "LocalObject["file:///C:/Users/",Get::noopen : Cannot open put.wl., and so on. What's important is that the ParallelTable slows the running speed down. Can you edit the program again to support parallelization? – likehust May 27 '19 at 13:25
  • 1
    @likehust parallelization does not work well together with memoization. Maybe at some point I'll give it another try. For now I'll just caution you that throwing parallelization at a problem is often less effective than working on improving the algorithm. – Roman May 27 '19 at 13:28