Plot a phase diagram

I am trying to draw this system of differentiable equations,

whose phase diagram is:

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?

• Please include the code you've tried already, and the error it generates. Jul 28, 2022 at 2:10
• 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. Jul 29, 2022 at 4:15

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[
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.

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[]}
] /.
Abs[Norm[First[p]] - 1]/2/ArcLength@Line@p}}], Arrow[p]}

• 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 Jul 28, 2022 at 4:33
• Oh, you did it the easy way, just solved the ode's :) Jul 28, 2022 at 4:39
• 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. Jul 28, 2022 at 4:43
• Here's an exercise for the reader: Each (nontrivial) trajectory approaches approaches an arc of length 4/3 radians — true or false? Jul 28, 2022 at 4:54
• @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. Jul 30, 2022 at 14:13

** 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"]


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


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


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.

• @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 :) Jul 28, 2022 at 4:50
• Thank you very much for your detailed answer, it was great! Jul 28, 2022 at 18:42

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]


• Evaluate[FromPolarCoordinates[{r[t], θ[t]}] /. sol[[1]]] Jul 28, 2022 at 5:53

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


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.

• 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.) Jul 29, 2022 at 17:08