# Plotting orbit of a discontinuous map

I'm trying to plot the orbit of a two dymensional discontinuous map. Could you help me please? I tried using Piecewise, but I didn't get the result I wanted. The map has this form:

$$\begin{pmatrix}x'\\y'\end{pmatrix}= \begin{pmatrix}a & -c \\1 & 0\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}+0.8 n\begin{pmatrix}1\\0\end{pmatrix} \quad \text{if} \quad (x,y)\in R_{1}$$ $$\begin{pmatrix}x'\\y'\end{pmatrix}=\begin{pmatrix}d & c \\1 & 0\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}+1.3 n\begin{pmatrix}1\\0\end{pmatrix}\quad \text{if} \quad (x,y)\in R_{2}$$ where $R_{1}$ and $R_{2}$ are : \begin{aligned}R_{1}&= \{-x_{t}\leq y_{t} \leq x_{t}\} &\cup & \quad \{x_{t}\leq y_{t} \leq-x_{t}\}\\R_{2}&=\begin{cases}y_{t}>x_{t} \\y_{t}>-x_{t}\end{cases} &\cup & \quad \begin{cases}y_{t}<-x_{t} \\y_{t}<x_{t}\end{cases} \end{aligned}

This is the code I used to define the map

a = 0.6; c = 1; d = 0.8; n = 0.4; x = 1.2; y = 1.3;
z[{x_, y_}] := {a x - c y + 0.8 n , x}
g[{x_, y_}] := {d x + c y + 1.3 n, x}
A = NestList[z, {1.25, 1.3}, 1500];
B = NestList[g, {1.25, 1.3}, 1500];
If[-x[t] <= y[t] <= x[t] || x[t] <= y[t] <= -x[t], ListPlot[A],
ListPlot[B]]


but it returns the following output

If[-x[t] <= y[t] <= x[t] || x[t] <= y[t] <= -x[t], ListPlot[A], ListPlot[B]]

• What values do the constants have? Also, please provide your equations in Mathematica format. – bbgodfrey Aug 3 '17 at 16:52
• Since you say "map", we can assume that the primes denote $x_{t+1}$, not $dx/dt$? – Chris K Aug 3 '17 at 17:15
• @bbgodfrey Thank you for your comments. I have updated the question with the mathematica code I have used. It also includes the values of the constants. – smashing6 Aug 3 '17 at 17:58
• Running your code: A is fine, ListPlot[A] plots as an ellipse. But B is unstable and diverges. You can plot the first few and watch it diverge using ListPlot[B[[1 ;; 30]]]. – bill s Aug 3 '17 at 17:59

I interpret the map to be

z[{x_, y_}] := Piecewise[{{{a x - c y + 0.8 n , x}, -Abs[x] <= y <= Abs[x]}},
{d x + c y + 1.3 n, x}]


in which case the plot is

ListPlot@NestList[z, {1.25, 1.3}, 1500] Is this what you had in mind?

Here is a similar map that does not blow up.

z[{x_, y_}] := Piecewise[{{{a x - c y + 0.8 n , x}, -1/5 Abs[x] <= y <= 5 Abs[x]}},
{d x + c y + 1.3 n, x}]
ListPlot[NestList[z, {1.25, 1.3}, 1500], AxesLabel -> {x, y}] • Yes, this is what I meant. But in your definition of the map region 2 is missing. Is there also a way to plot the orbit over time? – smashing6 Aug 3 '17 at 18:31
• @smashing6 Map region 2 is given by {d x + c y + 1.3 n, x} in Piecewise. – bbgodfrey Aug 3 '17 at 18:33
• I think that my problem derived from a bad usage of function Piecewise. This is really what I meant. Thank you very much. If you don't mind, I want to ask you if it is possible to plot the basins of attraction of this map. I've tried to use the ListDensityPlot without results. – smashing6 Aug 3 '17 at 18:56
• @smashing6 So that I am sure I understand, are you asking whether it is possible to plot all initial conditions that lead to, for instance, the second plot? If so, the answer is yes. See, for instance, 101400. Note, however, that the present problem would be more difficult, because the attractor is 2D rather than 1D. – bbgodfrey Aug 3 '17 at 19:00
• Yes, I agree with you. I'll follow the page you have suggested me. I hope to plot something similar. Anyway, thank you for your helpful suggestions. – smashing6 Aug 3 '17 at 19:16

Considerable insight can be gained by solving parts of the question symbolically before using numerical analysis. Each of the two recursions can be treated as second-order linear recurrence relations, which for generality can be written as

eq = {x[m + 1] == c1 x[m] + c2 x[m - 1] + c3;


and solved by

s = x[m] /. Collect[Flatten@RSolve[{eq, x == x0, x == x1}, x[m], m],
{(c1 + Sqrt[c1^2 + 4 c2])^m, (c1 - Sqrt[c1^2 + 4 c2])^m},
FullSimplify[#, (c1 | c2 | c3) ∈ Reals && m ∈ Integers] &]
(* -(c3/(-1 + c1 + c2)) + (2^(-1 - m) (c1 - Sqrt[c1^2 + 4 c2])^ m
((-2 + c1 + Sqrt[c1^2 + 4 c2]) c3 + (-1 + c1 + c2) ((c1 + Sqrt[c1^2 + 4 c2]) x0 - 2 x1)))
/((-1 + c1 + c2) Sqrt[c1^2 + 4 c2]) + (2^(-1 - m) (c1 + Sqrt[c1^2 + 4 c2])^ m
((2 - c1 + Sqrt[c1^2 + 4 c2]) c3 + (-1 + c1 + c2) (-c1 x0 +
Sqrt[c1^2 + 4 c2] x0 + 2 x1)))/((-1 + c1 + c2) Sqrt[c1^2 + 4 c2]) *)


with a fixed point,

fp = Flatten@Solve[{Coefficient[s, (c1 + Sqrt[c1^2 + 4 c2])^m ] == 0,
Coefficient[s, (c1 - Sqrt[c1^2 + 4 c2])^m ] == 0}, {x0, x1}] // Simplify
(* {x0 -> -(c3/(-1 + c1 + c2)), x1 -> -(c3/(-1 + c1 + c2))} *)


Now, insert parameters from the first recurrence relation.

a = 3/5; c = 1; d = 4/5; n = 2/5;
s /. {c1 -> a, c2 -> -c, c3 -> 4 n/5};
s1 = Collect[ComplexExpand[%, TargetFunctions -> {Re, Im}], {_Sin, _Cos}, Simplify]
(* 8/35 + (-(8/35) + x0) Cos[m ArcTan[Sqrt/3]] -
((8 + 15 x0 - 50 x1) Sin[m ArcTan[Sqrt/3]])/(5 Sqrt) *)
fp /. {c1 -> a, c2 -> -c, c3 -> 4 n/5}
(* {x0 -> 8/35, x1 -> 8/35} *)


Based on the expression above, curves of x[m-1] vs x[m] will resemble ellipses, encircling the fixed point just evaluated. For the second recurrence relation,

s2 = s /. {c1 -> d, c2 -> c, c3 -> 13 n/10}
(* -(13/20) + (25 2^(-4 - m) (4/5 - (2 Sqrt)/5)^m (13/25 (-(6/5) + (2 Sqrt)/5) +
4/5 ((4/5 + (2 Sqrt)/5) x0 - 2 x1)))/Sqrt + (25 2^(-4 - m)
(4/5 + (2 Sqrt)/5)^m (13/25 (6/5 + (2 Sqrt)/5) + 4/5 (-((4 x0)/5)
+ (2 Sqrt x0)/5 + 2 x1)))/Sqrt *)


Parhaps, it could be simplified in terms of hyperbolic functions, but FullSimplify seems unable to do so. Instead, consider the m-dependent quantities.

Cases[s2, z_^m -> z/2, Infinity] // Simplify
(* {1/5 (2 - Sqrt), 1/5 (2 + Sqrt)} *)
% // N
(* {-0.677033, 1.47703} *)


Consequently, s2 can be expected to grow exponentially. For completeness,

fp /. {c1 -> d, c2 -> c, c3 -> 13 n/10}
(* {x0 -> -(13/20), x1 -> -(13/20)} *)


A qualitative picture of the behavior of the combined recurrence relation given in the question now emerges. Concentric closed curves can exist entirely within the first region. However, a partial curve crossing the border into the second region will grow exponentially in magnitude, possible reentering the first region but more typically growing arbitrarily large. However, for the particular boundaries given in the question, -Abs[x] <= y <= Abs[x], the fixed point for the first region lies on the boundary, indicating that no closed curves are likely to exist. I have performed a thorough scan of {x1, x0} and found none. Instead, all solutions move toward infinity except those at the fixed points.

Consider, instead, the boundaries, -1/5 Abs[x] <= y <= 5 Abs[x], which do not pass through the stable fixed point. Consequently, a small region of closed solutions can exist. For instance,

divides = ContourPlot[{y == 5 Abs[x], y == -Abs[x]/5}, {x, -2, 2}, {y, -2, 2},
ContourStyle -> Black, FrameLabel -> {x, y},
FrameStyle -> Directive[12, Bold, Black], ImageSize -> Medium];
Show[divides,
ListPlot[Partition[Table[s1 /. {x0 -> .373, x1 -> .373}, {m, 2, 500}], 2, 1]]] Shown are the largest closed curve existing entirely within the first region, and the boundaries between the regions. Curves crossing boundaries must be determined numerically. Below is a plot of initial conditions that lead to closed curves.

zz[{x_, y_}] := Piecewise[{{{a x - c y + 0.8 n , x}, -1/5 Abs[x] <= y <= 5 Abs[x]}},
{d x + c y + 1.3 n, x}]
(dtaa = ParallelTable[{x, y, Length@Union[lst = NestWhileList[zz, {x, y},
Norm[#] < 1000 &, 1, 30000]; If[Length@lst < 1000, {}, lst[[29000 ;;]]],
SameTest -> (Norm[#1 - #2] < .6 &)]}, {y, -2, 2, .02}, {x, -2, 2, .02}];)
Show[divides, ArrayPlot[Map[Last, dtaa, {2}], ColorFunction -> "Rainbow",
DataReversed -> {True, False}, DataRange -> {{-2, 2},{-2, 2}}], FrameLabel -> {x, y},
Epilog -> {PointSize[Medium], Point[{{8/35, 8/35}, {-13/20, 13/20}}]},
FrameStyle -> Directive[12, Bold, Black], ImageSize -> Large] All curves originating from the red areas eventually find themselves in rougly elliptical orbits within the center-most red area with the stable fixed point (a black dot) at its center. (For completeness, the unstable fixed point for the second recurrence relation also is shown as a black dot.) Solutions not originating from red areas grow exponentially in magnitude and do not form closed curves.