1
$\begingroup$

I am trying to find the roots of a polynomial of eight degree $$\sin^8\varphi-2\lambda_i^2\sin^6\varphi+(6\lambda_i^2+\lambda_i^4)\sin^4\varphi-(4\lambda_i^2+2\lambda_i^4)\sin^2\varphi+\lambda_i^4=0$$ where $$\lambda _i=\frac{mg-F\frac i n}{2kl}.$$

I have to find the solutions to that eight degree polynomial for each $\lambda_i$ where $i=1...n$, which I did with the following code:

F = 1430;(*N*)
k = 1010;(*N/m*)
l = 0.4;(*m*)
m = 10;(*kg*)
g = 9.81;(*m/s^2*)
n = 100;

λ = Table[(m g - F*(i/n))/(2 k l), {i, 0, n, 1}];
resitve = Table[
    Solve[
     x^8 - 2 λ[[i]]^2 x^6 + (6 λ[[i]]^2 + λ[[i]]^4) x^4 - (2 λ[[i]]^4 + 
     4 λ[[i]]^2) x^2 + λ[[i]]^4 == 0, x, Reals
    ], 
   {i, 1, Length[λ], 1}
  ]

kotiMinusBig =   Table[ArcSin[resitve[[i, 1, 1, 2]]]*180/Pi, {i, 1, Length[λ], 1}];
kotiMinusSmall = Table[ArcSin[resitve[[i, 2, 1, 2]]]*180/Pi, {i, 1, Length[λ], 1}];
kotiPlusSmall =  Table[ArcSin[resitve[[i, 3, 1, 2]]]*180/Pi, {i, 1, Length[λ], 1}];
kotiPlusBig =    Table[ArcSin[resitve[[i, 4, 1, 2]]]*180/Pi, {i, 1, Length[λ], 1}];

tockeMinusBig =   Table[{(F/n)*(i - 1)/1000, kotiMinusBig[[i]]}, {i, 1, Length[λ], 1}];
tockeMinusSmall = Table[{(F/n)*(i - 1)/1000, kotiMinusSmall[[i]]}, {i, 1, Length[λ], 1}];
tockePlusSmall =  Table[{(F/n)*(i - 1)/1000, kotiPlusSmall[[i]]}, {i, 1, Length[λ], 1}];
tockePlusBig =    Table[{(F/n)*(i - 1)/1000, kotiPlusBig[[i]]}, {i, 1, Length[λ], 1}];

ListPlot[{tockeMinusBig, tockeMinusSmall, tockePlusBig, tockePlusSmall}]

This produces the following output: enter image description here

Now my question is: Does Mathematica somehow automatically order the solutions to Solve?

Or a question related to the example above: Starting at the left: Does the $y$ value of those green dots really start increasing for $x>0.1$, or do they actually continue to decrease continuously to the blue dots for $x>0.1$, but Mathematica ordered them weirdly and now it looks as if the green solutions starts increasing?

Which is the case, and how can I correct it?

NOTE: The problem explicitly says that the $F_i$ domain has to stay discrete.

$\endgroup$
6
  • 1
    $\begingroup$ It does, see e.g. this answer Finding parameters making real part of eigenvalues vanish. There are also different posts on this issue, e.g. Switching between solutions of a fourth order polynomial $\endgroup$
    – Artes
    Mar 14, 2016 at 15:39
  • $\begingroup$ You cn see them here resitve = Solve[x^8 - 2 l^2 x^6 + (2 l^2 + l^4 + 4 l^2) x^4 - (2 l^4 + 4 l^2) x^2 + l^4 == 0, x, Reals]; Plot[x /. # & /@ resitve, {l, 0, 10}, Evaluated -> True] $\endgroup$ Mar 14, 2016 at 15:41
  • 1
    $\begingroup$ @skrat You may have a typo in one of your expressions: the coefficients of the $\sin^2\varphi$ expressions in your formatted equation at the top and in your code do not seem to be the same ($4\lambda_i^2-2\lambda_i^4$ vs. (2 λ^4 + 4 λ^2)). Which one is correct? $\endgroup$
    – MarcoB
    Mar 14, 2016 at 16:17
  • $\begingroup$ @MarcoB: you are right. There was a typo. I have edited my OP now. The code was right. There has to be a "+" sign in the brackets. $\endgroup$
    – skrat
    Mar 14, 2016 at 16:25
  • $\begingroup$ @Dr.belisarius: You might have given me the answer, but due to my poor knowledge of mathematica I have no idea what you are trying to show me. :D $\endgroup$
    – skrat
    Mar 14, 2016 at 16:42

1 Answer 1

2
$\begingroup$

The numerical roots are sorted in increasing order.

For this problem a cleaner approach is to use Solve to get an analytic expression and plot that:

asol = x /. 
   Solve[x^8 - 
      2 \[Lambda]i^2 x^6 + (6 \[Lambda]i^2 + \[Lambda]i^4) x^4 -
     (2 \[Lambda]i^4 + 4 \[Lambda]i^2) x^2 + \[Lambda]i^4 == 0, x];
Show[{Plot[ArcSin[{asol[[1]], asol[[8]]}] /Degree /. 
        \[Lambda]i -> m g/(2 k l) - x,
        {x, 0, F/(2 k l)}, PlotStyle -> Red, PlotRange -> All], 
     Plot[ArcSin[{asol[[2]], asol[[7]]}] /Degree /.
        \[Lambda]i ->  m g/(2 k l) - x, 
        {x, 0, F/(2 k l)}, PlotStyle -> Blue, PlotRange -> All], 
     Plot[ArcSin[{asol[[3]], asol[[6]]}]/Degree /. 
        \[Lambda]i ->  m g/(2 k l) - x, 
        {x, 0, F/(2 k l)}, PlotStyle -> Green, PlotRange -> All], 
     Plot[ArcSin[{asol[[4]], asol[[5]] }]/Degree /.
        \[Lambda]i -> m g/(2 k l) - x, 
        {x, 0, F/(2 k l)},PlotStyle -> Black, PlotRange -> All]}]

enter image description here

Note there are 8 solutions each of which is real on half of the plot and I had to manually sort them so they look continuous. (This is much easier than sorting your list data however )

Edit: an improved version that splices the pairs of solutions together so it plots smoothly at the transition.

 f[x_, soli_, solj_] := Module[{res}, ArcSin[Which[
     Im[res = asol[[soli]] /. \[Lambda]i -> m g/(2 k l) - x] == 0, 
     res ,
     Im[res = asol[[solj]] /. \[Lambda]i -> m g/(2 k l) - x] == 0, 
     res ]]/Degree]
 Plot[{f[x, 1, 8], f[x, 2, 7], f[x, 3, 6], f[x, 4, 5]},
      {x, 0, F/(2 k l)}]

Edit: a discrete ListPlot version, just showing one curve for clarity.

data1 = Flatten[Select[#, Im[#] == 0 &] & /@
    ((ArcSin[{asol[[1]], asol[[8]]}] /Degree
            /. \[Lambda]i -> #) & /@ \[Lambda])];
Length@data1 == Length@\[Lambda]  (* True this will fail if both
        soln's are real or both complex at the same point *)
ListPlot[MapIndexed[{F/n (First@#2 - 1)/1000, #} &, data1], 
 PlotRange -> All]

enter image description here

$\endgroup$
6
  • $\begingroup$ That is great contribution, george2079. But my problem explicitly orders to avoid the continuous case, meaning the discretization of the $F_i$ can't simply disappear. This wasn't explicitly written in the OP here, because I didn't find it important - I will add it as a note now. $\endgroup$
    – skrat
    Mar 14, 2016 at 21:03
  • $\begingroup$ PS: You probably meant to convert the radians to degrees so I assume ArcSin[asol[[4]]] Pi/180 (and others...) is a small typo. $\endgroup$
    – skrat
    Mar 14, 2016 at 21:09
  • 1
    $\begingroup$ added a discrete plot version. $\endgroup$
    – george2079
    Mar 15, 2016 at 18:37
  • $\begingroup$ I have no ideat what to do to debug this code but this is the error I get: **Flatten::normal: Nonatomic expression expected at position 1 in Flatten[[Lambda]]. ** $\endgroup$
    – skrat
    Mar 16, 2016 at 15:46
  • 1
    $\begingroup$ Ok, there is a gap where the curves meet -- Added a fix.. $\endgroup$
    – george2079
    Mar 16, 2016 at 20:06

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.