10
$\begingroup$

I am trying to draw this system of differentiable equations,

enter image description here

whose phase diagram is:

enter image description here

found in the following article. I am new using Mathematica, I was guided by these posts: 1 and 2. But I get an error every time so I gave up. Any ideas please?

$\endgroup$
2
  • 7
    $\begingroup$ Please include the code you've tried already, and the error it generates. $\endgroup$
    – MarcoB
    Jul 28 at 2:10
  • 1
    $\begingroup$ BTW, one can get a symbolic solution, after some coaxing: {{r -> Function[t, 1 + (r0 - 1)/(1 + 4 (-1 + r0)^4 t)^(1/4)], \[Theta] -> Function[{t}, 1/3 (3 \[Theta]0 - 2 Cos[3/(2 (-1 + r0))]) + 2/3 Cos[(3 (1 + 4 (-1 + r0)^4 t)^(1/4))/(2 (-1 + r0))]]}}. It might be helpful for further analysis. $\endgroup$
    – Michael E2
    Jul 29 at 4:15

4 Answers 4

18
$\begingroup$

Since we're giving answers...

ndsol = ParametricNDSolveValue[{r'[t] == -(r[t] - 1)^5, θ'[t] ==
      Piecewise[{{-(r[t] - 1)^3 Sin[3/(2 (r[t] - 1))], r[t] != 1}}],
    r[0] == r0, θ[0] == θ0},
   {r, θ}, {t, 0, 2000}, {r0, θ0}];

polarTrajectory[
   ndsolution : {_InterpolatingFunction, _InterpolatingFunction}] := 
#[[1]] Transpose@Through[{Cos, Sin}[#[[2]]]] &@
   Through[ndsolution["ValuesOnGrid"]];

startingPoints = {{3, Pi/2}, {0.01, 0}, {2, Pi}};
Show[
 Graphics@{Red, Thick, Circle[]},
 ListLinePlot[
  polarTrajectory /@ ndsol @@@ startingPoints,
  InterpolationOrder -> 3, AspectRatio -> Automatic
  ]
 , PlotRange -> All, Frame -> True
 ]

Update: Response to comment. With arrowheads at about the places shown in the OP:

startingPoints = {{3, Pi/2}, {0.01, 0}, {2, Pi}};
Show[
 Graphics@{Red, Thick, Circle[]},
 ListLinePlot[
   polarTrajectory /@ ndsol @@@ startingPoints,
   InterpolationOrder -> 3, AspectRatio -> Automatic
   ] /. Line[
    p_] :> {Arrowheads[{{.04, 
       Abs[Norm[First[p]] - 1]/2/ArcLength@Line@p}}], Arrow[p]}
 , PlotRange -> All, Frame -> True
 ]

The scaled position Abs[Norm@First@p - 1]/2/ArcLength@Line@p is relative to the length of the trajectory: it won't change perceptibly if the length of the integration is changed (to, say, {t, 0, 10^6} as suggested in a comment). The arrowhead will be placed along the trajectory at an arc length equal to half the distance from the starting point to the circle.

Addendum

Here is an exact, symbolic solution (spoils the challenge in the comment, tho'):

exactsol = {
  {r -> Function[t, 
      1 + (r0 - 1)/(1 + 4 (-1 + r0)^4 t)^(1/4)],
   θ -> Function[{t}, 
      1/3 (3 θ0 - 2 Cos[3/(2 (-1 + r0))]) + 
       2/3 Cos[(3 (1 + 4 (-1 + r0)^4 t)^(1/4))/(2 (-1 + r0))]]}};
ParametricPlot[
  r[t] {Cos[θ[t]], Sin[θ[t]]} /. exactsol /. 
      Thread[{r0, θ0} -> #] & /@ startingPoints // Evaluate,
  {t, 0, 10^6}, PlotPoints -> 200, MaxRecursion -> 15, 
  Prolog -> {Red, Thick, Circle[]}
  ] /.
 Line[p_] :> {Arrowheads[{{.04, 
      Abs[Norm[First[p]] - 1]/2/ArcLength@Line@p}}], Arrow[p]}
$\endgroup$
9
  • $\begingroup$ Running the integration out to t = 10^6 illustrates the point in the caption to the OP figure better: i.stack.imgur.com/xQroV.png $\endgroup$
    – Michael E2
    Jul 28 at 4:33
  • $\begingroup$ Oh, you did it the easy way, just solved the ode's :) $\endgroup$
    – Nasser
    Jul 28 at 4:39
  • $\begingroup$ I see nothing wrong with f1, f2. I used same equations posted. Just converted $r$ to $x,y$ coordinates. I think may be stream plot sampling not hitting on the correct stream lines which will show the same ones posted. If you see something wrong please let me know what it is. I think it is hard to get same plot with SteamPlot, since the steam lines are very sensitive to where they start at for this system. $\endgroup$
    – Nasser
    Jul 28 at 4:43
  • 1
    $\begingroup$ Here's an exercise for the reader: Each (nontrivial) trajectory approaches approaches an arc of length 4/3 radians — true or false? $\endgroup$
    – Michael E2
    Jul 28 at 4:54
  • 1
    $\begingroup$ @SigisK Piecewise[f, default] has a default value if none of the conditions are met. The default default value is 0. The condition is there but hidden. Evaluate Piecewise[{{-(r[t] - 1)^3 Sin[3/(2 (r[t] - 1))], r[t] != 1}}] and examine the output; or examine its InputForm. Add the piece {0, r[t] == 0} to Piecewise after the {.., r[t] != 1} case and examine the InputForm again. You will see, in this simple case, that Piecewise eliminated it, since it's redundant. Add it in front of the r[t] != 0 case, and it does not eliminate it. So autosimplification is limited. $\endgroup$
    – Michael E2
    Jul 30 at 14:13
8
$\begingroup$

** corrected transformation of coordinate**

Thanks to Michael E2 and user293787 my earlier plot did not match well, as I had error in transformation from polar to $x,y$.

If you give the fixed points, then one can make the specific streamlines you show. I am not going to read the paper to find these.

But here is overall phase plot, which will contain many stream lines. The one you show has specific ones that pass through fixed points it seems. Once you know these, then StreamPlot can plot only those you specificed.

tf = TransformedField[
   "Polar" -> 
    "Cartesian", {(-(r - 1)^5), -(r - 1)^3*Sin[3/(2*(r - 1))]}, {r, 
     theta} -> {x, y}];
f1[x_, y_] := tf[[1]];
f2[x_, y_] := If[Sqrt[x^2 + y^2] == 1, 0, tf[[2]]];


p0 = {Red, PointSize[0.03], Point[{0, 0}]};
c = {Red, Circle[{0, 0}, 1]};

StreamPlot[Evaluate[{f1[x, y], f2[x, y]}], {x, -2, 2}, {y, -2, 2}, 
 Epilog -> {p0, c}, ImageSize -> 300, PerformanceGoal -> "Quality"]

Mathematica graphics

To plot stream lines that only passes through specific point, you do (as an example, change the point as you need

StreamPlot[Evaluate[{f1[x, y], f2[x, y]}], {x, -2, 2}, {y, -2, 2},
 StreamPoints -> {{{{.7, .7}, Red}, 1}},
 Epilog -> {p0, c}, ImageSize -> 300, PerformanceGoal -> "Quality"]

Mathematica graphics

This is a plot for larger range

StreamPlot[Evaluate[{f1[x, y], f2[x, y]}], {x, -4, 4}, {y, -4, 4},
 Epilog -> {p0, c}, ImageSize -> 300, PerformanceGoal -> "Quality"]

Mathematica graphics

The above shows more clearly the unit circle is a limit circle (i.e. solutions that start outside it, remain outside, and solutions that start from inside it, remain inside it.

$\endgroup$
2
  • 1
    $\begingroup$ @user293787 Yes, you are right, I was just looking at this now. Thanks. I will correct this. Let me close this for now and fix it :) $\endgroup$
    – Nasser
    Jul 28 at 4:50
  • $\begingroup$ Thank you very much for your detailed answer, it was great! $\endgroup$
    – Zaragosa
    Jul 28 at 18:42
7
$\begingroup$

For the beginner,maybe the easy way is use ToPolarCoordinates and ParametricPlot. ( We can change the initial point {x0, y0} = {0.1, 4}. to anothers points)

Clear[x0, y0, r0, θ0, sol];
{x0, y0} = {0.1, 4};
{r0, θ0} = ToPolarCoordinates[{x0, y0}];
sol = NDSolve[{r'[t] == -(r[t] - 1)^5, θ'[t] == 
     If[r[t] != 1, -(r[t] - 1)^3 Sin[3/(2 (r[t] - 1))], 0], 
    r[0] == r0, θ[0] == θ0}, {r[t], θ[t]}, {t, 0,
     2000}];
ParametricPlot[
 r[t] {Cos[θ[t]], Sin[θ[t]]} /. sol[[1]], {t, 0, 2000}, 
 PlotPoints -> 200, MaxRecursion -> 4, AspectRatio -> Automatic, 
 Prolog -> {{Red, Circle[]}}, PlotRange -> 4]

enter image description here

$\endgroup$
1
  • $\begingroup$ Evaluate[FromPolarCoordinates[{r[t], θ[t]}] /. sol[[1]]] $\endgroup$
    – cvgmt
    Jul 28 at 5:53
0
$\begingroup$

You can also find an analytical solution to this ODE system:

DSolve[{r'[t] == -(r[t] - 1)^5, \[Theta]'[t] == Piecewise[{{0, r[t] == 1}, {-(r[t] - 1)^3 Sin[3/(2 (r[t] - 1))], r[t] != 1}}]}, {r, \[Theta]}, t] // Simplify

enter image description here

Pay attention to the fact that the system has 4 general solutions, but the particular solution (depending on the initial conditions) will be only 1.

It is this solution that is built on the phase diagram.

General solutions may be needed for a deeper study of the properties of the ODE's system under consideration.

$\endgroup$
1
  • $\begingroup$ One can find the general solution at the end of my answer. (The process of reducing, selecting, and simplifying the DSolve solutions was more of a pain than using inspection and one's intelligence. It is easier to solve & reduce the r equation and then plug the solution into the θ equation and solve that. Still a bit of a pain, though.) $\endgroup$
    – Michael E2
    Jul 29 at 17:08

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.