Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
1  
do you have an actual example? –  Nasser Jan 3 at 17:14
    
Not good: i.imgur.com/aSjHES5.png?1 –  Steve Jan 3 at 17:21
    
good: i.imgur.com/vWeSYGc.png?1 –  Steve Jan 3 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 Jan 3 at 18:13
    
Try generating a Table, then working on an algorithm to sort the table, finally using ListPlot3D.. –  george2079 Jan 3 at 18:20
show 1 more comment

2 Answers 2

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

enter image description here

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

enter image description here

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]}] ]

enter image description here

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

enter image description here

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

http://forums.wolfram.com/mathgroup/archive/2006/Aug/msg00179.html

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]]
share|improve this answer
    
Haven't see it before. Very nice –  belisarius Mar 3 at 20:38
add comment
up vote 2 down vote accepted

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}] 
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.