# How to speed up the calculation of a multi-dimension matrix involving symbolic integral?

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}


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.

• 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 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 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 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 at 5:58
• @likehust that's difficult to answer in general. Maybe you can ask this as a separate question. – Roman May 26 at 6:02