# Switching between solutions of a fourth order polynomial

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?

Mathematica's "canonical" orderings of symbolic expressions do not always correspond to numerical orderings1, 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.

1 As an example:

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

• Exactly what I was looking for, many thanks! Just one small question: I have never seen these nice colours you have in your plot. They are different than the usual Blue/Red/Orange, right? How did you get them? Jun 16, 2015 at 15:50
• @Funzies They're the new V10 colors. Jun 16, 2015 at 15:56
• @MichaelSeifert very nice and +1 of course...I only posted my answer out of guilt for misreading the r range...:) Jun 17, 2015 at 7:08

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


• in the region r in [0.1,0.2] the solutions are indeed real, however for larger r the solutions become imaginary. Jun 16, 2015 at 12:01
• @Funzies I only based this on your own codes range for r...hence my post...thank you for clarication Jun 16, 2015 at 12:03
• @Funzies mea culpa I mistyped the range! My sincere apologies. Jun 16, 2015 at 12:06
• @Funzies I have corrected my misconception. It may be not what you want but perhaps it may be helpful. Again apologies for my error. Jun 16, 2015 at 12:52