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$
2

2 Answers 2

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$
3
  • $\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, 2015 at 15:50
  • $\begingroup$ @Funzies They're the new V10 colors. $\endgroup$
    – Michael E2
    Jun 16, 2015 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, 2015 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$
4
  • $\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, 2015 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, 2015 at 12:03
  • $\begingroup$ @Funzies mea culpa I mistyped the range! My sincere apologies. $\endgroup$
    – ubpdqn
    Jun 16, 2015 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, 2015 at 12:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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