5
$\begingroup$

I have an autonoous differential equation dx/dt = x(1-x)

dxdt[x_] := x*(1 - x)

I plotted the slope field with a few solutions (left) and the constant solutions (right)

SolsAll = VectorPlot[{1, dxdt[x]}, {t, 0, 2}, {x, -2, 2}, 
          StreamPoints -> Coarse, StreamColorFunction -> None, 
          StreamStyle -> Magenta, FrameLabel -> {"t", "x"}];

SolsPar = VectorPlot[{1, dxdt[x]}, {t, 0, 2}, {x, -2, 2}, 
          StreamPoints -> {{{{1, 0}, Magenta}, {{1, 1}, Magenta}}}, 
          StreamColorFunction -> None, FrameLabel -> {"t", "x"}];

GraphicsGrid[{{SolsAll, SolsPar}}]

enter image description here

I would like to add a phase line as in this picture. How can I do it?

enter image description here

Improve this question
$\endgroup$
3

1 Answer 1

Reset to default
8
$\begingroup$

Some years ago I started on something like this; then some years later, I made a package for my differential equations course, in which three plots were aligned to show three different views of equilibria in 1D (autonomous) systems (just as one finds in V. Arnold's Ordinary Differential Equations); some years after that, I started to extend and refactor that package but have not finished. All that is meant to prepare my apology for cutting code out of the package that is complicated and incompletely reworked. It seems to work just fine, which is one reason refactoring it is a low priority.

portrait1D[x'[t] == (1 - x[t]) (x[t] - 3), {x, 0, 4}, {t, 0, 3}]

Mathematica graphics

The code that draws the phase line portrait begins with a (*2*) in the code dump below and has another comment (*phase line field*), which may facilitate searching the long code. Here is the code with appropriate initializations to reproduce the middle plot above:

Block[{
  vf = -((-3 + x) (-1 + x)),
  fn = x, x1 = 0, x2 = 4,
  intervals1D = {{0.`, 1.`}, {1.`, 3.`}, {3.`, 4.`}},
  stationaryPts1D = {1, 3},
  init1DQ = True,
  plotrange, padding},
 plotrange = {All, {x1, x2}};
 padding = {0.1, (x2 - x1) {0.03, 0.02}};
 (*2*)
 Show[(*phase line field*)
  vectorPlot1D[vf, {fn, x1, x2}, "Vertical", 
   ImageSize -> {Automatic, 150}], 
  Graphics[{White, Point[{-0.7, x1}]}] (*extends left margin*), 
  PlotRange -> plotrange, PlotRangePadding -> padding, 
  Axes -> {False, True}, AxesLabel -> {None, fn /. x_[0] :> x}, 
  Ticks -> None]
 ]

There's a polar phase portrait in the package; I could have cut it out, but since it is built on some common utilities, I thought it generous to include it and hope it will be well-received.

polarPortrait1D[
 r'[t] == 0.25 r[t] (1 - r[t]) (r[t] - 3),
 {r, 0, 4}, {t, 0, 2 Pi}]

Mathematica graphics

(*
 * Phase portrait 
 *)

SetAttributes[odeVFComponents, Listable];
odeVFComponents[de : (Derivative[n_][f_][var_] -> nthDerivative_)] := 
  Table[Derivative[i][f][var], {i, n}] /. de;

derivativesToVariables[expr_, fns_, ___] := 
 expr /. {Derivative[n0_][f0 : Alternatives @@ Flatten[{fns}]][_] :> 
    f0[n0], (f1 : Alternatives @@ Flatten[{fns}])[_] :> f1[0]}

doOdeToVF[sols_List] := Flatten@odeVFComponents[sols]

odeToVF[eqn_, fns_, var_] := Module[{solution},
   doOdeToVF[solution] /; (
     (Head[eqn] === Equal || 
        VectorQ[eqn, 
         MatchQ[#, _Equal] &] || (Message[Solve::naqs, eqn]; False)) &&
            Check[
       solution = Solve[eqn, Last /@ Sort /@ GatherBy[
            Cases[eqn, 
             Derivative[_][Alternatives @@ Flatten[{fns}]][var], 
             Infinity]
            , # /. Derivative[_] -> Identity &]
         ]; True, False, {Solve::ivar}]
     )
   ];

(* need to convert ODE, List VFs to scalar VFs
*   to map var x[t] to x if necessary
*   to check vf is 1d func of x only
*   to initialize singularities and intervals
*   to add stationary point to de object
 *)
ClearAll[init1D, vectorPlot1D, portrait1D, polarPortrait1D];
SetAttributes[init1D, HoldAll];
(*
 * This is an important hook to call init1D on any 1D system plotter
*   init1D[] parses the ode/vector field and sets up stationary \
solutions
*     and the intervals between them
 *)
$plotters1D = vectorPlot1D | portrait1D | polarPortrait1D;
(call : #[args___] /; ! TrueQ[init1DQ] := 
     init1D[call]) & /@ $plotters1D;
(* * * * *)
zeros[vf_, {var_, a_, b_}] := 
  var /. Solve[vf == 0 && a <= var <= b, var, Reals] /. var -> {};
parseODE1D[ode_Equal, {x_[t_] | x_, x1_, x2_}, T : {t_, t1_, t2_}, 
  rest___] :=
 {odeToVF[ode, {x}, t] /. x[t] -> x, {x, x1, x2}, T, 
  rest}
parseODE1D[ode_Equal, {x_[t_], x1_, x2_}, rest___] /; ! 
    NumericQ[t] :=
  {odeToVF[ode, {x}, t] /. x[t] -> x, {x, x1, x2}, 
   rest};
parseODE1D[de : diffEq[asc_], rest___] := With[{
    x = getKey[de, "dependentVars", {}] /. {t1_} :> t1,
    t = getKey[de, "independentVars", {}] /. {t1_} :> t1},
   {odeToVF[getKey[de, "de", {}], {x}, t] /. x[t] -> x, rest} (** 
   TBD: Check rest **)
   ];
parseODE1D[___] := TBD["input error (diagnosis TBD)"]; (** TBD ??? **)
init1D[h_[ode : _Equal | _diffEq, args___]] := 
 init1D[h[##]] & @@ parseODE1D[ode, args];
init1D[h_[vf_, X : {x_, a_, b_}, args___]] := 
 Block[{init1DQ = True, stationaryPts1D, intervals1D}
     , stationaryPts1D = zeros[#, X]
     ; intervals1D = 
      Partition[Union[N[Join[{a, b}, stationaryPts1D]]], 2, 1]
     ; h[#, X, args]
     ] & /@ Flatten@vf /. {res_} :> res; (** Flatten@vf ??? **)
vectorPlot1D // Options = Join[{
    VectorPoints -> Automatic,
    PlotStyle -> Automatic,
    "VectorColors" -> {ColorData[1][1], ColorData[1][2]}},
    Options[Graphics]];
(* new diffEqClass` interfaces *)
vectorPlot1D[vf_, X : {_, _, _}, T : {_, _, _}, rest___] :=
  vectorPlot1D[vf, X, rest];
vectorPlot1D[{{vf_}}, {var_, a_, b_},  (* Needed ??? *)
    Optional[orientation : "Horizontal" | "Vertical", "Vertical"], 
    opts : OptionsPattern[]] /; True :=
  vectorPlot1D[vf, {var, a, b}, orientation, opts];
vectorPlot1D[vfs : {Repeated[{_}, {2, Infinity}]}, 
   Optional[orientation : "Horizontal" | "Vertical", 
    "Vertical"], {var_, a_, b_}, opts : OptionsPattern[]] :=
  vectorPlot1D[#, {var, a, b}, orientation, opts] & /@ Flatten@vfs;
(*vectorPlot1D[bad_,{var_,a_,b_},opts:OptionsPattern[]]:=(
Message[diffEq::CantHappen,"Plotting",Row[{"the 1D vector field is \
",bad}']];
Throw[$Failed]);*)
(* old math212` function *)
vectorPlot1D[vf_, {var_, a_, b_}, 
   Optional[orientation : "Horizontal" | "Vertical", "Vertical"], 
   opts : OptionsPattern[]] :=
  Module[{style, colorFn, coord}
   , If[orientation === "Horizontal"
    , coord[y_][x_] := {x, y}
    , coord[x_][y_] := {x, y}]
   ; style = OptionValue[PlotStyle] /. 
     Automatic -> Directive[Thickness[Medium], ColorData[1][1]]
   ; colorFn[{y1_, y2_}] :=
    If[(vf /. var -> Mean[{y1, y2}]) > 0
     , Sequence @@ OptionValue["VectorColors"]]
   ; Graphics[{
     Thickness[Medium], ColorData[1][1], style
     , {colorFn[#1], Line[Thread[coord[0][#]]]} & /@ intervals1D
     , {colorFn[#1], Arrowheads[Medium]
        , Arrow[If[(vf /. var -> Mean[#1]) > 0
           , Identity
           , Reverse
           ][Thread[coord[0.08 (b - a)][With[{scale = 0.9},
             {{scale, 1 - scale}, {1 - scale, scale}} . #1]]]
          ]
         ]
        } & /@ SortBy[N[intervals1D], Differences]
     , EdgeForm[style]
     , {Black, Text[Round[N[#1], 0.01], coord[-0.1][#], coord[1][0]]
        , White, Disk[coord[0][#], Offset[{2, 2}]]
        } & /@ stationaryPts1D
     }
    , FilterRules[{opts}, Options[Graphics]]
    ]
   ];

(* is fn/.x_[0]\[RuleDelayed]x still necessary? useful? *)
portrait1D[vf_, {fn_, x1_, x2_}, {time_, t1_, t2_}, opts___?OptionQ] :=
    Module[{plotrange = {All, {x1, x2}}, 
    padding = {0.1, (x2 - x1) {0.03, 0.02}}}
   , Grid[{
     {(*1*)
       ParametricPlot[ (* velocity plot *)
       {vf, fn}, {fn, x1, x2}
       , Mesh -> {{0}}, MeshFunctions -> {#1 &}, Ticks -> None, 
       ImageSize -> {Automatic, 150}, PlotRange -> plotrange, 
       PlotRangePadding -> padding , AspectRatio -> 2, 
       AxesLabel -> {HoldForm[Global`v], fn /. x_[0] :> x},
       Epilog -> {White, EdgeForm[Black], 
         Disk[{0, #}, Offset[{2, 2}]] & /@ stationaryPts1D}
       ]
      ,(*2*)
       Show[ (* phase line field *)
        vectorPlot1D[vf, {fn, x1, x2}, ImageSize -> {Automatic, 150}]
       , Graphics[{White, Point[{-0.7, x1}]}] (* 
       extends left margin *)
       , PlotRange -> plotrange, 
       PlotRangePadding -> padding, Axes -> {False, True}, 
       AxesLabel -> {None, fn /. x_[0] :> x}, Ticks -> None
       ]
      ,(*3*)
       Module[{eqns, x0} (* typical solutions *)
       , 
       eqns = {fn'[time] == (vf /. fn -> fn[time]), fn[0] == x0}
       ; With[{psol = ParametricNDSolveValue[eqns
           , fn, {time, t1, t2}, {x0}
           , 
           "ExtrapolationHandler" -> {Indeterminate &, 
             "WarningMessage" -> False}
           ]}
        , Plot[Evaluate[Quiet[
           (psol[#1][time] &) /@
            Join[stationaryPts1D
             , 
             If[(vf /. fn -> Mean[#1]) > 
                 0, {0.9, 0.1} . #1, {0.1, 0.9} . #1] & /@ 
              intervals1D
             ]]]
                          , {time, t1, t2}
                          , PlotRange -> plotrange, 
         PlotRangePadding -> padding , ImageSize -> {Automatic, 150}, 
         Ticks -> None, ExclusionsStyle -> None, 
         AxesLabel -> {time, fn /. x_[0] :> x}, 
         Epilog -> {White, EdgeForm[Black], 
           Disk[{0, #}, Offset[{2, 2}]] & /@ stationaryPts1D}
                          ]
                        ]
                  ]
                },
     (Row[#1, BaseStyle -> "Label"] &) /@ {{fn /. x_[0] :> x, " vs. ",
         HoldForm[Global`v]}, {"vector field"}, {fn /. x_[0] :> x, 
        " vs. " , time}}
            }]
   ];

(*
 * Polar direction field portrait
 * alternative APIs TBD ???
 *)
polarPortrait1D[vf_, {fn_, r1_, r2_}, {time_, t1_, t2_}, 
   opts : OptionsPattern@PolarPlot] := 
  Block[{plotrange = r2, padding = {0.1, (r2 - r1) {0.03, 0.02}}}
   , Module[{eqns, x0}
    , eqns = {
      fn'[time] == (vf /. fn -> fn[time]),
      fn[t1] == x0, (* IC *)
      Table[
       With[{r0 = r0}, 
        WhenEvent[Abs[fn[time] - r0] < 1*^-3, "StopIntegration"]],
       {r0, stationaryPts1D}]}
    ; With[
     {psol = ParametricNDSolveValue[eqns, fn, {time, t1, t2}, {x0},
        "ExtrapolationHandler" -> {Indeterminate &, 
          "WarningMessage" -> False}]}
     , PolarPlot[
       Evaluate@Flatten[List[
           ConditionalExpression[#, t1 <= time <= t1 + 2 Pi] & /@ 
           DeleteDuplicates@stationaryPts1D,
           Quiet[(psol[#1]@(time) &) /@ (If[(vf /. fn -> Mean[#1]) > 
                 0, {0.9, 0.1} . #1, {0.1, 0.9} . #1] &) /@ 
             intervals1D]]],
       {time, t1, t2},
       opts,
       PlotRange -> plotrange, PlotRangePadding -> padding, 
       Ticks -> None, ExclusionsStyle -> None, 
       PlotLabel -> First@eqns] /. 
      l_Line :> {Arrowheads[
         PadRight[{0.}, 1 + Ceiling[2/r2*ArcLength@l], 0.04]], 
        Arrow @@ l}
     ]
    ]
   ];
Improve this answer
$\endgroup$
2
  • $\begingroup$ Very nice, Indeed! :) $\endgroup$
    – E. Chan-López
    May 26, 2022 at 23:11
  • $\begingroup$ Amazing thank you! $\endgroup$
    – NC520
    May 27, 2022 at 11:58

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.