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.