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$
  • 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 '16 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$ – Dr. belisarius Mar 14 '16 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 '16 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 '16 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 '16 at 16:42
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$
  • $\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 '16 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 '16 at 21:09
  • 1
    $\begingroup$ added a discrete plot version. $\endgroup$ – george2079 Mar 15 '16 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 '16 at 15:46
  • 1
    $\begingroup$ Ok, there is a gap where the curves meet -- Added a fix.. $\endgroup$ – george2079 Mar 16 '16 at 20:06

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.