# Plotting eigenvalues with different colours

Bug introduced in version 11

I would like to plot all the eigenvalues (functions of Bz) of a matrix in the same plot, so as for Mathematica to understand that they are different functions and to plot in different colours.

Basically I used to do:

Plot[Evaluate[eigenvalues],{Bz, 0, 10}]


which used to work until I changed to Mathematica 11.0. Now it plots nothing.

Now I have to work around this by doing:

   Plot[eigenvalues /. Bz -> a, {a, 0, 10}]


which gives all the eigenvalues of the same colour. I can work with this but I'd like to understand why it's not plotting anymore.

Plotting the eigenvalues one by one I noticed that not all eigenvalues are not plottable, but only the ones where random '#' appeared in the expression.

For example, after running the code underneath, try to look at

eigenvalues[[7]]


It has #'s in it(why???), but if you do

 eigenvalues[[7]]/.Bz->100 //N


it still gives a finite answer!

--

FULL CODE:

MHz = 10^6;
h = 6.62607004*10^-34;  (* J\[Times]s *)
hbar = 1.054571800*10^-34; (* J\[Times]s *)

\[Mu]B = 9.27400968*10^-24; (* J\[Times]T^-1 *)
\[Mu]B = \[Mu]B/h*10^-4/10^6; (* G/MHz *)

splus[0] = {{0}} // SparseArray;
splus[S_] :=
splus[S] =
SparseArray[
Band[{1, 2}] ->
Table[Sqrt[S (S + 1) - M (M + 1)], {M, -S, S - 1, 1}], {2 S + 1,
2 S + 1}]
sminus[S_] :=
Transpose[
splus[S]] (* S- is just the transpose of S+ because have put the i \
in the definition of sx etc. *)
sx[S_] := (splus[S] + sminus[S])/2
sy[S_] := (splus[S] - sminus[S])/(2 I)
sz[S_] := SparseArray[
Band[{1, 1}] -> Range[-S, S, 1], {2 S + 1, 2 S + 1}]
SparseIdentityMatrix[n_Integer /; n >= 1] :=
SparseArray[Band[{1, 1}] -> 1, {n, n}]
id[S_] := id[S] = SparseIdentityMatrix[2 S + 1]

Ahfs = 6.093 ; (* MHz *)
Bhfs = 2.786 ; (* MHz *)
gJ = 4/3; (* ground *)
gI = -0.00014193489;
gS = 2.0023193043737;
gL = 0.99998627;

Ix = KroneckerProduct[sx[3/2], id[1/2], id[1]] ;(* tensor products *)

Iy = KroneckerProduct[sy[3/2], id[1/2], id[1]];
Iz = KroneckerProduct[sz[3/2], id[1/2], id[1]];
Itot = {Ix, Iy, Iz};

Sx = KroneckerProduct[id[3/2], sx[1/2], id[1]];
Sy = KroneckerProduct[id[3/2], sy[1/2], id[1]];
Sz = KroneckerProduct[id[3/2], sz[1/2], id[1]];
Stot = {Sx, Sy, Sz};

Lx = KroneckerProduct[id[3/2], id[1/2], sx[1]];
Ly = KroneckerProduct[id[3/2], id[1/2], sy[1]];
Lz = KroneckerProduct[id[3/2], id[1/2], sz[1]];
Ltot = {Lx, Ly, Lz};

Jtot = {Jx, Jy, Jz} = Stot + Ltot;
Ftot = {Fx, Fy, Fz} = Itot + Jtot;

Fs = Fx.Fx + Fy.Fy + Fz.Fz;
Js = Jx.Jx + Jy.Jy + Jz.Jz;

Hhfs = Ahfs*Sum[Itot[[i]].Jtot[[i]], {i, 1, 3}] + \[Mu]B*
Bz*(gI*Iz + gS*Sz + gL*Lz); (* Bz in Gauss *)

{eigenvalues, eigenvectors} =
Eigensystem[
SetPrecision[{Hhfs}, \[Infinity]]];(* Set that to \[Infinity] \
otherwise won't work *)

eigenvectors = Normalize /@ eigenvectors;

{eigenvalues,
eigenvectors} = {{eigenvalues[[1]], eigenvalues[[2]],
eigenvalues[[3]], eigenvalues[[5]], eigenvalues[[6]],
eigenvalues[[8]], eigenvalues[[9]], eigenvalues[[11]],
eigenvalues[[13]], eigenvalues[[14]], eigenvalues[[16]],
eigenvalues[[18]], eigenvalues[[19]], eigenvalues[[20]],
eigenvalues[[22]], eigenvalues[[24]]}, {eigenvectors[[1]],
eigenvectors[[2]], eigenvectors[[3]], eigenvectors[[5]],
eigenvectors[[6]], eigenvectors[[8]], eigenvectors[[9]],
eigenvectors[[11]], eigenvectors[[13]], eigenvectors[[14]],
eigenvectors[[16]], eigenvectors[[18]], eigenvectors[[19]],
eigenvectors[[20]], eigenvectors[[22]], eigenvectors[[24]]}};

• A tip: IdentityMatrix[n, SparseArray] is a shorter way to generate a sparse identity matrix. I haven't tried running your code, but I suspected something when you were complaining about #; have you already seen this? – J. M. will be back soon Jan 26 '17 at 14:08
• Works well on 10.4.1 (that's another reason for me to not switch to v11 yet). – corey979 Jan 26 '17 at 18:20
• As @corey979 mentioned, I've had a chance to try your original code on MMA 10.4 at work, and it works fine there. This must be a bug introduced in version 11. – MarcoB Jan 27 '17 at 17:48

At this time, this is only an extended comment.

When playing around with your problem, I noticed that the first 8 eigenvalues from the full eigenvalues list returned by Eigensystem plot just fine, with colors, using the expected syntax:

Plot[Evaluate@eigenvalues[[1 ;; 8]], {Bz, 0, 10}]


All the other ones don't plot at all:

MapIndexed[
Plot[
#1, {Bz, 0, 10},
Epilog -> Inset[Style[First@#2, Red, 18], Scaled[{0.9, 0.4}]]
] &,
eigenvalues
]


But there seems to be nothing inherently wrong with these expressions, nor does this behavior depend on the presence of Root objects: those are also present in the first 8, and those plot fine!

In fact, one can successfully reconstruct the full plot point by point:

ListPlot@
Table[{Bz, eigenvalues[[i]]}, {i, Length@eigenvalues}, {Bz, 0, 10}]


After confirming that the original code works with no issues in MMA 10.4, I also found that Plot seem to calculate at least some of the points to display; it just doesn't display them:

Function[{arg},
ListPlot @ #[[-1, -1]]& @
Reap@
Plot[
arg, {Bz, 0, 10},
EvaluationMonitor :> Sow[{Bz, arg}],
PlotRange -> All
]
] /@ eigenvalues


Something funny is still going on, though: notice that some of the plots do not cover the whole requested range of Bz; a manual larger PlotPoints setting will force the calculation of all points, but nonetheless something seems off.

• So you convene with me that nothing seems to be wrong? The script by the way worked fine until I upgraded to 11.0, that's all that was different. – SuperCiocia Jan 26 '17 at 16:33
• @SuperCiocia Well, something is definitely wrong, since the plots don't work. But it may not be your data: my cursory check seems to indicate that Plot is behaving funny here. This is also supported by the fact that it used to work before. I'd submit a support request to Wolfram Support; perhaps they will recognize a known problem. If you do so, though, it would be best if you could find a simpler, minimal example that highlights this behavior. – MarcoB Jan 26 '17 at 16:40