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.

I try to solve a ODE system of equations and i want to plot the intersection of the 3D solutions with a costant plane.

To do this i use the follow code:

Clear["Global`*"]
tfin = 500;
T = p1^2/2 + p2^2/2;
V = x^2/2 + y^2/2 + x^2 y - y^3/3;
H = T + V;
eq1 = x'[t] == D[H, p1] /. p1 -> p1[t];
eq2 = y'[t] == D[H, p2] /. p2 -> p2[t];
eq3 = p1'[t] == -D[H, x] /. {x -> x[t], y -> y[t]};
eq4 = p2'[t] == -D[H, y] /. {x -> x[t], y -> y[t]};
energy = (H /. {x -> -.1, y -> -.2, p2 -> -.05}) == 0.06;
sol = Solve[energy, p1];
p10 = p1 /. sol[[2]];
sol2 = NDSolve[{eq1, eq2, eq3, eq4,
    x[0] == -.1,
    y[0] == -.2,
    p1[0] == p10,
    p2[0] == -.05},
   {x[t], y[t], p1[t], p2[t]}, {t, 0, tfin}, MaxSteps -> 5000];
a = ParametricPlot[Evaluate[{x[t], y[t]} /. sol2], {t, 0, 20}];
b = ParametricPlot[Evaluate[{x[t], p1[t]} /. sol2], {t, 0, 20}];
c = ParametricPlot[Evaluate[{y[t], p2[t]} /. sol2], {t, 0, 20}];


gr1 = ParametricPlot3D[Evaluate[{x[t], y[t], p2[t]} /. sol2],
   {t, 0, tfin}, 
   PlotPoints -> 4000, BoxRatios -> {1, 1, 1}, ViewPoint -> {1, 0, -2},
   DisplayFunction -> Identity];

This is the 'plane' of intersection. Maybe I must use something like PLot3D[0.,{y,ymin,ymax},{z,zmin,zmax}] because i want a plane that intersect the orbit in the plane x=0.

gr2 = Graphics3D[{Green, Opacity[0.9], 
    Polygon[{{0, -.31, -.31}, {0, -.35, .35}, {0, .35, .35}, {0, .35, \
-.31}}]}, DisplayFunction -> Identity];

This is the 'intersection' i only plot the solutions with a very thin x-axis but this is not a plot of the points of intersection between the solutions and the plane (pl1):

ParametricPlot3D[Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin},
 PlotPoints -> 4000, PlotStyle -> Directive[Red, Thick], 
 PlotRange -> {{0, 0.002}, {-.4, .4}, {-.4, .4}},
 ViewPoint -> {1, 0, 0}, AxesLabel -> {"", "y", "p2"}, 
 ImageSize -> {500, 500}]

Show[gr1, gr2, DisplayFunction -> $DisplayFunction, 
 ImageSize -> {500, 500}]

Can someone help me please?

share|improve this question

2 Answers 2

up vote 4 down vote accepted

Use Mesh and MeshFunctions:

gr1 = ParametricPlot3D[
  Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin}, 
  PlotPoints -> 4000, MeshFunctions -> {#1 &}, Mesh -> {{0}}, 
  MeshStyle -> {Directive[PointSize[Medium], Red]}, 
  BoxRatios -> {1, 1, 1}, ViewPoint -> {1, 0, -2}, 
  DisplayFunction -> Identity]

Mathematica graphics

share|improve this answer
getOneCluster[pts_List, maxDist_?NumericQ] :=(*Returns a cluster*)
 Module[{f},
  f = Nearest[pts];
  FixedPoint[Union@Flatten[f[#, {Infinity, maxDist}] & /@ #, 1] &, {First@pts}]]
clusters[data_] := 
 Module[{f, dist},(*Some Characteristic Distance, assuming no isolated points*)
  f = Nearest[data];
  dist = 3 Max[EuclideanDistance[Last@f[#, 2], #] & /@ data];
  Flatten[ Reap[NestWhile[Complement[#, 
           Sow@getOneCluster[#, dist]] &,  data, # != {} &]][[2]], 1]]

p1 = ParametricPlot3D[ Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin}, 
                       MeshFunctions -> {#1 &}, Mesh -> {{0}}, PlotStyle -> None, 
                       MeshStyle -> Red];
gc = Cases[p1, GraphicsComplex[__], Infinity];
data = First[(Normal@gc)[[1, 2, 1, 2 ;;]] /. Point -> Sequence][[All, 2 ;;]];
clusters[data];
Show[ListPlot@clusters@data, ListLinePlot[#[[FindCurvePath[#][[1]]]] & /@ clusters[data]]]

Mathematica graphics

share|improve this answer
    
Thanks a lot, but if i want another section respect the one i found with mesh how can i modify Your code? –  Panichi Pattumeros PapaCastoro Mar 4 at 20:13

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.