Switching between solutions of a fourth order polynomial

Question

I am plotting the roots of a fourth-order polynomial as a function of some parameter $r$ with the following code:

mat = {{\[Lambda]^2 - 2 \[Lambda]/r + 3/(4 r^2) + 1, 2}, 
       {2, \[Lambda]^2 + 2 \[Lambda]/r + 3/(4 r^2) + 3}};
sol1 = \[Lambda] /. Solve[Det[mat]== 0, \[Lambda]][[1]];
sol2 = \[Lambda] /. Solve[Det[mat]== 0, \[Lambda]][[2]];
sol3 = \[Lambda] /. Solve[Det[mat]== 0, \[Lambda]][[3]];
sol4 = \[Lambda] /. Solve[Det[mat]== 0, \[Lambda]][[4]];

Plot[{Im[sol1], Im[sol2], Im[sol3], Im[sol4]}, {r, .1, 2}, 
     Frame -> True, PlotStyle -> styles, Exclusions -> None]
Plot[{Re[sol1], Re[sol2], Re[sol3], Re[sol4]}, {r, .1, 2}, 
     Frame -> True, PlotStyle -> styles, Exclusions -> None]

The plot of the real part looks fine, but for the imaginary part it looks like Mathematica is switching between a solution and its complex conjugate for different values of the parameter $r$. This is confirmed by plotting the absolute value of the imaginary part of the four different solutions.

Hence, my question is: is there a simple way to make sure that Mathematica is plotting the same solution for every value of $r$, without switching between conjugate solutions?

Answer 1

Mathematica's "canonical" orderings of symbolic expressions do not always correspond to numerical orderings¹, which is (I think) why you're running into problems here. One way to get around this is to tell Mathematica to calculate the roots numerically instead:

mat = {{\[Lambda]^2 - 2 \[Lambda]/r + 3/(4 r^2) + 1, 2}, 
       {2, \[Lambda]^2 + 2 \[Lambda]/r + 3/(4 r^2) + 3}};
sol1 := \[Lambda] /. NSolve[Det[mat] == 0, \[Lambda]][[1]];
sol2 := \[Lambda] /. NSolve[Det[mat] == 0, \[Lambda]][[2]];
sol3 := \[Lambda] /. NSolve[Det[mat] == 0, \[Lambda]][[3]];
sol4 := \[Lambda] /. NSolve[Det[mat] == 0, \[Lambda]][[4]];

Plot[{Im[sol1], Im[sol2], Im[sol3], Im[sol4]}, {r, .1, 2}, 
      Frame -> True, PlotStyle -> styles, Exclusions -> None]
Plot[{Re[sol1], Re[sol2], Re[sol3], Re[sol4]}, {r, .1, 2}, 
      Frame -> True, PlotStyle -> styles, Exclusions -> None]

Note the use of SetDelayed (aka :=) in the definitions of the solutions; if you use a standard Set (=), Mathematica will complain that you haven't defined r and that it can't solve it numerically.

This does take more time to run than the original code. You could probably speed things up somewhat by defining something like

sols := \[Lambda] /. NSolve[Det[mat] == 0, \[Lambda]]

as a list, and then plotting Im[sols] and Re[sols]. That way you're only calling NSolve once for each value of r, rather than four times.

---

¹ As an example:

Sort[{Sqrt[2], 2}]
(* {2, Sqrt[2]} *)    

Answer 2

Correction

I correct my error in relation to range r. I post this as a way of visualizing the roots of polynomial (with this correction).

mat[\[Lambda]_, 
   r_] := {{\[Lambda]^2 - 2 \[Lambda]/r + 3/(4 r^2) + 1, 
    2}, {2, \[Lambda]^2 + 2 \[Lambda]/r + 3/(4 r^2) + 3}};
func[\[Lambda]_, r_] := Det[mat[\[Lambda], r]]

Visualizing:

Manipulate[With[{res = x /. Solve[func[x, r] == 0, x]},
  max = 1.2 Max[Abs /@ res];
  pts = Through[{Re, Im}[#]] & /@ res;
  Row[{Plot[func[x, r], {x, -max, max}, MeshFunctions -> (#2 &), 
     Mesh -> {{0.}}, MeshStyle -> Directive[Red, PointSize[0.04]], 
     ImageSize -> 200],
    Show[ContourPlot[
      Abs[func[x + I y, r]], {x, -max, max}, {y, -max, max}, 
      Contours -> (# Abs[func[max/2 + max I/2, r]] & /@ {0.1, 0.05, 
          0.01, 0.001}), 
      ColorFunction -> (ColorData["Pastel"][#^(1/4)] &), Axes -> True,
       AxesStyle -> White], 
     Graphics[{Red, Text[Style["\[Times]", 20, Bold], #] & /@ pts}], 
     ImageSize -> 200]}]], {r, 0.1, 2}]