8
$\begingroup$

Suppose we have the following simplified system of two ordinary differential equations:

$$\dot{x}(t)=x(t)^2+2y(t)\\ \dot{y}(t)=3x(t)$$

The system has a hyperbolic fixed point the origin. Hence there exits a stable and an unstable manifold going through the origin.

I would like to plot these exact sets. I am able to plot basins of attraction for regions in phase space but I fail to do so once the attracting region is a one-dimensional line (as is the case here with the stable manifold).

streamplot

Edit: With MMM's help I was able to roughly plot the unstable manifold. However, what is more challenging is to plot the stable one.

Eq1 = x'[t] == x[t]^2 + 2*y[t];
Eq2 = y'[t] == 3*x[t];
splot = StreamPlot[{x^2 + 2*y, 3*x}, {x, -6, 6}, {y, -6, 6}];
Show[splot, 
 ParametricPlot[
  Evaluate[First[{x[t], y[t]} /. 
  NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {-0.01, 0}]}, {x, 
   y}, {t, 0, 3}]]], {t, 0, 3}, PlotStyle -> Red], 
 ParametricPlot[
  Evaluate[First[{x[t], y[t]} /. 
     NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {0.01, 0}]}, {x, 
       y}, {t, 0, 2.4}]]], {t, 0, 2.4}, PlotStyle -> Red]]

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Share your try in Mathematica. $\endgroup$ – zhk Jan 13 '17 at 10:15
5
$\begingroup$

The trick to getting the unstable manifold is to run backwards in time. Here I added two more solutions starting near the equilibrium but with negative t:

Eq1 = x'[t] == x[t]^2 + 2*y[t];
Eq2 = y'[t] == 3*x[t];
splot = StreamPlot[{x^2 + 2*y, 3*x}, {x, -6, 6}, {y, -6, 6}];
Show[splot, 
  ParametricPlot[Evaluate[First[{x[t], y[t]} /. 
  NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {-0.01, 0}]}, {x, y}, {t, 0, 4}]]],
  {t, 0, 4}, PlotStyle -> Red, PlotRange -> {{-6, 6}, {-6, 6}}], 
  ParametricPlot[Evaluate[First[{x[t], y[t]} /. 
  NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {0.01, 0}]}, {x, y}, {t, 0, 2.6}]]],
  {t, 0, 2.6}, PlotStyle -> Red, PlotRange -> {{-6, 6}, {-6, 6}}],
  ParametricPlot[Evaluate[First[{x[t], y[t]} /. 
  NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {-0.01, 0}]}, {x, y}, {t, 0, -2.6}]]],
  {t, 0, -2.6}, PlotStyle -> Green, PlotRange -> {{-6, 6}, {-6, 6}}], 
  ParametricPlot[Evaluate[First[{x[t], y[t]} /.
  NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == {0.01, 0}]}, {x, y}, {t, 0, -4}]]],
  {t, 0, -4}, PlotStyle -> Green, PlotRange -> {{-6, 6}, {-6, 6}}]
]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thank you, this is exactly what I wanted. $\endgroup$ – ParetoWilli Jan 13 '17 at 13:26
6
$\begingroup$

Here's one way to get StreamPlot to show the desired manifolds and fill in the rest of the plot around them.

We can use the eigenvectors of the Jacobian at the equilibrium to approximate a point on each manifold. From such a point, another point further away from the equilibrium along the manifold may be constructed with NDSolve. Using these for StreamPoints allows StreamPlot to fill in the rest of the phase portrait.

(* some useful quantities/expression for playing *)    
sys = {x'[t] == x[t]^2 + 2*y[t], y'[t] == 3*x[t]};
vars = {x, y};
dvars = D[Through[vars[t]], t];                   (* {x'[t], y'[t]} *)
vel = dvars /. First@Solve[sys, dvars];           (* expressions for {x'[t], y'[t]} *)
equil = First@Solve[sys /. Thread[dvars -> 0]];   (* equilibrium *)
phasefield = vel /. var_[t] :> var;               (* strip [t] from expressions *)
linearized = D[vel, {Through[vars[t]]}] /. equil; (* its Eigensystem[] describes equil. *)

(* compute stream points for the separatrix manifolds *)
sp = Function[{λ,      (* eigenvalue *)
               v,      (* eigenvector *)
               side},  (* ±1 = which side of the equilibrium: ± v *)
   Module[{tfinal},
    {NDSolveValue[{sys, Through[vars[0]] == 10^-4 v*side, 
       WhenEvent[Norm[#] == 1, tfinal = t; "StopIntegration"] &@
        Through[vars[t]]},
      Through[vars[tfinal]], {t, 0, 100 Sign@λ}],
     If[λ < 0, Green, Red]}
    ]
   ] @@@ Append @@@ Tuples[{Thread@Eigensystem@linearized, {-1, 1}}]
(*
  {{{ 0.599966, -0.800025}, RGBColor[0, 1, 0]},
   {{-0.664941,  0.746896}, RGBColor[0, 1, 0]},
   {{-0.599966, -0.800025}, RGBColor[1, 0, 0]},
   {{ 0.664941,  0.746896}, RGBColor[1, 0, 0]}}
*)

splot = StreamPlot[phasefield, {x, -6, 6}, {y, -6, 6},
  StreamPoints -> {Append[sp, Automatic]}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Nice! I think there may be a typo or two in the code, though...it's not working for me. $\endgroup$ – Simon Rochester Jan 15 '17 at 6:01
  • $\begingroup$ @SimonRochester Thanks for letting me know. I had tried to simplify a couple of things and forgot to change all references to old variables. $\endgroup$ – Michael E2 Jan 15 '17 at 13:29
2
$\begingroup$

Your question is not clear. But here is a starting for you.

Just adopting Chris answer to same sort of question,

Eq1 = x'[t] == x[t]^2 + 2*y[t];
Eq2 = y'[t] == 3*x[t];
splot = StreamPlot[{x^2 + 2*y, 3*x}, {x, -6, 6}, {y, -6, 6}, 
  StreamColorFunction -> "Rainbow"]
Manipulate[
 Show[splot, 
  ParametricPlot[
   Evaluate[
    First[{x[t], y[t]} /. 
      NDSolve[{Eq1, Eq2, Thread[{x[0], y[0]} == point]}, {x, y}, {t, 
        0, T}]]], {t, 0, T}, PlotStyle -> Red]], {{T, 1}, 0.1, 
  1}, {{point, {0.01, 0.0}}, Locator}, SaveDefinitions -> True]

enter image description here

$\endgroup$
  • $\begingroup$ more specifically: I want to plot the two (unique) trajectories that go through the saddle at the origin. One of them is the stable manifold, the other one is unstable. By studying the vector field we can get a good idea of how they look like, but i want to plot them. $\endgroup$ – ParetoWilli Jan 13 '17 at 11:27
  • 1
    $\begingroup$ @ParetoWilli I think you want something exactly like this. mathematica.stackexchange.com/questions/80284/… $\endgroup$ – zhk Jan 13 '17 at 12:43

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.