2
$\begingroup$

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?

$\endgroup$
5
$\begingroup$

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

$\endgroup$
4
$\begingroup$
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

$\endgroup$
  • $\begingroup$ Thanks a lot, but if i want another section respect the one i found with mesh how can i modify Your code? $\endgroup$ – Panichi Pattumeros PapaCastoro Mar 4 '14 at 20:13

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.