4
$\begingroup$

I am trying to find the solution of cubic equation as a function of $kz$

-1.05976 + λ (-4.3872 + (-3.9 - 
       1. λ) λ) + (1.10144 + 0.624 λ) Cos[
   1.31 kz] + 0.035152 Cos[2.62 kz]==0
{e1[kz_], e2[kz_], 
  e3[kz_]} = λ /. 
   NSolve[-1.05976 + λ (-4.3872 + (-3.9 - 
       1. λ) λ) + (1.10144 + 0.624 λ) Cos[
   1.31 kz] + 0.035152 Cos[2.62 kz]==0, λ] // FullSimplify

If I plot e1[kz] as function of kz,

Plot[e1[z], {z, 0, π/1.31}]

enter image description here

I get a function that is not smooth. However, if I plot another solution alongside then it seems two solution got mixed up.

Plot[{e1[z], e3[z]}, {z, 0, π/1.31}]

enter image description here

Does anyone know why is this the case? Also, can you fix this?

$\endgroup$
3
  • $\begingroup$ It is not a bug, it is a feature: blue line is the lower root, orange live is the next one. If you do not like this numerical order, it is possible to resort according to your preferences, which ones? $\endgroup$
    – yarchik
    Commented Sep 6, 2021 at 22:04
  • $\begingroup$ @yarchik I get what you are saying and make sense. But what if I want a solution that is smooth regardless of its value being higher or lower. Is this possible? $\endgroup$
    – sslucifer
    Commented Sep 6, 2021 at 22:31
  • 2
    $\begingroup$ The issue we encounter here is well known on the basis of Root objects, which appear behind the scene of NSolve, see e.g. find where 3 inequalities are simultaneously greater than zero but we could find many more posts here regarding this issue. For this reason I recommend changing the title of the question since this might lead to confusion, a better example : "Aren't roots of polynomial equations smooth functions of their parameters?" $\endgroup$
    – Artes
    Commented Sep 7, 2021 at 23:13

2 Answers 2

3
$\begingroup$
Clear["Global`*"]

The equation can be solved exactly

eqn = -1.05976 + λ (-4.3872 + (-3.9 - 
          1. λ) λ) + (1.10144 + 0.624 λ) Cos[
      1.31 kz] + 0.035152 Cos[2.62 kz] == 0;

sol[kz_] = 
  Solve[eqn // Rationalize // Simplify, λ, Reals] // Simplify;

EDIT: For the revised plot range,

Plot[Evaluate[λ /. sol[kz]], {kz, 0, Pi/1.31},
 Frame -> True,
 FrameLabel -> (Style[#, 12, Bold] & /@ {kz, λ}),
 PlotLegends -> Placed[{λ1, λ2, λ3}, {.25, .5}]]

enter image description here

Looking at only the first two solutions

Plot[Evaluate[Most[λ /. sol[kz]]], {kz, 0, Pi/1.31},
 Frame -> True,
 FrameLabel -> (Style[#, 12, Bold] & /@ {kz, λ}),
 PlotLegends -> Placed[{λ1, λ2}, {.25, .6}]]

enter image description here

From the condition in the ConditionalExpression of the second root, the switchover occurs at

{so = Reduce[{sol[kz][[2, 1, -1, -1, 2]], 1 < kz < 2}, kz][[1, -1]], so // N}

(* {200/131 ArcTan[(33 Sqrt[3/23])/7], 1.58739} *)

EDIT: Alternatively, to find the switchover use

{so = kz /. 
   Solve[{Equal @@ (λ /. Most[sol[kz] // Normal]), 1 < kz < 2}, 
     kz][[1]], so // N}

(* {200/131 ArcTan[(33 Sqrt[3/23])/7], 1.58739} *)

Using Piecewise, the smooth curves are then

Format[lambda[n_]] = Subscript[λ, n];

lambda[1][kz_] = Piecewise[{
    {λ /. sol[kz][[1]], kz <= so},
    {λ /. sol[kz][[2]], so < kz < Pi/1.31}}];

lambda[2][kz_] = Piecewise[{
    {λ /. sol[kz][[2]], kz <= so},
    {λ /. sol[kz][[1]], so < kz < Pi/1.31}}];

Plot[{lambda[1][kz], lambda[2][kz]}, {kz, 0, Pi/1.31},
 Frame -> True,
 FrameLabel -> (Style[#, 12, Bold] & /@ {kz, λ}),
 PlotLegends -> Placed[{lambda[1], lambda[2]}, {.25, .6}]]

enter image description here

$\endgroup$
3
  • $\begingroup$ When the plots extend out farther to the right, the problem remains, no? $\endgroup$
    – LouisB
    Commented Sep 7, 2021 at 3:55
  • $\begingroup$ Yes. I used the plot range indicated in the OP's code, i.e., {kz, 0, Pi/3}. The smooth curves switch between branches of the solution. This is caused by meeting the conditions of the ConditionalExpressions. $\endgroup$
    – Bob Hanlon
    Commented Sep 7, 2021 at 4:52
  • $\begingroup$ @BobHanlon I made a mistake, the range of kz is from 0 to Pi/1.31. So LouisB is correct if you extend the plots, the two solutions switches. I made the necessary edits in the question. $\endgroup$
    – sslucifer
    Commented Sep 7, 2021 at 7:53
6
$\begingroup$

With a tip of the hat to @MichaelE2 (from whom I learned this trick), we can differentiate the implicit relation between $\lambda$ and $k_z$ and then solve it as an ODE:

neweq = eq /. λ -> λ[kz];
startkz = π/6;
startval = λ[startkz] /. Solve[neweq /. kz -> startkz, λ[startkz]]
soln = NDSolve[{D[neweq, kz], λ[startkz] == startval}, λ, {kz, 0, π/1.31}, Method -> "StiffnessSwitching"]
Plot[Evaluate[Table[Indexed[λ[kz] /. First[soln], i], {i, 1, 3}]], {kz, 0, Pi/1.31}]

enter image description here

This method may not be entirely robust, but it does seem to work in this case. Notes:

  • The initial values for the ODE are specified at startkz = Pi/6 in order to ensure that all three roots are distinct at that point and therefore that Mathematica produces three different solutions.
  • The StiffnessSwitching option is necessary to integrate past the points where the roots cross.
  • The use of Table and Indexed is required because Mathematica is returning a vector-valued InterpolatingFunction, and some chicanery is needed to get Plot to plot the three components with distinct colors (among other things you might want to do with the result.)
$\endgroup$

Your Answer

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

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