4
$\begingroup$

How would it be possible to plot in Mathematica the Poincaré surface of section, say on the plane y=0 of the streamlines of the so-called Hopf fibration which has a tangent vector field with components

u_x = (2A/d^2)*(x*z - r*y);
u_y = (2A/d^2)*(r*x + y*z);
u_z = (A/d^2)*(r^2 - x^2 - y^2 + z^2);

where

d=r^2+x^2+y^2+z^2;

and $A=const.$ and $r=const.$

The differential equations whose solutions are the trajectories of the Hopf-fibration are 

x'[t] = u_x/Abs[u];
y'[t] = u_y/Abs[u];
z'[t] = u_z/Abs[u];
Improve this question
$\endgroup$
6
  • $\begingroup$ What is u on the right hand sides of the differential equations? $\endgroup$
    – zhk
    Commented Dec 14, 2016 at 12:03
  • $\begingroup$ What are the initial conditions? $\endgroup$
    – zhk
    Commented Dec 14, 2016 at 12:37
  • $\begingroup$ @MMM u is the tangent vector with components u_x, u_y, u_z. The initial conditions are x(0)=x_0, y(0)=y_0, z(0)=z_0, with x_0, y_0, z_0 constants. $\endgroup$
    – SpaceChild
    Commented Dec 14, 2016 at 13:03
  • $\begingroup$ So, you meant to say that x'[t]=u_x/Abs[u_x]? Your system of differential equations is coupled nonlinear one, thus, we need specific numerical values for x0, y0 and z0 to be able to find a numerical solution. $\endgroup$
    – zhk
    Commented Dec 14, 2016 at 14:03
  • 2
    $\begingroup$ Have you seen this or this ? $\endgroup$
    – gpap
    Commented Dec 14, 2016 at 14:45

1 Answer 1

Reset to default
6
$\begingroup$

We can plot the components of u by using ContourPlot3D like this

A = 1; r = 1; d = r^2 + x^2 + y^2 + z^2;
ContourPlot3D[{(2 A/d^2)*(x*z - r*y), (2 A/d^2) (r*x + y*z), (A/
     d^2)*(r^2 - x^2 - y^2 + z^2)}, {x, -3, 3}, {y, -3, 3}, {z, -3, 
  3}]

enter image description here

For the above plot, I took random values for $A$ and $r$.

Since your system of differential equation is coupled and nonlinear, so, I will straight away go for NDSolve.

d = r^2 + x[t]^2 + y[t]^2 + z[t]^2;
ux = (2*A/d^2)*(x[t]*z[t] - r*y[t]);
uy = (2*A/d^2)*(r*x[t] + y[t]*z[t]);
uz = 1/2*(A/d^2)*(r^2 - x[t]^2 - y[t]^2 + z[t]^2);
u = Sqrt[ux^2 + uy^2 + uz^2];
soln = NDSolve[{x'[t] == ux/Abs[u], y'[t] == uy/Abs[u], 
    z'[t] == uz/Abs[u], x[0] == 0, y[0] == 1, z[0] == 0}, {x, y, 
    z}, {t, 0, 50}];

Finally ploting the results as a 3D, 

ParametricPlot3D[{x[t], y[t], z[t]} /. soln, {t, 0, 50}, 
 PlotRange -> All, MaxRecursion -> 8, AxesLabel -> {"x", "y", "z"}, ViewPoint -> Front]

enter image description here

Edit

In responce to your xz-plane view comment, @J.M. suggested ParametericPlot,

ParametricPlot[{x[t], z[t]} /. soln, {t, 0, 250}, PlotRange -> All, 
 MaxRecursion -> 8, AxesLabel -> {"x", "z"}]

enter image description here

For set of different initial conditions, you can do something like this,

sol[x0_?NumericQ, y0_?NumericQ, z0_?NumericQ] :=NDSolve[{x'[t] == ux/Abs[u], y'[t] == uy/Abs[u], z'[t] == uz/Abs[u],
         x[0] == x0, y[0] == y0, z[0] == z0}, {x, y, z}, {t, 0, 250}];
    ParametricPlot3D[
     Evaluate[{x[t], y[t], z[t]} /. sol[#, #, #] & /@ 
       Range[2, 5, .1]], {t, 0, 250}, PlotRange -> All, MaxRecursion -> 8,
      AxesLabel -> {"x", "y", "z"}]

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ if we multiply the component u_z by 1/2 we obtain a fibration. How can we take the map of this plot on the x-z plane? $\endgroup$
    – SpaceChild
    Commented Dec 14, 2016 at 16:21
  • $\begingroup$ Also for different initial values we take different curves. How can we depict those curves for various initial points in the same plot? $\endgroup$
    – SpaceChild
    Commented Dec 14, 2016 at 16:40
  • $\begingroup$ @DK13 To view just xz-plane use this ViewPoint -> Front and to label the axes use AxesLabel -> {"x", "y", "z"} $\endgroup$
    – zhk
    Commented Dec 14, 2016 at 16:43
  • $\begingroup$ ...or just use ParametricPlot[] on the x and z components. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Dec 14, 2016 at 16:48
  • $\begingroup$ My mistake, I want to take the piercings of the curve with tha y=0 plane (the Poincare map). I suppose that one has to use the routine WhenEven[]... $\endgroup$
    – SpaceChild
    Commented Dec 14, 2016 at 16:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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