# Draw phase line associated to slope field

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}}]


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

• – Moo
May 26, 2022 at 14:53
• Yes the second ink does the job. I was just hoping to find a more "automated" process. May 26, 2022 at 16:24
• This one is similar to what you are looking for: mathematica.stackexchange.com/a/115676/731 May 26, 2022 at 19:58

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}]


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 = {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*),
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}]


(*
* 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] /; (
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;
Message[diffEq::CantHappen,"Plotting",Row[{"the 1D vector field is \
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
, 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,
AxesLabel -> {HoldForm[Globalv], 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,
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,
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[Globalv]}, {"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,