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

Question

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

Answer 1

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] &) ]]

Answer 2

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]

Answer 3

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