How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters) - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-09-22T23:43:00Z https://mathematica.stackexchange.com/feeds/question/39747 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/39747 7 How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters) Steve https://mathematica.stackexchange.com/users/7479 2014-01-03T17:09:24Z 2016-03-08T18:00:04Z <p>I've been working on this research problem for months with some success by no global/proper solution to the problem.</p> <p>So I have square, Hermitian matrices of various sizes ($8,20,38$ dimensions) with functional dependence on two variables, call them <code>x</code> and <code>y</code>. I am interested in plotting all the eigenvalues as functions of the $2$ variables <code>x</code> and <code>y</code>. </p> <p>Let's call <code>M</code> the $8 \times 8$ hermitian matrix.</p> <pre><code>Plot3D[ Eigenvalues[M], {x,-1,1}, {y,-1,1}] </code></pre> <p>gives $8$ sheets corresponding to the eigenvalues of <code>M</code> as functions of <code>x</code> and <code>y</code>. </p> <p>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, <em>Mathematica</em> organizes the $8$ values differently, so the different sheets swap places, so we get tons of random garbage vertical bits between sheets. </p> <p>I have spent months trying to resolve this in different ways, mostly involving methods to force <em>Mathematica</em> into a consistent sorting method. Here are some examples of failed methods: </p> <pre><code>Plot3D[ Sort[Eigenvalues[M]], {x,-1,1}, {y,-1,1}] Plot3D[ Sort[Eigenvalues[M], Abs[#1] &lt;= Abs[#2]&amp;], {x,-1,1}, {y,-1,1}] Block[{set=Ordering[ Eigenvalues[ M/.x-&gt;0 /. y-&gt;0], All, Abs[#1] &lt;= Abs[#2]&amp;]}, Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]] Block[{set=Ordering[ Ordering[Eigenvalues[ M/.x-&gt;0 /.y-&gt;0], All, Abs[#1] &lt;= Abs[#2]&amp;]]}, Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]] Block[{set = Ordering[Eigenvalues[ M /.x-&gt;0 /.y-&gt;0], All, Abs[#1] &lt;= Abs[#2]&amp;]}, Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]] Block[{set = Ordering[ Ordering[ Eigenvalues[ M /.x-&gt;0 /.y-&gt;0], All, Abs[#1] &lt;= Abs[#2]&amp;]]}, Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]] </code></pre> <p>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. </p> <p>Please help me understand the proper solution to this problem!</p> https://mathematica.stackexchange.com/questions/39747/-/39754#39754 11 Answer by george2079 for How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters) george2079 https://mathematica.stackexchange.com/users/2079 2014-01-03T19:37:50Z 2015-05-29T16:02:19Z <blockquote> <p>note see bottom of answer for a required fix for version 10+</p> <p>see <a href="https://mathematica.stackexchange.com/a/84669/2079">here</a> for an alternate approach</p> </blockquote> <p>A 2D example of what i suggested in my comment:</p> <p>Based on Belisarius example..we get three lines obviously not connected properly:</p> <pre><code> m = SparseArray[{i_, j_} -&gt; 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 -&gt; True, PlotRange -&gt; All, PlotStyle -&gt; Hue[First@#2/3]] &amp;, Range]] </code></pre> <p><img src="https://i.stack.imgur.com/xWZzD.png" alt="enter image description here"></p> <p>Now construct a sort of surface curvature penalty function and minimize by selecting from each of the possibilities..</p> <pre><code>v = Table[V[i], {i, Length[alle]}]; cons = Table[1 &lt;= V[i] &lt;= 3, {i, Length[alle]}]; Off[Part::pspec] SetOptions[NMinimize, MaxIterations -&gt; 500] v1 = v /. Last@ Minimize[ {Total[((#[] + #[] - 2 #[])^2 &amp; /@ Partition[MapIndexed[#[[v[[First@#2]]]] &amp;, alle] , 3, 1, {1, 3}, {}])], cons}, v, Integers] </code></pre> <p>Nicely finds a smooth line..</p> <p><img src="https://i.stack.imgur.com/wbLM7.png" alt="enter image description here"></p> <pre><code>Show[original, ListPlot[MapIndexed[#[[v1[[First@#2]]]] &amp;, alle], Joined -&gt; True, PlotStyle -&gt; {Thick, Black, Dashed}]] </code></pre> <p>Repeat the process adding to the constraints like this:</p> <pre><code>cons = Join[cons, Table[V[i] != v1[[i]], {i, Length[alle]}] ] </code></pre> <p><img src="https://i.stack.imgur.com/UlXH5.png" alt="enter image description here"></p> <p>Works well even if the input ordering is a random mess:</p> <p><img src="https://i.stack.imgur.com/yAPfg.png" alt="enter image description here"></p> <p>Quite likely horribly slow in 3d...though you might tackle it as a series of 2d sorts and join them together.</p> <p>See here for that <code>SetOptions[NMinimise]</code> trick..</p> <p><a href="http://forums.wolfram.com/mathgroup/archive/2006/Aug/msg00179.html" rel="nofollow noreferrer">http://forums.wolfram.com/mathgroup/archive/2006/Aug/msg00179.html</a> </p> <p>Edit - a more compact version, reordering the points instead of adding more constraints at each step: </p> <pre><code> v = Table[V[i], {i, Length[alle]}]; Off[Part::pspec]; SetOptions[NMinimize, MaxIterations -&gt; 500]; Do[ alle = MapThread[ Prepend[ Drop[ #2 , {#1}] , #2[[#1]] ] &amp;, { (v /. Last@Minimize[{Total[((#[] + #[] - 2 #[])^2 &amp; /@ Partition[MapIndexed[#[[v[[First@#2]]]] &amp;, alle], 3, 1, {1, 3}, {}])], Table[k &lt;= V[i] &lt;= 3, {i, Length[alle]}]}, v, Integers]) , alle}] , {k, 2}]; Show[original, ListPlot[Flatten@Take[alle, All, {#}], Joined -&gt; True, PlotStyle -&gt; {Thick, Black, Dashed}] &amp; /@ Range] </code></pre> <h2>v10+ fix</h2> <p>for mathematica versions with the <code>Indexed</code> function, replace :</p> <pre><code> MapIndexed[#[[v[[First@#2]]]] &amp;, alle] </code></pre> <p>with </p> <pre><code> Indexed @@@ Transpose[{alle, v}] </code></pre> <p>(and setting <code>Off[Part::pspec];</code> is no longer needed )</p> https://mathematica.stackexchange.com/questions/39747/-/39756#39756 4 Answer by Steve for How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters) Steve https://mathematica.stackexchange.com/users/7479 2014-01-03T20:05:03Z 2014-03-03T20:39:51Z <p>Thanks to everyone for the help. Here's my working solution:</p> <p>======<br> 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</p> <pre><code>Sort[(*function*),Abs[#1]&lt;=Abs[#2]&amp;] </code></pre> <p>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. </p> <p>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<br> ====== </p> <p>"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</p> <pre><code>lowb=Min[Drop[Ordering[ Sort[Eigenvalues[M/. x -&gt; 0 /. y -&gt; 0]], All, Abs[#1] &lt;= Abs[#2] &amp;], -(Dimensions[M][] - lownumber)]]]; lowt=Max[Drop[Ordering[ Sort[Eigenvalues[M/. x -&gt; 0 /. y -&gt; 0]], All, Abs[#1] &lt;= Abs[#2] &amp;], -(Dimensions[M][] - lownumber)]]]; </code></pre> <p>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: </p> <pre><code>Plot3D[ Table[Sort[Eigenvalues[hamiltonian]][[s]], {s, lowb, lowt}] </code></pre> <p>or </p> <pre><code>Plot3D[ Table[RankedMin[Eigenvalues[hamiltonian],s], {s, lowb, lowt}] </code></pre> https://mathematica.stackexchange.com/questions/39747/-/109444#109444 4 Answer by Michael for How to plot several functions without jumping? (multiple eigenvalues of a system as functions of 2 parameters) Michael https://mathematica.stackexchange.com/users/35947 2016-03-08T17:33:34Z 2016-03-08T18:00:04Z <p>I know the question was answered long ago but I ran into the same problem and was going to implement the solution above when I may have stumbled on a much simpler solution. I would run the following code for a 16x16 matrix (which is block diagonal with the largest blocks being 4x4)</p> <pre><code>Plot[Eigenvalues[FullHamiltonian], {x, -1, 1}] </code></pre> <p>and would obtain the horrible graph</p> <p><a href="https://i.stack.imgur.com/fX5Sb.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/fX5Sb.jpg" alt="A lot of jumps between branches."></a></p> <p>I then realised that when I tried to plot a single eigenvalue I could use the horrible <code>Root["polynomial"&amp;,1]</code> expression for a given eigenvalue to obtain just one of these lines in a smooth way. I therefore assumed it might be possible to solve it using</p> <pre><code>EV=Eigenvalues[FullHamiltonian] Plot[EV, {x,-1,1}] </code></pre> <p>And in fact this gives me the following plot</p> <p><a href="https://i.stack.imgur.com/6Z1PR.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/6Z1PR.jpg" alt="The jumps are all gone."></a></p> <p>I repeated the experiment using two variables (on a 4x4 matrix) and it still seems to work. <a href="https://i.stack.imgur.com/DaTfn.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DaTfn.jpg" alt="Jumps in 3D."></a><a href="https://i.stack.imgur.com/TIw4b.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/TIw4b.jpg" alt="Smooth in 3D."></a></p> <p>My understanding is that the <code>Root</code> function has a well defined ordering whereas the <code>Eigensystem</code> function doesn't if all coefficients are numerical. When we do the diagonalisation analytically first, we are guaranteed to be using <code>Root</code> with the same argument every time. This should probably be slower than doing the diagonalisation of the fully numerical matrix but has the advantage of making the code much more legible.</p> <p>As I have only tested it on an effectively 4x4 matrix it would be interesting to see what happens for larger matrices but I cannot see any obvious reason why it would fail.</p> <p>EDIT: I tested on a 5x5 (first one without analytic solutions) and it still seems to work, however computing the eigenvectors analytically is very slow. There might be a way to combine the analytic eigenvalues, which are fast, with a fully numeric finding of the eigenvectors to combine speed and code sobriety.</p>