4
$\begingroup$

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?

$\endgroup$
7
$\begingroup$

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]

enter image description here

enter image description here

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]} *)    
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Funzies Jun 16 '15 at 15:50
  • $\begingroup$ @Funzies They're the new V10 colors. $\endgroup$ – Michael E2 Jun 16 '15 at 15:56
  • $\begingroup$ @MichaelSeifert very nice and +1 of course...I only posted my answer out of guilt for misreading the r range...:) $\endgroup$ – ubpdqn Jun 17 '15 at 7:08
1
$\begingroup$

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

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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