3
$\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?

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
6
$\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

Improve this answer
$\endgroup$
5
$\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

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks a lot, but if i want another section respect the one i found with mesh how can i modify Your code? $\endgroup$ Commented Mar 4, 2014 at 20:13

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.