Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Here I have a problem that probably needs two loops, but I am not sure how to implement them together. The code calculates M for various values of parameter G. Now I want to find M which crosses the real axis (Imaginary part is zero) for each α1 = (1.1, 2.2). So the code is working fine and for example for α1 = 1.1, I got the curve which crosses the real axis, so I have one real solution (for chosen segment G), now I want to find these real points M for various α1 and to plot M[α1].

 \[Alpha]1 = 1.1;
 \[Nu] = 0.001;
 Tt = 0;
  poly = -g^2 + 0.25 x^4 - (Tt + \[Alpha]1^2) x^2 + 2 \[Alpha]1 g x + 
 I \[Nu] g - I \[Nu] \[Alpha]1 x + 1;

 LM = {};
 Do[Lroots = x /. NSolve[(poly /. {g -> G}) == 0, x];
 res = 0;
 Do[root = Lroots[[i]];
 If[Im[root] >= 0, L1 = Drop[Lroots, {i}];
 res = res + 1/(Times @@ (root - L1));];, {i, 1, Length[Lroots]}];
 M = 1/((4*I*G^2)*res);
  AppendTo[LM, M];, {G, 0.01, 1.5, 0.001}];


   a1 = Min@Select[Table[Im[LM[[i]]], {i, 1, Length[LM]}], Positive];
   a2 = Max@Select[Table[Im[LM[[i]]], {i, 1, Length[LM]}], Negative];

   b1 = Table[LM[[i]], {i, 1, Length[LM]}] - a1*I;
   b2 = Table[LM[[i]], {i, 1, Length[LM]}] - a2*I;

   c1 = Select[b1, Im[#] == 0 &][[1]];
   c2 = Select[b2, Im[#] == 0 &][[1]];

   d1=Re[(c1 + c2)/2]

The curve crosses the real axis, but just for α1 = 1.1; now I need to extract for various α1 = 1.1, ..., 2.2, all points which cross the real axis.

 p1 = ListPlot[Table[{Re[LM[[i]]], Im[LM[[i]]]}, {i, 1, Length[LM]}], 
 PlotRange -> {{-1, 18}, {-10, 10}}]
share|improve this question
    
Yes Öskå, from the code I need just one point when the Imaginary part is zero and then to calculate it for various alpha –  Pipe May 6 at 11:58
    
not this point, just one which cross the real axis, if doesn't exist it is infinity –  Pipe May 6 at 12:02
    
Yes and maybe I should change the G and change the step to see the crossing real axis –  Pipe May 6 at 12:04
    
Öskå please take a look now. I made a mistake because of segment. Now for each alpha the plotted curve cross the real axis, I need to extract these points for each alpha and to plot the M(alpha) –  Pipe May 6 at 17:26
    
Öskå it is not important to be exact, this is because of chosen step, it is enough just to choose the point with smallest imaginary part and for this point to extract the real value for each alpha on the segment. For all values alpha which I changed the curve cross the real axis for sure, but I just need to do this automatically, not to calculate for each alpha manually –  Pipe May 6 at 22:36

1 Answer 1

up vote 1 down vote accepted

Tell me if the following fits you. Note that I haven't tried to optimize your code. The idea is simply to put everything under the same function getNearZero which depends on α1. Then you are free to do every operations you want and do your second loop over α1.

getNearZero[α1_] := 
  Module[{ν, Tt, Poly, g, G, x, M, LM, poly, res, root, Lroots, L1},
    ν = 0.001; Tt = 0;
    poly = -g^2 + 0.25 x^4 - (Tt + α1^2) x^2 + 2 α1 g x + I ν g - I ν α1 x + 1;
    LM = {};
    Do[Lroots = x /. NSolve[(poly /. {g -> G}) == 0, x];
      res = 0;
      Do[root = Lroots[[i]];
        If[Im[root] >= 0, L1 = Drop[Lroots, {i}];
        res = res + 1/(Times @@ (root - L1));];, {i, 1, Length[Lroots]}];
      M = 1/((4*I*G^2)*res);
      AppendTo[LM, M];, {G, 0.01, 1.5, 0.001}]; 
    (* return the Nearest to zero and LM *)
    {Extract[LM, Position[Im /@ LM, First@Nearest[Im /@ LM, 0]]], LM}]

Checking if that's the closest point to zero:

With[{LM = Last@getNearZero@1.3, nearZero = First@getNearZero@1.3}, 
  ListPlot[Table[{Re[LM[[i]]], Im[LM[[i]]]}, {i, 1, Length[LM]}], 
  PlotRange -> {{-1, 18}, {-10, 10}}, 
  Epilog -> {PointSize@0.02, Red, Point[{First@Re@nearZero, First@Im@nearZero}]}]]

enter image description here

Then you can run the getNearZero from 1.1 to 2.2 every di = 0.1 and ListLinePlot the result (note the Tooltip showing the corresponding α1).

data = {#, First@getNearZero@#} & /@ Range[1.1, 2.2, .1];
ListLinePlot[
  Tooltip[#, "α1=" <> ToString@#2] & @@@ 
  Thread[{Thread[{Flatten@Re[Last /@ data], Flatten@Im[Last /@ data]}], First /@ data}],
  Frame -> True, PlotMarkers -> Automatic]

enter image description here


If you want to plot the real part of the points nearest to zero in terms of α1:

ListLinePlot[Thread[{First /@ data, Flatten[Re /@ Last /@ data]}], 
  PlotMarkers -> Automatic, AxesLabel -> {"α1"}]

enter image description here

And finally, if you want to work with your newly edited ai, bi, ci and d1 the following is shorter:

d1 = First @ Mean[Re /@ Extract[LM, Position[Im /@ LM, #]] & /@ Nearest[Im /@ LM, 0, 2]]

where Nearest[Im /@ LM, 0, 2] replaces your ai, the bi are not needed anymore, and Re /@ Extract[LM, Position[Im /@ LM, #]] & replaces your ci.

share|improve this answer
    
Thank you very much Öskå. This is it. I should learn through your code –  Pipe May 7 at 10:20

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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