2
$\begingroup$

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 enter image description here

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

enter image description here

|improve this question
$\endgroup$
  • $\begingroup$ From where you get these equations and graphs? Any weblink? $\endgroup$ – zhk Dec 6 '18 at 5:57
  • $\begingroup$ @zhk this paper $\endgroup$ – DMH16 Dec 6 '18 at 6:00
4
$\begingroup$

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

enter image description here

enter image description here

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

enter image description here

|improve this answer
$\endgroup$
  • $\begingroup$ is there a way to make it start at (1,0) as shown in the picture? $\endgroup$ – DMH16 Dec 6 '18 at 4:43
  • $\begingroup$ $a=-1.1$ and $b=-2$ should produce chaotic behavior, however it does not seem so from your animation $\endgroup$ – DMH16 Dec 6 '18 at 4:49
  • $\begingroup$ @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. $\endgroup$ – kglr Dec 6 '18 at 4:50
  • $\begingroup$ @DMH16, try increasing the bounds on t in ParametricNDSolveValue and ParametricPlot, e.g. {t, 0,500}. $\endgroup$ – kglr Dec 6 '18 at 4:55
  • $\begingroup$ I tried but it doesn't show what I am looking for. Can you take a look at the addendum? $\endgroup$ – DMH16 Dec 6 '18 at 5:12
2
$\begingroup$

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]

enter image description here

|improve this answer
$\endgroup$
1
$\begingroup$

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

fig1 fig2

|improve this answer
$\endgroup$

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.