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 is the corresponding code

First the equations setup

Vcl[x_, y_, z_] := (-G*(Mcl/a))/Sqrt[(x^2 + y^2 + z^2)/a^2 + c^2];
V[x_, y_, z_] := Vcl[x, y, z] + 1/2*(κ2 - 4*ω^2)*x^2 + 1/2*v2*z^2;

G = 1;
Mcl = 2.2; a = 0.182; c = 1;
ω = 1; κ2 = 1.8; v2 = 7.6;
cl = -3.264444506;
ch = 0.1;
E0 = cl*(1 - ch);

Then the DO loop with the Solve

data = {};
Do[
  Vx = V[x, 0, z0];
  sol = Solve[Vx == E0 && -1.01 < x < 1.01,x];
  xmin = x /. sol[[1]];
  xmax = x /. sol[[2]];
  AppendTo[data, {xmin, xmax, z0}],
  {z0, 0, 0.6, 0.001}
  ]

And finally the plot

neg = data[[All, {1, 3}]];
pos = data[[All, {2, 3}]];
L1 = ListPlot[neg, Joined -> True, 
     PlotStyle -> {Black, Thickness[0.003]}];
L2 = ListPlot[pos, Joined -> True, 
     PlotStyle -> {Black, Thickness[0.003]}]; 
L0 = Show[{L1, L2}, Axes -> False, Frame -> True, 
     PlotRange -> {{-1, 1}, {0, 0.6}}, ImageSize -> 550]

which produces this enter image description here

The plot has several issues: (i) it should be a continuous line from -1 to 1, while we observe a gap near 0 and also it never reaches +1, (ii) the lower horizontal line it should not be there; probably Joined joins two extreme points. I suspect, that both issues are related to the Solve inside the DO loop. There z0 belongs to the interval [0, 0.6] but the equation does not have solutions for all these values. The data list should somehow store only the correct solutions and reject the cases where there are no solutions at all.

Many thanks in advance.

share|improve this question
    
There are some warnings when solving your equation in the Do loop. You'd better correct the code. By the way, the solutions to your equations contains complex numbers, so the condition -1.01 < x < 1.01 might give you the void set. –  Z-Y.L Apr 10 at 8:28
    
@Kuba I think that the flood of error messages is because for some values of z0 there are no solutions. This is exactly what I want to know; how to get reject these cases from the list. The Solve syntax is corrected. BTW, thanks for the shorter version of neg! –  Vaggelis_Z Apr 10 at 8:36
    
@Z-Y.L I think that the warnings when solving the equations in the Do loop are mainly because there are no solutions for every value of z0; –  Vaggelis_Z Apr 10 at 8:38
    
@Kuba So is there no solution to the problem at all? –  Vaggelis_Z Apr 10 at 8:54
3  
Why not use ContourPlot[V[x, 0, z0] == E0, {x, -2, 2}, {z0, 0, 0.6}] –  Simon Woods Apr 10 at 10:05
show 3 more comments

3 Answers 3

up vote 3 down vote accepted

You could avoid the problem altogether by using ContourPlot:

ContourPlot[V[x, 0, z0] == E0, {x, -2, 2}, {z0, 0, 0.6}]

enter image description here

share|improve this answer
add comment

Using a finer grid in the loop (e.g. 0.0001) helps to reduce the gap in the middle.

If we do not join the data points and give different colors to your neg and pos lists, we see that something strange happens for small $z_0$: Your equation has 4 solutions (two negative, two positive), so just using the first or the second argument wont pick the smallest or largest.

 neg = data[[All, {1, 3}]];
 pos = data[[All, {2, 3}]];
 L1 = ListPlot[neg, PlotStyle -> {Blue, Thickness[0.002]}];
 L2 = ListPlot[pos, PlotStyle -> {Red, Thickness[0.002]}];
 L0 = Show[{L1, L2}, Axes -> False, Frame -> True, 
 PlotRange -> {{-1, 1}, {0, 0.6}}, ImageSize -> 550]

enter image description here

We can also see this if we plot $V$ and $E_0$ together for different $z_0$ (just at the beginning you can see it)

 Export["~/test.gif", Table[Plot[{V[x, 0, z0], E0}, {x, -1.01, 1.01}, PlotRange -> {-4.5, -1.5}, Frame -> True], {z0, 0.2, 0.55, 0.01}], "GIF"]

enter image description here

If you want xmin and xmax you will always have a gap close to +/- 1, but you can do it like this:

 data = {};
 Do[Vx = V[x, 0, z0];
    sol = Solve[Vx == E0 && -1.01 < x < 1.01, x];
    xmin = x /. First@sol;
    xmax = x /. Last@sol;
    AppendTo[data, {xmin, xmax, z0}], {z0, 0.25, 0.55, 0.0001}]

enter image description here

If you want all real solutions between -1.01 and 1.01 you can do this:

data = {};
Do[Vx = V[x, 0, z0];
  sol = Solve[Vx == E0 && -1.01 < x < 1.01, x];
  allx = {x, z0} /. sol;
  If[NumberQ[First@Flatten@allx],
     AppendTo[data, Flatten[{allx}, 1]]], 
{z0, 0.25, 0.55, 0.0001}]
ListPlot[Flatten[data, 1]]

where it seems a little hacky with all the Flatten and I am checking if there is a solution before I append it to the list, but it works and you get (maybe) what you want (where using a finer grid for z0 reduces the gap even more):

enter image description here

share|improve this answer
    
This would be improved if you use Reap/Sow instead of AppendTo –  george2079 Apr 10 at 18:21
    
sure, I just didn't want to obfuscate the important changes, that help to solve the main problem. –  jenson Apr 11 at 5:22
add comment

Here is a 'cheat':

p1 = Plot3D[V[x, 0, z], {x, -1.01, 1.01}, {z, 0, 0.6}, 
  MeshFunctions -> (#3 &), Mesh -> {{E0}}, MeshStyle -> {Red, Thick}, 
  PlotPoints -> 100]

enter image description here

You can then extract the desired mesh points from the object:

lns = Cases[p1, Line[x_] :> x, Infinity];
gr = Graphics[{Red, Thick, 
   Line[(p1[[1, 1]][[lns[[2]]]])[[All, {1, 2}]]]}, Axes -> True, 
  Frame -> True, AspectRatio -> 1, AxesOrigin -> {0, 0}]

enter image description here

share|improve this answer
    
See above a much simpler solution by @Simon Woods! –  Vaggelis_Z Apr 10 at 10:25
    
@Vaggelis_Z thank you...should have thought of that...and should have read the comments... –  ubpdqn Apr 10 at 10:28
    
@ubpdqn at least you're keeping it funky :) –  Jacob Akkerboom Apr 10 at 16:18
    
Ingenuous yet horribly convoluted. I like it! +1 –  Pillsy Apr 10 at 17:12
add comment

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.