25
$\begingroup$

Consider the system:

\begin{align*} x'&=(1-x-y)x\\ y'&=(4-7x-3y)y \end{align*}

The system has a saddle point at (1/4,3/4). How can I plot the separatrices on the phase portrait having domain $\{(x,y):\ 0\le x\le 1,\ 0\le y\le 2\}$?

Here is my attempt. Start with a vector plot.

Clear[x, y];
Clear[Derivative];
f[x_, y_] = (1 - x - y) x;
g[x_, y_] = (4 - 7 x - 3 y) y;
vp = VectorPlot[{f[x, y], g[x, y]}, {x, 0, 1}, {y, 0, 2},
   VectorScale -> {0.045, 0.5, None},
   VectorStyle -> {GrayLevel[0.8]},
   VectorPoints -> 16,
   Axes -> True, AxesLabel -> {x, y}];

Get the equilibrium points.

eqpts = Solve[{f[x, y] == 0, g[x, y] == 0}, {x, y}]

There is a saddle point at (1/4,3/4). Set an eps value.

eps=1/10000;

Define a function.

sol[{x0_, y0_, tmin_, tmax_}] := 
 NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], 
   x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, tmin, tmax}]

Spend time (brute force) adjusting arguments to sol function to plot the separatrices.

sep1 = ParametricPlot[
   Evaluate@sol[{1/4 + eps, 3/4, 0, 40}], {t, 0, 40},
   PlotStyle -> {Red, Thick}];
sep2 = ParametricPlot[
   Evaluate@sol[{1/4 - eps, 3/4, 0, 40}], {t, 0, 40},
   PlotStyle -> {Red, Thick}];
sep3 = ParametricPlot[
   Evaluate@sol[{1/4 + eps, 3/4 + eps, -2.8, 0}], {t, -2.8, 0},
   PlotStyle -> {Red, Thick}];
sep4 = ParametricPlot[
   Evaluate@sol[{1/4 - eps, 3/4 - eps, -40, 0}], {t, -40, 0},
   PlotStyle -> {Red, Thick}];

Show everything together.

Show[vp, Graphics[{Black, PointSize[Large], 
   Point[{x, y}] /. eqpts}], sep1, sep2, sep3, sep4]

Result:

enter image description here

Which worked. Just wondering if there would be a simpler approach, one that beginning students can easily understand.

P.S. All code pasted below for an easy copy and paste.

Clear[x, y];
Clear[Derivative];
f[x_, y_] = (1 - x - y) x;
g[x_, y_] = (4 - 7 x - 3 y) y;
vp = VectorPlot[{f[x, y], g[x, y]}, {x, 0, 1}, {y, 0, 2},
   VectorScale -> {0.045, 0.5, None},
   VectorStyle -> {GrayLevel[0.8]},
   VectorPoints -> 16,
   Axes -> True, AxesLabel -> {x, y}];
eqpts = Solve[{f[x, y] == 0, g[x, y] == 0}, {x, y}];
eps = 1/10000;
sol[{x0_, y0_, tmin_, tmax_}] := 
  NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], 
    x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, tmin, tmax}];
sep1 = ParametricPlot[
   Evaluate@sol[{1/4 + eps, 3/4, 0, 40}], {t, 0, 40},
   PlotStyle -> {Red, Thick}];
sep2 = ParametricPlot[
   Evaluate@sol[{1/4 - eps, 3/4, 0, 40}], {t, 0, 40},
   PlotStyle -> {Red, Thick}];
sep3 = ParametricPlot[
   Evaluate@sol[{1/4 + eps, 3/4 + eps, -2.8, 0}], {t, -2.8, 0},
   PlotStyle -> {Red, Thick}];
sep4 = ParametricPlot[
   Evaluate@sol[{1/4 - eps, 3/4 - eps, -40, 0}], {t, -40, 0},
   PlotStyle -> {Red, Thick}];
Show[vp, Graphics[{Black, PointSize[Large], 
   Point[{x, y}] /. eqpts}], sep1, sep2, sep3, sep4]
Improve this question
$\endgroup$
1
  • $\begingroup$ Take a look at this post, might get you started. $\endgroup$
    – ciao
    Apr 19, 2015 at 5:59

2 Answers 2

Reset to default
23
$\begingroup$

We can solve (approximately) for the initial conditions of solutions that approach an equilbrium by comparing the displacement vector from the equilibrium with the vector field of the ODE. Such trajectory is characterized by the condition that these two vectors become parallel as the solution nears the equilibrium. I used a similar idea before, which is buried in this answer.

sys = {(1 - x - y) x, (4 - 7 x - 3 y) y};
vars = {x, y};
equilibria = Solve[sys == {0, 0}, vars, Reals]
(*  {{x -> 0, y -> 4/3}, {x -> 1/4, y -> 3/4}, {x -> 1, y -> 0}, {x -> 0, y -> 0}}  *)

saddles = Pick[equilibria, Sign@Det@D[sys, {vars}] /. equilibria, -1]
(*  {{x -> 1/4, y -> 3/4}}  *)

Here is a function to get initial conditions for the separatrices.

sepICS[p0_, eps_] := With[{p1 = p0 + eps * Norm[p0] {Cos[t], Sin[t]}},
   p1 /. NSolve[Det[{p1 - p0, sys /. Thread[vars -> p1]}] == 0 && 0 <= t < 2 Pi]
   ];

We get parametrizations for the separatrices, stopping the integration when the solution gets close to the saddle and when it leaves the plot domain. One problem with approaching a saddle point is that the initial condition, as well as the subsequent integration, is approximate. If the solution is pushed too far, it will ricochet off along another separatrix (approximately).

separatrices = Flatten[
   Module[{eps = 10^-7,         (* tunable distance from equilibrium *)
     X0, dX},
    With[{Xa = 0, Xb = 1, Yc = 0, Yd = 2,  (* plot domain boundaries *)
        p0 = vars /. #,         (* equilibrium *)
        X = Through[vars[t]]},  (* variables at t *)
       X0 = X /. t -> 0;        (* initial values *)
       dX = D[X, t];            (* derivatives *)
       With[{X1 = X[[1]], X2 = X[[2]]},
        First@NDSolve[{
             dX == (sys /. v : Alternatives @@ vars :> v[t]), 
             X0 == #,
             (*stop when close to saddle*)
             WhenEvent[Norm[X - p0] < 0.5 eps * Norm[p0], "StopIntegration"],
             (*stop when solution leaves plot domain*)
             WhenEvent[Abs[X1 - (Xa + Xb)/2] > (Xb - Xa)/2, "StopIntegration"], 
             WhenEvent[Abs[X2 - (Yc + Yd)/2] > (Yd - Yc)/2, "StopIntegration"]},
            vars, {t, -100, 100}] & /@ sepICS[p0, eps]
        ]] & /@ saddles
    ],
 1];


sepPlots = ParametricPlot @@@ ({{x[t], y[t]}, 
       Hold[Flatten][{t, x["Domain"]}]} /. separatrices // ReleaseHold);

Show[
 background,
 vp,
 sepPlots,
 PlotRange -> All, Frame -> True, Axes -> False, AspectRatio -> 1]

Mathematica graphics

The background is some attempt at mimicking ubpqdn's. One can use the lines from the plots of the separatrices to construct poiygon's to illustrate the regions created by them. It is a bit awkward to add the corner points. 

sepLines = First@Cases[#, _Line, Infinity] & /@ sepPlots;
background = Graphics[
  Riffle[
   Lighter[#, 0.6] & /@ {Red, Purple, Yellow, Green},
   MapAt[Reverse, 
     Partition[sepLines, 2, 1, 1], {{2}, {4}}] /. {Line[p1_], 
      Line[p2_]} :> Polygon[
      Join[p1, p2,
       Nearest[Tuples[{{0, 1}, {0, 2}}], {p2[[-1, 1]], p1[[1, 2]]}, {1, 0.01}],
       Nearest[Tuples[{{0, 1}, {0, 2}}], {p1[[1, 1]], p2[[-1, 2]]}, {1, 0.01}]]]
   ],
  PlotRange -> All, Frame -> True, Axes -> False, AspectRatio -> 1
  ]

Addendum: Notes on the code.

1. v : x | y :> v[t]: The : is short for Pattern; | is short for Alternatives. So v : x | y defines the pattern symbol v to represent x or y. The whole v : x | y :> v[t], means replace x or y by x[t] or y[t] respectively. The rules {x -> x[t], y -> y[t]} are equivalent.

2. ParametricPlot @@@ ...: The main problem here are the domains of the solutions. The variable separatrices contains a list of solutions of the form

{{x -> x1ifn, y -> y1ifn}, {x -> x2ifn, y -> y2ifn},...}

where x1ifn, y1ifn etc. are interpolating functions. Each pair x1ifn, y1ifn has the same domain, but another pair x2ifn, y2ifn will have a different domain. So what is an easy way to plot all of the solutions? If ifn is an InterpolatingFunction, then ifn["Domain"] returns a list of domains for each input; in this case, it will have the form {{tmin, tmax}}. Flatten[{t, x["Domain"]}] will have the form {t, tmin, tmax} as needed for ParametricPlot. The problem is that the x in x["Domain"] has to be replaced by an InterpolatingFunction and evaluated before Flatten is evaluated. Hence the Hold[Flatten], to prevent flattening until after the /. separatrices has been executed; the ReleaseHold then lets the domains be evaluated and flattened. Since separatrices is a list of solutions (each of which is a list of Rules), the replacement yields a list of the form:

{{{x1ifn[t], y1ifn[t]}, Hold[Flatten][{t, {{tmin1, tmax1}}}]},
 {{x2ifn[t], y2ifn[t]}, Hold[Flatten][{t, {{tmin2, tmax2}}}]},
 ...}

After ReleaseHold, these elements will be ready to have ParametricPlot applied to them with @@@. This replaces the {} around each element with ParametricPlot[]:

{ParametricPlot[{x1ifn[t], y1ifn[t]}, {t, tmin1, tmax1}],
 ParametricPlot[{x2ifn[t], y2ifn[t]}, {t, tmin2, tmax2}],
 ...}

These automatically evaluate to the plots of each separatrix.

Improve this answer
$\endgroup$
7
  • $\begingroup$ this is beatiful and correct answer...I will direc OP towards your answer...I was merely trying to point but am glad have learned a lot from your answer in process $\endgroup$
    – ubpdqn
    Apr 20, 2015 at 2:18
  • $\begingroup$ @ubpdqn Thanks! I also realized I forgot to say how the displacement and the vector field are to be compared. $\endgroup$
    – Michael E2
    Apr 20, 2015 at 2:38
  • $\begingroup$ @MichaelE2 Absolutely wonderful answer. But I am having difficulty figuring it out. My best first question is what does (sys /. v : x | y :> v[t]) mean? I know that sys is {x (1 - x - y), (4 - 7 x - 3 y) y}. Then (sys /. v : x | y :> v[t]) returns {x[t] (1 - x[t] - y[t]), (4 - 7 x[t] - 3 y[t]) y[t]}. So it seems to replace x with x[t] and y with y[t], but why? What does v:x | y:>v[t] mean? $\endgroup$
    – David
    Apr 20, 2015 at 4:11
  • $\begingroup$ @MichaelE2 Excellent answer, easily understood by me. May I now ask my next question? This command is puzzling: sepPlots = ParametricPlot @@@ ({{x[t], y[t]}, Hold[Flatten][{t, x["Domain"]}]} /. separatrices // ReleaseHold); Could you add a nice explanation of this above? $\endgroup$
    – David
    Apr 20, 2015 at 18:27
  • 1
    $\begingroup$ @Axel You could ask your question as a question. Refer to this Q&A and include the code for your system. $\endgroup$
    – Michael E2
    Jun 16, 2019 at 16:31
9
$\begingroup$

Perhaps this would be helpful:

e1 = (1 - x - y) x;
e2 = (4 - 7 x - 3 y) y;
sp = StreamPlot[{e1, e2}, {x, 0, 1}, {y, 0, 2}, 
   Epilog -> {Red, PointSize[0.02], Point[{{0, 0}, {1/4, 3/4}}]}];
{a, b, c, d} = 
  ImplicitRegion[#, {x, 
      y}] & /@ ({#1[e1, 0] && #2[e2, 0] && 0 < x < 1 && 0 < y < 2} & @@@
      Tuples[{Less, Greater}, 2]);
lab = "x'(t)" <> #1 <> ", " <> "y'(t)" <> #2 & @@@ 
   Tuples[{"<0", ">0"}, 2];
With[{col = {Red, Green, Blue, Yellow}}, 
 Legended[Show[sp, 
   Sequence @@ 
    MapThread[
     RegionPlot[DiscretizeRegion[#, PrecisionGoal -> 10], 
       PlotStyle -> {Opacity[0.3], #2}] &, {{a, b, c, d}, col}]], 
  SwatchLegend[Lighter[#, 0.7] & /@ col, lab, 
   LegendFunction -> "Panel"]]]

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ Beautiful! I definitely can use this in class. However, what I need is the same image, but with the separatrices, the solutions moving toward and away from the equilibrium point. $\endgroup$
    – David
    Apr 19, 2015 at 15:02
  • $\begingroup$ @David I would look at MichaelE2 answer and accept it as my answer merely partitions region of interest and not what you want and there is a lot of instructive code :) $\endgroup$
    – ubpdqn
    Apr 20, 2015 at 2:21
  • $\begingroup$ You were both equally helpful. I really appreciate your answers. I'm going to have to give a little effort to understand the coding. $\endgroup$
    – David
    Apr 20, 2015 at 6:22
  • 1
    $\begingroup$ @David MichaelE2's answer is beautiful, instructive and correct...mine was half-baked (if well intentioned)...done without deep thought but aimed to motivate diving in...glad the MichaelE2's answer emerged and I have also learned a lot...this is joy of this site... $\endgroup$
    – ubpdqn
    Apr 20, 2015 at 6:55

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.