# How to plot the phase portrait of this non-linear system of ODEs

I am trying to plot the phase portrait of the following differential equation $$\begin{cases} \dot{R}= aR+bJ\vert1-J\vert+f(t)\\ \dot{J}= cR\vert1-R\vert+dJ, \end{cases} \quad f(t)=5\sin(\pi t).$$

for $$a=-7,\,b=-2,\,c=1,\,d=1$$. The plot should look something like

but using the code

s[a_, b_] := a*R[t] + b*J[t]*Abs[1 - J[t]] + 5 Sin[\[Pi]*t];
h[c_, d_] := c*R[t]*Abs[1 - R[t]] + d*J[t];
StreamPlot[{s[-7, -2], h[1, 1]}, {R[t], 0, 5}, {J[t], 0, 5}]


does not give the desired result.

Setting the paramters to be a = −1.1, b = −2, c = 1, and d = 1 should give chaotic behavior, as shown below

• From where you get these equations and graphs? Any weblink?
– zhk
Commented Dec 6, 2018 at 5:57
• @zhk this paper Commented Dec 6, 2018 at 6:00

Perhaps you can use ParametricNDSolveValue and ParametricPlot?

ClearAll[a, b, c, d, pndsv]
pndsv = ParametricNDSolveValue[{R'[t] == a*R[t] + b*J[t]*Abs[1 - J[t]] + 5 Sin[π*t],
J'[t] == c*R[t]*Abs[1 - R[t]] + d*J[t], R[0] == r, J[0] == j}, {R, J},
{t, 0, 100}, {a, b, c, d, r, j}];

Manipulate[ParametricPlot[Evaluate[Through@pndsv[a, b, c, d, r, j][t]], {t, 0, tmax},
AspectRatio -> 1, Frame -> True, Axes -> False,
PlotRange -> {{-10, 10}, {-10, 10}}],
{{a, -1.1}, -10, 10}, {{b, -2}, -10, 10}, {{c, 1}, -10, 10}, {{d, 1}, -10,  10},
{{r, 1}, 0, 1}, {j, 0, 1}, {{tmax, 10}, 1, 200}]


ParametricPlot[Evaluate[Join @@ Table[Through@pndsv[-7, -2, 1, 1, r, j][t],
{r, 0, 1, .25}, {j, 0, 1, .25}]], {t, 0, 10},
AspectRatio -> 1, Frame -> True, Axes -> False,
PlotRange -> {{-3, 3}, {-5, 5}}, FrameLabel -> {{J, None}, {R, None}},
PlotLegends -> LineLegend[97, Join @@ Table[{r, j}, {r, 0, 1, .25}, {j, 0, 1, .25}],
LegendLayout -> {"Column", 2},
LegendFunction -> (Labeled[Panel[#], Style["{R[0], J[0]}", 16], Top] &) ]]


• is there a way to make it start at (1,0) as shown in the picture? Commented Dec 6, 2018 at 4:43
• $a=-1.1$ and $b=-2$ should produce chaotic behavior, however it does not seem so from your animation Commented Dec 6, 2018 at 4:49
• @DMH16, if you set r=1 and j=0 the trajectory starts at {1,0}; but figuring out the exact configuration of {a,b,c,d} to get the picture you posted is a challenge.
– kglr
Commented Dec 6, 2018 at 4:50
• @DMH16, try increasing the bounds on t in ParametricNDSolveValue and ParametricPlot, e.g. {t, 0,500}.
– kglr
Commented Dec 6, 2018 at 4:55
• I tried but it doesn't show what I am looking for. Can you take a look at the addendum? Commented Dec 6, 2018 at 5:12

StreamPlot is not producing anything because you have Sin[Pi*t].

Here is something for you to experiment with,

a = -7; b = -2; c = 1; d = 1;

sol[R0_?NumericQ] :=
First@NDSolve[{R'[t] == a*R[t] + b*J[t]*Abs[1 - J[t]] + 5 Sin[\[Pi]*t],
J'[t] == c*R[t]*Abs[1 - R[t]] + d*J[t], R[0] == R0,J[0] == R0}, {R, J}, {t, 0, 10}];

ParametricPlot[Evaluate[{R[t], J[t]} /. sol[#] & /@ Range[-15, 15, 1]], {t, 0, 10},
Frame -> True]


I think we will not be able to accurately reproduce the data from Youngchul Bae, Nonlinear Behavior in Love Model with Discontinuous External Force, Int. J. Fuzzy Log. Intell. Syst. 2016;16(1):64-71., since the author does not provide details about f[t]= 5 Sin[π*t1], only that it is a discontinuous external force. I picked up a suitable function, but it does not allow to reproduce the solution in all its details.

s[a_, b_, t_, R_, J_] := a*R + b*J*(1 - J) + 5*Sin[2*\[Pi]*f[t]];
h[c_, d_, t_, R_, J_] := c*R*(1 - R) + d*J;
f[t_] := Piecewise[
Flatten[Table[{{t, 10*n <= t <= 10*n + 2}, {0,
10*n + 2 <= t <= 10*(n + 1)}}, {n, 0, 10}], 1]]

sol = NDSolve[{R'[t] == s[-7, -2, t, R[t], J[t]],
J'[t] == h[1, 1, t, R[t], J[t]], R[0] == 1, J[0] == 1}, {R, J}, {t,
0, 100}]

{Grid[{{"a=", -7}, {"b=", -2}, {"c=", 1}, {"d=", 1}}],
Plot[Evaluate[{R[t], J[t]} /. sol], {t, 0, 100},
PlotLegends -> {"R", "J"}],
ParametricPlot[{R[t], J[t]} /. sol, {t, 0, 100}, PlotRange -> All,
AspectRatio -> 1, Frame -> True, FrameLabel -> {"R", "J"},
Axes -> False]}