I want to improve my algorithm of knowing eigenvalues of the matrix H(k)

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. 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. 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.

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] This is a quick answer, maybe someone else can offer other improvements.

• Thank you. But I have an error and it does not work. "Table::itform: 位置2における引数nは反復演算の正しい形式ではありません．" I use Mathematica 8. – Sakurai.JJ Jan 19 '18 at 9:19
• 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. – Chris K Jan 19 '18 at 16:50
• Thank you!!! It works!! – Sakurai.JJ Jan 20 '18 at 3:21

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[ 