# How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters)

I've been working on this research problem for months with some success by no global/proper solution to the problem.

So I have square, Hermitian matrices of various sizes ($8,20,38$ dimensions) with functional dependence on two variables, call them x and y. I am interested in plotting all the eigenvalues as functions of the $2$ variables x and y.

Let's call M the $8 \times 8$ hermitian matrix.

Plot3D[ Eigenvalues[M], {x,-1,1}, {y,-1,1}]


gives $8$ sheets corresponding to the eigenvalues of M as functions of x and y.

However, there will be many "jumps" between the sheets which totally obscure the plots. Basically, I think that during the plotting of the eigenvalues, for certain regions of the parameter space, Mathematica organizes the $8$ values differently, so the different sheets swap places, so we get tons of random garbage vertical bits between sheets.

I have spent months trying to resolve this in different ways, mostly involving methods to force Mathematica into a consistent sorting method. Here are some examples of failed methods:

Plot3D[ Sort[Eigenvalues[M]], {x,-1,1}, {y,-1,1}]

Plot3D[ Sort[Eigenvalues[M], Abs[#1] <= Abs[#2]&], {x,-1,1}, {y,-1,1}]

Block[{set=Ordering[ Eigenvalues[ M/.x->0 /. y->0], All, Abs[#1] <= Abs[#2]&]},
Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]]

Block[{set=Ordering[ Ordering[Eigenvalues[ M/.x->0 /.y->0], All,
Abs[#1] <= Abs[#2]&]]},
Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]]

Block[{set = Ordering[Eigenvalues[ M /.x->0 /.y->0], All, Abs[#1] <= Abs[#2]&]},
Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]]

Block[{set = Ordering[ Ordering[ Eigenvalues[ M /.x->0 /.y->0], All,
Abs[#1] <= Abs[#2]&]]},
Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]]


as well as other permutations of similar ideas and more. Some of these methods have worked somewhat, reducing the jumping or working in certain regimes of the parameter domain. But nothing works consistently.

Please help me understand the proper solution to this problem!

-
do you have an actual example? – Nasser Jan 3 '14 at 17:14
Not good: i.imgur.com/aSjHES5.png?1 – Steve Jan 3 '14 at 17:21
– Steve Jan 3 '14 at 17:22
Example: m = SparseArray[{i_, j_} -> Sin[i j x y], {3, 3}]; Plot3D[ Eigenvalues[m][[1]], {x, 0, 1}, {y, 0, 1}, PlotRange -> All] – belisarius has settled Jan 3 '14 at 18:13
Try generating a Table, then working on an algorithm to sort the table, finally using ListPlot3D.. – george2079 Jan 3 '14 at 18:20

note see bottom of answer for a required fix for version 10+

see here for an alternate approach

A 2D example of what i suggested in my comment:

Based on Belisarius example..we get three lines obviously not connected properly:

 m = SparseArray[{i_, j_} -> Sin[i j 9 /10 y], {3, 3}];
alle = Table[Eigenvalues[m], {y, 0, 1, .1}];
original = Show[
MapIndexed[ListPlot[Flatten[Take[alle, All, {#}] , 2],
Joined -> True, PlotRange -> All,
PlotStyle -> Hue[First@#2/3]] &, Range[3]]]


Now construct a sort of surface curvature penalty function and minimize by selecting from each of the possibilities..

v = Table[V[i], {i, Length[alle]}];
cons  = Table[1 <= V[i] <= 3, {i, Length[alle]}];
Off[Part::pspec]
SetOptions[NMinimize, MaxIterations -> 500]
v1 = v /. Last@
Minimize[ {Total[((#[[1]] + #[[3]] - 2 #[[2]])^2 & /@
Partition[MapIndexed[#[[v[[First@#2]]]] &, alle] , 3,
1, {1, 3}, {}])], cons}, v, Integers]


Nicely finds a smooth line..

Show[original,
ListPlot[MapIndexed[#[[v1[[First@#2]]]] &, alle], Joined -> True,
PlotStyle -> {Thick, Black, Dashed}]]


Repeat the process adding to the constraints like this:

cons  = Join[cons, Table[V[i] != v1[[i]], {i, Length[alle]}] ]


Works well even if the input ordering is a random mess:

Quite likely horribly slow in 3d...though you might tackle it as a series of 2d sorts and join them together.

See here for that SetOptions[NMinimise] trick..

Edit - a more compact version, reordering the points instead of adding more constraints at each step:

 v = Table[V[i], {i, Length[alle]}];
Off[Part::pspec];
SetOptions[NMinimize, MaxIterations -> 500];
Do[
alle = MapThread[ Prepend[  Drop[ #2 , {#1}]  , #2[[#1]] ]  &,
{ (v /. Last@Minimize[{Total[((#[[1]] + #[[3]] - 2 #[[2]])^2 & /@
Partition[MapIndexed[#[[v[[First@#2]]]] &, alle], 3, 1, {1, 3}, {}])],
Table[k <= V[i] <= 3, {i, Length[alle]}]}, v, Integers]) ,
alle}] , {k, 2}];
Show[original,
ListPlot[Flatten@Take[alle, All, {#}], Joined -> True,
PlotStyle -> {Thick, Black, Dashed}] & /@ Range[3]]


## v10+ fix

for mathematica versions with the Indexed function, replace :

 MapIndexed[#[[v[[First@#2]]]] &, alle]


with

 Indexed @@@ Transpose[{alle, v}]


(and setting Off[Part::pspec]; is no longer needed )

-
Haven't see it before. Very nice – belisarius has settled Mar 3 '14 at 20:38
Sounds great. However, with mathematica v10, I get the following error: Part::pkspec1: The expression V[1] cannot be used as a part specification. >> Is there something wrong? after the instuction v1 = v /. Last@ Minimize[ ..... – altroware Aug 13 '14 at 13:52
i guess the warning i surpressed got elevated to an error. i dont have 10 to work on it though – george2079 Aug 13 '14 at 14:19
@altroware fixed that issue i think - verified only with 9 which it turns out has the Indexed function (undocumented). – george2079 Aug 19 '14 at 17:57
@george2079 would you care to help me out with my question where I try to use your above answer (the questions are strongly related): mathematica.stackexchange.com/questions/83906/… – Mr Z May 27 at 8:01

Thanks to everyone for the help. Here's my working solution:

======
First, the spirit of the underlying problem and the nature of the ultimate solution. Simply sorting by the function values (executed with the default one-argument Sort command) is enough to separate the functions and prevent the jumping. However, it does not give the user control over the functions and it is difficult to do anything further (e.g. remove several of the functions from the plot, apply colors in a predictable way, etc.). If you sort instead by the two argument Sort with a pure absolute value function as

Sort[(*function*),Abs[#1]<=Abs[#2]&]


there is a pairwise degeneracy problem when two functions cross planes equidistant above and below the zero plane, and you get tons of pairwise swapping at those places and a garbage plot. The spirit of this attempt, however, would be to gain control over the smallest absolute value function bands, for example.

The solution that follows uses a different technique. We instead use a reference point (the origin) to compute a single list of eigenvalues, perform simple operations there to compute which bands in the final plot will correspond to meeting a desired criteria (such as the smallest n bands in terms of absolute value) and then plot the functions sorted by a simple Sort[] with these indices
======

"lownumber" is the desired number of plotted eigenvalues, the smallest lownumber bands in terms of absolute value. lowb and lowt are the lower and upper bound indices for the plotted functions

lowb=Min[Drop[Ordering[
Sort[Eigenvalues[M/. x -> 0 /. y -> 0]], All,
Abs[#1] <= Abs[#2] &], -(Dimensions[M][[1]] - lownumber)]]];

lowt=Max[Drop[Ordering[
Sort[Eigenvalues[M/. x -> 0 /. y -> 0]], All,
Abs[#1] <= Abs[#2] &], -(Dimensions[M][[1]] - lownumber)]]];


After this, there are two possible solutions. I have tested them both to be of the same speed to within about 3%, for matrices up to size 38. For much larger matrices, I am not sure, one might be preferable. Here they are:

Plot3D[
Table[Sort[Eigenvalues[hamiltonian]][[s]], {s, lowb, lowt}]


or

Plot3D[
Table[RankedMin[Eigenvalues[hamiltonian],s], {s, lowb, lowt}]

-