3
$\begingroup$

I am putting all of my code in one block for easy copy and paste.

sys = {u (1 - u + a v), r v (1 - v + b u)};
sys1 = sys /. {a -> 2, b -> 3, r -> 1};
equilibria = Solve[sys1 == {0, 0}];
equilibria = equilibria[[{2, 3, 4}]];
sol[{u0_, v0_}] := 
  NDSolveValue[{{u'[t], 
      v'[t]} == (sys1 /. {u -> u[t], v -> v[t]}), {u[0], v[0]} == {u0,
       v0},
    WhenEvent[Abs[u[t] - 1] > 1, "StopIntegration"],
    WhenEvent[Abs[v[t] - 1] > 1, "StopIntegration"]}, {u[t], 
    v[t]}, {t, -4, 4},
   "ExtrapolationHandler" -> {Indeterminate &, 
     "WarningMessage" -> False}];
initconds = {{0.1, 1.5}, {1.5, 0.1}, {0.1, 0.5}, {0.5, 0.1}};
toplot = Map[sol, initconds];
plt = ParametricPlot[toplot, {t, -4, 4}, PlotStyle -> Blue];
Show[vp, cp, plt, 
 Graphics[{Red, PointSize[Large], Point[{u, v}] /. equilibria,
   Text[Style["(0,0) nodal source", Black, Background -> White], {0.1,
      0.1}],
   Text[Style["(1,0) saddle", Black, Background -> White], {1, 0.1}],
   Text[Style["(0,1) saddle", Black, Background -> White], {0.1, 
     1.1}]}]]

The resulting figure follows.

enter image description here

I've only tested these concepts one time, on this one system. I am wondering what folks think. Will I get in trouble in general?

Grisha's Suggestion:

Clear["Global`*"];
sys = {u (1 - u + a v), r v (1 - v + b u)};
sys1 = sys /. {a -> 2, b -> 3, r -> 1};
vp = VectorPlot[sys1, {u, -0.1, 2}, {v, -0.1, 2},
   VectorScale -> {0.045, 0.9, None},
   VectorStyle -> {GrayLevel[0.8]},
   Axes -> True, AxesLabel -> {x, y}, AxesOrigin -> {0, 0}];
cp = ContourPlot[sys1, {u, -0.1, 2}, {v, -0.1, 2},
   ContourStyle -> {Orange, Green}];
equilibria = Solve[sys1 == {0, 0}];
equilibria = equilibria[[{2, 3, 4}]];
initconds = {{0.1, 1.5}, {1.5, 0.1}, {0.1, 0.5}, {0.5, 0.1}};
sp = StreamPlot[sys1, {u, -0.5, 2}, {v, -0.5, 2}, StreamScale -> None,
    StreamPoints -> initconds,
   StreamStyle -> {Blue, "Line"}];
Show[vp, cp, sp, 
 Graphics[{Red, PointSize[Large], Point[{u, v}] /. equilibria,
   Text[Style["(0,0) nodal source", Black, Background -> White], {0.1,
      0.1}],
   Text[Style["(1,0) saddle", Black, Background -> White], {1, 0.1}],
   Text[Style["(0,1) saddle", Black, Background -> White], {0.1, 
     1.1}]}]]

enter image description here

One can see that Grisha's suggestion worked nicely here.

Case where ab<1

Clear["Global`*"];
sys = {u (1 - u + a v), r v (1 - v + b u)};
sys2 = sys /. {a -> 1/3, b -> 1/4, r -> 1};
vp2 = VectorPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2},
   VectorScale -> {0.045, 0.9, None},
   VectorStyle -> {GrayLevel[0.8]},
   Axes -> True, AxesLabel -> {x, y}, AxesOrigin -> {0, 0}];
cp2 = ContourPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2},
   ContourStyle -> {Orange, Green}];
equilibria = Solve[sys2 == {0, 0}];
initconds2 = {{0.1, 1.5}, {1.5, 0.1}, {0.1, 0.5}, {0.5, 0.1}, {0.7, 
    2}, {2, 2}, {2, 0.7}};
sp2 = StreamPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2}, 
   StreamScale -> None, StreamPoints -> initconds2,
   StreamStyle -> {Blue, "Line"}];
Show[vp2, cp2, sp2, 
 Graphics[{Red, PointSize[Large], Point[{u, v}] /. equilibria,
   Text[Style["(0,0) nodal source", Black, Background -> White], {0.1,
      0.1}],
   Text[Style["(1,0) saddle", Black, Background -> White], {1, 0.1}],
   Text[Style["(0,1) saddle", Black, Background -> White], {0.1, 1.1}],
   Text[Style["(16/11,15/11) nodal sink", Black, 
     Background -> White], {1.5, 1.5}]}]]

Produces the following image.

enter image description here

Notice that it failed to use three initial conditions. Changing the initial condition {0.7, 2} to {0.7,1.9} fixes the problem.

Clear["Global`*"];
sys = {u (1 - u + a v), r v (1 - v + b u)};
sys2 = sys /. {a -> 1/3, b -> 1/4, r -> 1};
vp2 = VectorPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2},
   VectorScale -> {0.045, 0.9, None},
   VectorStyle -> {GrayLevel[0.8]},
   Axes -> True, AxesLabel -> {x, y}, AxesOrigin -> {0, 0}];
cp2 = ContourPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2},
   ContourStyle -> {Orange, Green}];
equilibria = Solve[sys2 == {0, 0}];
initconds2 = {{0.1, 1.5}, {1.5, 0.1}, {0.1, 0.5}, {0.5, 0.1}, {0.7, 
    1.9}, {2, 2}, {2, 0.7}};
sp2 = StreamPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2}, 
   StreamScale -> None, StreamPoints -> initconds2,
   StreamStyle -> {Blue, "Line"}];
Show[vp2, cp2, sp2, 
 Graphics[{Red, PointSize[Large], Point[{u, v}] /. equilibria,
   Text[Style["(0,0) nodal source", Black, Background -> White], {0.1,
      0.1}],
   Text[Style["(1,0) saddle", Black, Background -> White], {1, 0.1}],
   Text[Style["(0,1) saddle", Black, Background -> White], {0.1, 1.1}],
   Text[Style["(16/11,15/11) nodal sink", Black, 
     Background -> White], {1.5, 1.5}]}]]

Giving the following image.

enter image description here

Not bad, not bad at all. Curious why I had to make the initial condition point change I made.

However, note that now all the stream plots don't push very close to the nodal sink. Is there a way you can force them to move a bit closer?

My Technique Again My technique seemed to work again.

Clear["Global`*"];
sys = {u (1 - u + a v), r v (1 - v + b u)};
sys2 = sys /. {a -> 1/3, b -> 1/4, r -> 1};
vp2 = VectorPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2},
   VectorScale -> {0.045, 0.9, None},
   VectorStyle -> {GrayLevel[0.8]},
   Axes -> True, AxesLabel -> {x, y}, AxesOrigin -> {0, 0}];
cp2 = ContourPlot[sys2, {u, -0.1, 2}, {v, -0.1, 2},
   ContourStyle -> {Orange, Green}];
equilibria = Solve[sys2 == {0, 0}];
sol[{u0_, v0_}] := 
  NDSolveValue[{{u'[t], 
      v'[t]} == (sys2 /. {u -> u[t], v -> v[t]}), {u[0], v[0]} == {u0,
       v0},
    WhenEvent[Abs[u[t] - 1] > 1, "StopIntegration"],
    WhenEvent[Abs[v[t] - 1] > 1, "StopIntegration"]}, {u[t], 
    v[t]}, {t, -4, 4},
   "ExtrapolationHandler" -> {Indeterminate &, 
     "WarningMessage" -> False}];
initconds2 = {{0.1, 1.5}, {1.5, 0.1}, {0.1, 0.5}, {0.5, 0.1}, {0.7, 
    1.9}, {2, 2}, {2, 0.7}};
toplot2 = Map[sol, initconds2];
plt2 = ParametricPlot[toplot2, {t, -4, 4}, PlotStyle -> Blue];
Show[vp2, cp2, plt2, 
 Graphics[{Red, PointSize[Large], Point[{u, v}] /. equilibria,
   Text[Style["(0,0) nodal source", Black, Background -> White], {0.1,
      0.1}],
   Text[Style["(1,0) saddle", Black, Background -> White], {1, 0.1}],
   Text[Style["(0,1) saddle", Black, Background -> White], {0.1, 1.1}],
   Text[Style["(16/11,15/11) nodal sink", Black, 
     Background -> White], {1.5, 1.5}]}]]

Producing this image:

enter image description here

However, this time I got a couple warning messages:

enter image description here

I'd love to hear suggestions from folks.

$\endgroup$
9
  • $\begingroup$ @MichaelE2 What do you think? Is this a good idea for my students? I am amazed that I found that ExtrapolationHandler thing. $\endgroup$
    – David
    Apr 22, 2015 at 6:43
  • 1
    $\begingroup$ Why do you not use StreamPlot[sys1, {u, -0.5, 2}, {v, -0.5, 2}, StreamScale -> Full, StreamPoints -> initconds, StreamStyle -> {Blue, "Line"}] for the selected initial conditions? $\endgroup$ Apr 22, 2015 at 7:34
  • $\begingroup$ Why not plot in 4 quadrants to see arrows on all sides of saddle points? $\endgroup$
    – Narasimham
    Apr 22, 2015 at 10:29
  • $\begingroup$ @GrishaKirilin Used your suggestion and ran into a small problem. And I am looking for a few more suggestions. See the edits to my post. Thanks for your idea. It is really a good one that I might be able to easily use with my student. I will note however, that StreamStyle -> {Blue, "Line"} doesn't really work. Try StreamStyle -> {Blue, "ArrowArrow"}. Think we've seen these errors before. $\endgroup$
    – David
    Apr 22, 2015 at 20:47
  • $\begingroup$ (1) I often use Grisha's approach, but as you've discovered, it does many things automatically. I find it hard to control in some cases. It stops, deletes or moves streamlines that get too close to one another. (2) The messages occur because the initial points {2,2} and {2,0.7} start on the boundary of the WhenEvent conditions; in such a case, the WhenEvent is ignored and NDSolve tries to integrate to infinity. WhenEvent only detects crossing a boundary. If you change each 2 to 1.99, no messages are generated. (3) I don't quite understand the purpose of the ContourPlot. $\endgroup$
    – Michael E2
    Apr 22, 2015 at 21:48

0

Your Answer

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

Browse other questions tagged or ask your own question.