2
$\begingroup$

I want to improve my algorithm. I want to know the eigenvalues of below n-dimensional matrix H. When n=10, the matrix H is like below. enter image description here

Now, $To$ and $M$ are constants and $T1$ and $T2$ are functions as k.

$T1=2t_1\cos{\frac{k}{2}}$

$T2=2t_2\cos{\frac{k}{2}}$

$t1$ and $t2$ are constants. I want to know the eigenvalues of matrix H(k) each k numerically. And I want to plot the eigenvalues of H(k) as the function of k. So I wrote

n = 40;
M = 0.4;
Clear[X]
X = {{0, 0}};
For[k = -3.14, k < 6.28, k = k + 0.02,
To = 1.86;
T1 = 2*0.63*Cos[k/2];
T2 = 2*0.35*Cos[k/2];
H = 
Table[If[i == j + 1 && Mod[i, 2] == 0, T1, 0], {i, n}, {j, n}] +
Table[If[i + 1 == j && Mod[i, 2] == 1, T1, 0], {i, n}, {j, n}]+
Table[If[i == j && (Mod[i, 4] == 1 || Mod[i, 4] == 2), +M, 0], {i, n}, {j, n}] +
Table[If[i == j && (Mod[i, 4] == 3 || Mod[i, 4] == 0), -M, 0], {i, n}, {j, n}] +
Table[If[i == j + 1 && Mod[i, 2] == 1, To, 0], {i, n}, {j, n}] +
Table[If[i + 1 == j && Mod[i, 2] == 0, To, 0], {i, n}, {j, n}] +
Table[If[i == j + 3 && Mod[i, 2] == 1, T2, 0], {i, n}, {j, n}] +
Table[If[i + 3 == j && Mod[i, 2] == 0, T2, 0], {i, n}, {j, n}] ;
For[q = 1, q < n + 1, q++,X = Prepend[X, {k, Eigenvalues[N[H]][[q]]}]]]

Now, I can get eigenvalues of H each k numerically. I have a list of {k, eigenvalues}. So,

ListPlot[X]

finally, we can get the graph like below.

enter image description here

I want to improve my algorithm. When n is 100, it takes time. How should I improve my algorithm? I want to speed up the calculation.

I use Mathematica 8.

$\endgroup$
4
$\begingroup$

The biggest speed-up comes from constructing H only once, not every time in the main loop. I also avoided the For loop in saving your results, which repeatedly called Eigenvalues.

n = 100;
M = 0.4;
Clear[X, T1, T2]
X = {{0, 0}};
H = Table[If[i == j + 1 && Mod[i, 2] == 0, T1, 0], {i, n}, {j, n}] + 
 Table[If[i + 1 == j && Mod[i, 2] == 1, T1, 0], {i, n}, {j, n}] + 
 Table[If[i == j && (Mod[i, 4] == 1 || Mod[i, 4] == 2), +M, 0], {i, n}, {j, n}] + 
 Table[If[i == j && (Mod[i, 4] == 3 || Mod[i, 4] == 0), -M, 0], {i, n}, {j, n}] + 
 Table[If[i == j + 1 && Mod[i, 2] == 1, To, 0], {i, n}, {j, n}] + 
 Table[If[i + 1 == j && Mod[i, 2] == 0, To, 0], {i, n}, {j, n}] + 
 Table[If[i == j + 3 && Mod[i, 2] == 1, T2, 0], {i, n}, {j, n}] + 
 Table[If[i + 3 == j && Mod[i, 2] == 0, T2, 0], {i, n}, {j, n}];
For[k = -3.14, k < 6.28, k = k + 0.02, To = 1.86;
 T1 = 2*0.63*Cos[k/2];
 T2 = 2*0.35*Cos[k/2];
 ev = Eigenvalues[N[H]];
 X = Join[X, Transpose[{Table[k, {n}], ev}]]
]

takes only 0.57sec, with n=100.

ListPlot[X]

Mathematica graphics

This is a quick answer, maybe someone else can offer other improvements.

$\endgroup$
  • 1
    $\begingroup$ Thank you. But I have an error and it does not work. "Table::itform: 位置2における引数nは反復演算の正しい形式ではありません." I use Mathematica 8. $\endgroup$ – Sakurai.JJ Jan 19 '18 at 9:19
  • $\begingroup$ Whoops, I guess that the form Table[k, n] was added after v.8. I changed it to Table[k, {n}] above, which might work better. $\endgroup$ – Chris K Jan 19 '18 at 16:50
  • $\begingroup$ Thank you!!! It works!! $\endgroup$ – Sakurai.JJ Jan 20 '18 at 3:21
3
$\begingroup$

First, generate the appropriate number of values for k and the associated parameters namely T1 and T2.

Next, in order to overcome the fact that the implemented approach re-constructs matrix H for every instance of parameter k, the For 'loop' will be replaced by an appropriate and more efficient choice.

Namely, it seems reasonable to allow the main evaluation loop of Mathematica to take care of the various conditions required to construct matrix H (another solution would be to take advantage of the structure of H but that requires more effort and won't be pursued here, even though it might have potential performance benefits).

This solution, will generate the Eigenvalues of matrix H sequentially (as is the case in the code presented in the question) and finally will plot the list of eigenvalues.

Matrix H can be generated with one of the following equivalent ways (the last one seems to evaluate a little faster contrary to my presumption that the main loop would be more efficient):

Define separately the various cases that make up matrix H:

f[i_, j_] /; i == j + 1 && Mod[i, 2] == 0 := T1
f[i_, j_] /; i + 1 == j && Mod[i, 2] == 1 := T1
f[i_, j_] /; i == j && (Mod[i, 4] == 1 || Mod[i, 4] == 2) := M
f[i_, j_] /; i == j && (Mod[i, 4] == 3 || Mod[i, 4] == 0) := -M
f[i_, j_] /; i == j + 1 && Mod[i, 2] == 1 := To
f[i_, j_] /; i + 1 == j && Mod[i, 2] == 0 := To
f[i_, j_] /; i == j + 3 && Mod[i, 2] == 1 := T2
f[i_, j_] /; i + 3 == j && Mod[i, 2] == 0 := T2
f[i_, j_] := 0

or-more or less-equivalently:

g[i_, j_] := T1 /; i == j + 1 && Mod[i, 2] == 0
g[i_, j_] := T1 /; i + 1 == j && Mod[i, 2] == 1
g[i_, j_] := M /; i == j && (Mod[i, 4] == 1 || Mod[i, 4] == 2)
g[i_, j_] := -M /; i == j && (Mod[i, 4] == 3 || Mod[i, 4] == 0)
g[i_, j_] := To /; i == j + 1 && Mod[i, 2] == 1
g[i_, j_] := To /; i + 1 == j && Mod[i, 2] == 0
g[i_, j_] := T2 /; i == j + 3 && Mod[i, 2] == 1
g[i_, j_] := T2 /; i + 3 == j && Mod[i, 2] == 0
g[i_, j_] := 0

or by using Which:

h[i_, j_] := Which[
  i == j + 1 && Mod[i, 2] == 0, T1,
  i + 1 == j && Mod[i, 2] == 1, T1,
  i == j && (Mod[i, 4] == 1 || Mod[i, 4] == 2), M,
  i == j && (Mod[i, 4] == 3 || Mod[i, 4] == 0), -M,
  i == j + 1 && Mod[i, 2] == 1, To,
  i + 1 == j && Mod[i, 2] == 0, To,
  i == j + 3 && Mod[i, 2] == 1, T2,
  i + 3 == j && Mod[i, 2] == 0, T2,
  True, 0
 ]

Finally, use the eigenvalues of H to produce the required plot:

With[{n = 40},
  Module[{rng = Range[n], evs},

     (* generate matrix H *once* *)
     H[T1_, T2_] = Outer[h[##] &, rng, rng] /. Thread[{M -> 0.4, To -> 1.86}];

     (* for successive values of k, replace appropriate values for 
        T1 an T2, obtain the eigenvalues and assign k as the 
        x-coordinate *)
     evs = Table[
       Thread[
         {k, Eigenvalues[H[T1, T2] /. {T1 -> 2 0.63 Cos[k/2], T2 -> 2 0.35 Cos[k/2]}]}], 
      {k, -3.14, 6.28, 0.02}];

     (* produce separate plots for the n-th eigenvalue *)
     Grid[
       Partition[
         MapIndexed[
           ListPlot[#1, PlotLabel -> Row[{"eigenvalue", Null, #2[[-1]]}]] &,
           Transpose[evs]], 5, 5, {1, 1}, Null]]
    ]
  ]

enter image description here

$\endgroup$

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.