1
$\begingroup$

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}
$\endgroup$
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – likehust May 24 at 16:19
  • $\begingroup$ Careful, it gives l-m, not 1-m. It's a bad idea to use l as a variable for exactly this reason. $\endgroup$ – Roman May 24 at 17:23
  • $\begingroup$ 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. $\endgroup$ – Roman May 24 at 17:25
  • $\begingroup$ 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! $\endgroup$ – likehust May 26 at 5:58
  • $\begingroup$ @likehust that's difficult to answer in general. Maybe you can ask this as a separate question. $\endgroup$ – Roman May 26 at 6:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.