2
$\begingroup$

I am trying to solve an eigenvalue problem of a large matrix. The issue is that it takes too much time to construct this large matrix. This is the code that I built:

msa = 1.711*10^6;
msb = 0;
mu0 = 1.256*10^-6;
h0 = 0.1/mu0;
aa = 8.3*10^-12;
ab = 0;
alpa = 0.0019;
alpb = 0;
a = 10*10^-9;
rcyl = Sqrt[a^2/(2*Pi)];

nmax = 10;
n2 = 4;
ncap = (2*n2 + 1)^2;

gx[n_] := 2*n*Pi/a; 
gy[m_] := 2*m*Pi/a;
gxy = Table[{gx[n], gy[m]}, {n, -n2, n2}, {m, -n2, n2}];
g = ArrayReshape[gxy, {(2*n2 + 1)^2, 2}];

msg[{x_, y_}] := If[x == 0 && y == 0, (msa + msb)/2, (msa - msb)* BesselJ[1, (Sqrt[x^2 + y^2]*rcyl)]/(Sqrt[x^2 + y^2]*rcyl)];

q[{x_, y_}] := If[x == 0 && y == 0, ((2*aa/(mu0*h0*msa)))/2, ((2*aa/(mu0*h0*msa)))* BesselJ[1, (Sqrt[x^2 + y^2]*rcyl)]/(Sqrt[x^2 + y^2]*rcyl)];

sum1[kx_, ky_, i_, j_] := Sum[msg[g[[i]] - g[[l]]]* q[g[[l]] - g[[j]]]*(({kx, ky} + g[[j]]).({kx, ky} + g[[l]]) - (g[[i]] - g[[j]]).(g[[i]] - g[[l]])), {l, (2*n2 + 1)^2}];

bxy[kx_, ky_] := Table[If[(g[[j]] + {kx, ky}) == {0, 0}, KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]/(2*h0) + sum1[kx, ky, i, j], KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]*((g[[j, 2]] + ky)*(g[[j, 2]] + ky))/(h0*((g[[j]] + {kx, ky}).(g[[j]] + {kx, ky}))) + sum1[kx, ky, i, j]], {i, 1, (2*n2 + 1)^2}, {j, 1, (2*n2 + 1)^2}];

Next, I used arbitrary kx and ky to construct this bxy matrix. However, because of sum1, it takes a lot of time to construct it. When I used $n2=4$ and $kx=ky=\dfrac{2\pi}{3a}$,

k1=Pi*2/(3*a);
k2=Pi*2/(3*a);

Absolutetiming[bxy[k1,k2]]
(*{501.605,{{9.55609*10^8,4.84554*10^8...}}*)

It took about 502 seconds. It may not look bad, but later, I want to use $n2$ more than 10. Considering that the time would increase exponentially, I cannot even imagine how long it will take with $n2=10$.

After some investigation, I realized that it is caused by sum1 function.

bxy[kx_, ky_] := Table[If[(g[[j]] + {kx, ky}) == {0, 0}, KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]/(2*h0), KroneckerDelta[i, j] + msg[g[[i]] - g[[j]]]*((g[[j, 2]] + ky)*(g[[j, 2]] + ky))/(h0*((g[[j]] + {kx, ky}).(g[[j]] + {kx, ky})))], {i, 1, (2*n2 + 1)^2}, {j, 1, (2*n2 + 1)^2}];

k1 = Pi*2/((3)*a);
k2 = Pi*2/((3)*a);

AbsoluteTiming[bxy[k1, k2]]
(*{2.29623,{{6.37254,1.62134...}}}*)

So without sum1 function, it only takes about 2 seconds. I wonder if there is any way to speed up the summation. Thank you in advance!

P.S. I am not really familiar with this forum, so if I did not follow any rules, please let me know.

Improve this question
$\endgroup$
0
$\begingroup$

At the moment, a lot of the computations are performed in exact arithmetic although you want finite precision results in the end. Making this tiny change 

 a = N[10 10^-9];

along with defining gby

 g = Tuples[{gx /@ Range[-n2, n2], gy /@ Range[-n2, n2]}];

should speed up computations a lot.

Moreover, you compute many quantities involving the Bessel function multiple times. The Bessel function is rather expensive, so one should its evaluation to a minimum. In the code below, I do that by setting the matrices msgmatrix and qmatrix.

Moreover, you can use replace the Table in bxy by ParallelTable.

Here is the modified code; it executes in about 1.32 seconds on my Quad Core machine:

msa = 1.711 10^6;
msb = 0.;
mu0 = 1.256 10^-6;
h0 = 0.1/mu0;
aa = 8.3 10^-12;
alpa = 0.0019;
a = N[10 10^-9];
rcyl = Sqrt[a^2/(2 Pi)];

nmax = 10;
n2 = 4;
ncap = (2 n2 + 1)^2;

gx[n_] := 2  Pi/a n;
gy[m_] := 2 Pi/a m;
g = Tuples[{gx /@ Range[-n2, n2], gy /@ Range[-n2, n2]}];

msg[{x_, y_}] := If[Sqrt[x^2 + y^2] < 1. 10^-8,
   (msa + msb)/2,
   (msa - msb) BesselJ[
      1, (Sqrt[x^2 + y^2] rcyl)]/(Sqrt[x^2 + y^2] rcyl)
   ];

q[{x_, y_}] := If[Sqrt[x^2 + y^2] < 1. 10^-8,
   (2 aa/(mu0 h0 msa))/2,
   (2 aa/(mu0 h0 msa)) BesselJ[
      1, (Sqrt[x^2 + y^2] rcyl)]/(Sqrt[x^2 + y^2] rcyl)
   ];
msgmatrix = Table[msg[gi - gj], {gi, g}, {gj, g}];
qmatrix = Table[q[gi - gj], {gi, g}, {gj, g}];
sum1[kx_, ky_, i_, j_] := Sum[
   msgmatrix[[i, l]] msgmatrix[[l, j]] (({kx, ky} + g[[j]]).({kx, ky} + g[[l]]) - (g[[i]] - g[[j]]).(g[[i]] - g[[l]])),
   {l, ncap}];

bxy[kx_, ky_] := ParallelTable[
   If[Norm[(g[[j]] + {kx, ky})] < 10^-8,
    KroneckerDelta[i, j] + msgmatrix[[i, j]]/(2 h0) + sum1[kx, ky, i, j], 
    KroneckerDelta[i, j] + msgmatrix[[i, j]] (g[[j, 2]] + ky)^2/(h0 ((g[[j]] + {kx, ky}).(g[[j]] + {kx, ky}))) + sum1[kx, ky, i, j]],
   {i, 1, ncap}, {j, 1, ncap}];

k1 = N[Pi 2/(3*a)];
k2 = N[Pi 2/(3*a)];
A = bxy[k1, k2];

Of course, a lot more can be done. For example, sum1 or bxy could be compiled with Compile...

Improve this answer
$\endgroup$
3
  • $\begingroup$ I did not know that Bessel functions are heavy in mathematica. Thanks for your help! Now it only takes 5 secs to build the matrix which is really great! Also, I have one question. Is there any reason why you replaced 0 in if statement to a very small number (10^-8)? $\endgroup$
    – RIck
    Feb 7 '19 at 18:14
  • $\begingroup$ You're welcome. I used the tolerance 10^-8 because it is usually not a good idea to compare inexact number by Equal (due to rounding errors). One should always use certain tolerances (in fact, Equal also uses such a tolerance for inexact numbers, but may be too small for certain applications IIRC). So I used a tolerance just out of a habit. The number 10^-8 is somewhat arbitrary but it is certainly suffient to separate your grid points (since your grid is much coarser than that). $\endgroup$
    – Henrik Schumacher
    Feb 7 '19 at 18:29
  • $\begingroup$ That makes perfect sense. Thanks again! $\endgroup$
    – RIck
    Feb 8 '19 at 7:35

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.