2
$\begingroup$

I have a system of equations: \begin{equation} \dot{x}=y-\alpha xz^2 \\ \dot{y}=x+ \alpha yz^2 \\ \dot{z}=-z(x^2+y^2) \end{equation} How can I create a 3D phase plane of this system that I can use for global analysis? I want to be able to see how this system behaves for all initial conditions. Here is what I tried:


FunSystem[\[Alpha]_, i_, l_, k_] = {x'[t] == y[t] - \[Alpha] x[t] z[t]^2,y'[t] == x[t] + \[Alpha] y[t] z[t]^2,z'[t] == -z[t] (x[t]^2 + y[t]^2), x[0] == i, y[0] == l, z[0] == k};
sol = NDSolve[FunSystem[\[Alpha], 0, 0, 0], {x, y, z}, {t, 0, 100}];
ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 100},ColorFunction -> (Hue[#4] &), BoxRatios -> 1,AxesLabel -> {x[t], y[t], z[t]}, ImageSize -> 300]

Only a little dot is showing up because I set the initial conditions to all be 0. Is there any way I can create a manipulate plot to work with this?

$\endgroup$
2
  • $\begingroup$ Here's my initial scheme (but I don't have time to code or test it): Create list of lines/curves (= {{x1,y1,z1}, {x2,y2,z2}...}) as follows: Start with user defined ICs and get the lists from NDSolve. Make a pool of initial ICs (say a 3D grid covering the plot region). Repeat: (1) Remove points from the grid that within a given radius of the points in the current lines; (2) Pick a point from the pool of ICs, integrate with NDSolve, and add the line to your list; repeat (1) and (2) until there are no more ICs. -- You tend to get a plate of spaghetti, though. $\endgroup$
    – Michael E2
    Commented Mar 21, 2020 at 20:29
  • $\begingroup$ Do you have an image of what you would like the result to be? $\endgroup$
    – Michael E2
    Commented Mar 21, 2020 at 22:47

2 Answers 2

3
$\begingroup$

You have not given enough info on what range of parameters is expected. But I would use ParametricPlot3D for plotting the solution and then combine this with Graphics3D to show how the solution (red dot) changes in time over the solution trajectory.

Something like this. Feel free to modify as needed.

enter image description here

Manipulate[
 Module[{p, x, y, z, t, init, xs, ys, zs, eq},
  eq = {x'[t] == y[t] - a *x[t]*z[t]^2, y'[t] == x[t] + a*y[t]*z[t]^2,
     z'[t] == -z[t] ( x[t]^2 + y[t]^2)}; 
  init = {x[0] == x0, y[0] == y0, z[0] == z0};
  {xs, ys, zs} = 
   NDSolveValue[{eq, init}, {x, y, z}, {t, 0, tMax}, 
    Method -> "StiffnessSwitching"]; 
  p = ParametricPlot3D[{xs[t], ys[t], zs[t]}, {t, 0, tMax}, 
    PlotRange -> All,
    PlotStyle -> Thin,
    AxesLabel -> {"x", "y", "z"}, BaseStyle -> 14, ImageSize -> 400];
  Grid[{{Row[{"time : ", NumberForm[t0, {4, 2}]}]},
    {
     Show[p,
      Graphics3D[
       {
        Red, PointSize[0.03], Point[{xs[t0], ys[t0], zs[t0]}]
        }
       ]
      , PlotRange -> {{0, 100}, {0, 100}, {0, 50}}
      ]
     }}
   ]
  ]
 ,
 {{t0, 0, "time"}, 0, tMax, .01, Appearance -> "Labeled"},
 {{a, 0.01, "alpha"}, -3, 3, .01, Appearance -> "Labeled"},
 {{x0, 0.01, "x(0)"}, 0, 30, .01, Appearance -> "Labeled"},
 {{y0, 0.01, "y(0)"}, 0, 30, .01, Appearance -> "Labeled"},
 {{z0, 40, "z(0)"}, 0, 50, .01, Appearance -> "Labeled"},
 {{tMax, 5}, None},
 TrackedSymbols :> {t0, a, x0, y0, z0}
 ]
$\endgroup$
3
  • $\begingroup$ Is it possible to use reference.wolfram.com/language/ref/VectorPlot3D.html and also show flows for different ICs? $\endgroup$
    – Moo
    Commented Mar 21, 2020 at 3:03
  • 2
    $\begingroup$ @Moo I have not used vectorPlot3D much myself. But please feel free to post that as an alternative solution if you like, it looks like it should work. $\endgroup$
    – Nasser
    Commented Mar 21, 2020 at 3:32
  • $\begingroup$ mathematica.stackexchange.com/questions/249167/… Can you tell me how to build a vector field for a given trajectory? $\endgroup$
    – ayr
    Commented Jun 4, 2021 at 9:29
3
$\begingroup$

What about

With[{\[Alpha] = -1.4}, 
VectorPlot3D[{y - \[Alpha] x  z ^2,x + \[Alpha] y  z^2, -z  (x ^2 + y ^2)}, {x, 0, 100}, {y, 0,100}, {z, 0, 100}]]

enter image description here

This vectorplot might be combined with ParametricPlot3D for special solutions you're interested in.

$\endgroup$

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.