Discrete-time dynamical system orbit diagram

This article by Dan Polansky complements the Wikipedia article on "bifurcation diagram". The term to be used in the title was chosen to be "orbit diagram" (rather than "bifurcation diagram"), used by multiple decent sources and fitting the thing plotted well, the values visited by the orbits. This article is concerned only with discrete-time dynamical systems; the topic of bifurcation diagrams could cover continuous-time dynamical systems as well, but whether it really does needs more research.

We will show images of orbit diagrams for various one-dimensional real-valued maps, where map is a function that is to be iterated, that is, applied repeatedly to its previous output. The map defines the transition function of the dynamical system, that is, it defines the next state of the system given the current state. One kind of orbit diagram plots, for each value of a parameter of the system, the values visited by the orbit starting from a fixed initial value. A more typical orbit diagram is a modification of this: instead of plotting all visited values, it only plots a certain number of last visited values; therefore, if the dynamical system quickly converges near one or more attracting values, the result is more concentrated around these.

An orbit diagram often shows shades of gray rather than black and white. This tracks the density of points in the pixel: a pixel (a square) in which more orbit points lie gets a darker color.

An orbit is the sequence of values x₀, f(x₀), f(f(x₀)), f(f(f(x₀)))), etc. for a given map f, where the map is parametrized with one real-valued parameter. A fixed initial value of x₀ is often chosen to plot the diagram. Alternatively, one can choose the initial value randomly in a certain range; thus for the logistic map, one could choose x₀ as random value in [0, 1].

An orbit diagram is usually presented as a raster image, e.g. a PNG. However, one can also think of the orbit diagram as an infinitely precise mathematical object of which the PNG is a mere rasterization. Such a diagram is a set of pairs (c, x) such that the orbit for parameter c (or lambda, etc.) visits the value x within first n iterations. This diagram is a union of n curves that represent the first n iterates. Thus, each point in the diagram lies on at least one curve (iterate) that brought it into the diagram rather than being isolated. Moreover, one can consider a diagram (a mathematical object rather than rasterized plot) that contains n-iterates for all n rather than only for n up to a limit.

Along with the orbit diagram, we also show cobweb plots where available since these show the shape of the map, and remind how the iteration works and impacts values.

Universality

Having a look at orbit diagrams of these dynamical systems reveals certain common features, e.g. bifurcating points at which the set of values attracted to doubles, and regions of chaotic orbits or orbits with very high size (even if finite) of the sets of attracting values.

Logistic map

The logistic map is f(x) = r * x * (1 - x).^[1] The typically plotted range of r values is [0, 4]. The initial x value is e.g. 0.5.

Orbit diagrams of the logistic map:

Cobweb plot:

Orbit diagram showing negative values of r as well:

Orbit diagram zoomed on the leftmost bifurcation period-doubling region:

Orbit diagram zoomed on the period-3 window:

Further reading:

Logistic map, wikipedia.org
Logicedges - Bifurcation diagrams and the Feigenbaum constant, sites.google.com

Mandelbrot set quadratic map

The quadratic map used in the Mandelbrot set is f(x) = x**2 + c.^[2] The typically plotted range of c values is [-2, 0.5]. The initial x value is e.g. 0.

Orbit diagrams of this quadratic map:

Critical curves of this quadratic map:

The above diagram lets one realize that the orbit diagram is but a union of what the above diagram calls "critical curves".

Further reading:

Commons:Category:Real quadratic map
The Road to Chaos is Filled with Polynomial Curves by Richard D. Neidinger and R. John Annen, The American Mathematical Monthly, 1996

Tent map

The tent map is $f_{\mu }(x):=\mu \min\{x,\,1-x\}$ . (MathWorld uses the formula mu * (1 - 2 * abs(x - 0.5)).^[3])

Orbit diagram of this map:

Cobweb diagram:

Orbit diagram of this map for a broader range of parameter values:

Further reading:

Tent map, wikipedia.org

Sine map 1

Let the term "sine map" denote f(x) = c * sin(pi * x). The function's overall shape on the [0, 1] interval is similar to the logistic function. The initial value can be taken to be 0.5; if one takes 0, the system will stay at 0 regardless of the value of c.

Orbit diagram of this map:

Orbit diagram of this map on a broader set of parameter values:

Orbit diagram of this map with the initial value being near zero yet non-zero:

Further reading:

Chaos from maps, physics.ox.ac.uk
Feigenbaum Constant, mathworld.wolfram.com -- mentions and defines the "sine map"

Sine map 2

We considered one parametrization of the sine function as a "sine map" above. However, one can take a different formula, f(x) = sin(c * pi * x), and this can also reasonably be called "sine map". Since the range of the function is kept in [-1, 1], plotting the orbit diagram for this formula on a larger interval of the parameter c does not lead to linear stretching of the range of orbit values.

Orbit diagram of this map:

Orbit diagram of this map on a broader set of parameter values:

Further reading:

Dynamics, Chaos, and Fractals (part 2): Dynamics of One-Parameter Families by Evan Dummit, 2015, northeastern.edu

Circle map

The circle map is f(x) = (x + omega - K / (2 * pi) * sin(2 * pi * x)) mod 1. The mod 1 means that the values of x represent points on a circle by the angles pointing to the points; the values are not in radians but in such units that 1 makes the complete circle. (In radians, 2 * pi would make the complete circle.)

Orbit diagram of this map:

Further reading:

Circle Map, mathworld.wolfram.com
Arnold tongue, wikipedia.org

Tanh map

Orbit diagram of this map:

The diagram documentation does not exactly identify the map; tanh (hyperbolic tangent) probably plays a role.

Ricker map

Alleged orbit diagram of this alleged map:

The term "Ricker map" is based on the description of File:Ricker bifurcation.png and does not trace to a reliable source at this point.

Bell et al. 2010 use the term "Ricker's population model" to refer to the formula x_n+1 = x_n^delta * e^r-x_n−k.^[4] It is not clear how to obtain a map from this formula since x_n+1 depends not only on x_n but also on x_n-k. However, in later slides, Bell et al. 2010 use the formula x_n+1 = x_n^delta * e^r-x_n (the "-k" part is gone), which is amenable to be interpreted as a map and to have an orbit diagram plotted as long as one chooses one of the two parameters delta and r to be fixed, and the other to be varied along the x-axis.

Wikipedia gives the formula $N_{t+1}=N_{t}e^{r\left(1-{\frac {N_{t}}{k}}\right)}.\,$ for what is calls "Ricker model". File:Ricker bifurcation.png does not state any formula or plotting code. To plot an orbit diagram, we need to choose one parameter to be varied along the x-axis, and fix any other parameter to a chosen value; the formula from Wikipedia features parameters r and k. The plot lets r vary along the x-axis; the Wikipedia caption for the image states "Bifurcation diagram of the Ricker model with carrying capacity of 1000", which could mean that k is the carrying capacity and was set to 1000. This needs clarification. As an exercise, one can try to reproduce the plot using the Python code in this article as a starting point to modify using the information just presented.

Further reading:

Ricker model, wikipedia.org

Other maps

Orbit diagrams of other one-dimensional maps studied in the subject of dynamical systems can be found in the linked Commons category. Other orbit diagrams can be plotted: one can take the Python code from the following section, and change the iterative function and plotting parameters.

Plotting in Python

The following code plots an orbit diagram for the Mandelbrot set quadratic map (adapted from Wikipedia article on the logistic function). Tweaking the code should lead to other maps being plotted, by changing the iterative function and other parameters. The code plots the orbits using individual pixels (each pixel hit by a point is dark) and the result is two-colored rather than greyscale.

import numpy as np
import matplotlib.pyplot as plt

saveFigure = True
markerSize = 0.001 # orig. 0.02
markerStyle = "," # . - point; , - pixel
preset = 1
if preset == 1: # Complete diagram
  cInterval = (-2.1, 0.26)
  #cInterval = (-1.4011551890, -1.3940462)
  cResolution = 0.0003
  mapIterationCount = 1000  # (orig. 600)
  lastMapIterationCountToPlot = 1000 # (orig. 200)
  startingIterationValue = 0
  fileName = "MandelbrotMapOrbitDiag.png"
elif preset == 2: # Zoom to a region near -2
  cInterval = (-2, -1.99)
  cResolution = 1E-6
  mapIterationCount = 1000
  lastMapIterationCountToPlot = 1000
  startingIterationValue = 0
  fileName = "MandelbrotMapOrbitDiag_ZoomNearM2.png"

if cInterval[0] < -2 or cInterval[1] > 0.25:
  yLim = (-2.2, 2.5)
else:
  yLim = (-2, 2)

lims = np.zeros(mapIterationCount)

fig, biax = plt.subplots()
fig.set_size_inches(16, 9)

lims[0] = startingIterationValue
for c in np.arange(cInterval[0], cInterval[1], cResolution):
  for i in xrange(mapIterationCount - 1):
    lims[i + 1] = lims[i] * lims[i] + c
    if lims[i] > 2:
        lims[i + 1] = lims[i] # Allow plotting of c values outside of Mandelbrot set, diverging ones

  biax.plot([c] * lastMapIterationCountToPlot,
            lims[(mapIterationCount - lastMapIterationCountToPlot):], "b.",
            markersize=markerSize,
            marker=markerStyle)

biax.set(xlabel="c", ylabel="x",
         title="Mandelbrot set quadratic map for real values - bifurcation diagram/orbit diagram" + "\n" +
               "c resolution: " + str(cResolution) +
               ", iteration count: " + str(mapIterationCount) +
               ", last values to plot count: " + str(lastMapIterationCountToPlot))

biax.set_ylim(yLim[0], yLim[1])
if saveFigure:
  plt.savefig(fileName, dpi=600)
else:
  plt.show()

The most performance intensive part seems to be the plotting, not the iteration. Experimentation suggests it is quite feasible to iterate, say, 5000 values but plot only the last 500.

Some of the images used in this article were created by modification of the above code.

Other Python code is in the further reading, for the logistic map.

Further reading:

Coding a bifurcation diagram, stackoverflow.com

Plotting in Python and chaotic-maps library

chaotic-maps seems to be a pure-Python library installable via pip from PyPI. It can plot orbit diagrams but also cobweb plots. It comes with a set of predefined maps, seen in map_library.py.

Further reading:

chaotic-maps, pypi.org
chaotic_maps, github.com/eowagner
map_library.py, github.com/eowagner

Plotting in PyPy

PyPy is a just-in-time compiler implementation of Python, standing in contrast to CPython. As long as one does only floating-point calculations and no plotting, it can be many times faster than CPython. A problem with PyPy is successfully installing matplotlib for it (numpy can easily be done without). It seems one can get matplotlib installed for PyPy, but it seems to require having a supporting C compiler installed, with required libraries.

To address the matplotlib problem, one can output the items to be plotted to a text file and plot it using e.g. Gnuplot.

Alternatively, one can have a single Python file with two parts, one for PyPy and one for CPython. The file would be called two times, each time with different parameters. The PyPy call would take care of the initial iterates that are not be plotted, outputting, say, the 100,000th iterate to a file. The CPython call would start the iteration where PyPy finished, reading the 100,000th iterates from the file. Similar idea is explored below in the Plotting in C language section.

The above idea is implemented below. The following script is to be called two times, first with "pypy" argument and via pypy, then with "py" argument and from normal Python. The speed difference between py and pypy concerning the first pass is huge. The second pass only works in normal Python unless one succeeds in installing matplotlib for pypy. The setup of parameters/presets is shared between the two passes, which would not be the case in the approach in which the first pass is coded in the C language and the second in Python.

Since we are now using high number of iterations--100,000 or even more--the question raised in section Calculation precision becomes more acute: should we use a high-precision floating-point library rather than Float64? Are we plotting largely noise created by chaotic-system amplification of numerical errors?

from __future__ import print_function

# To be called twice, first with "pypy", then "py".
# This version of presets skips early iterations (has high mapIterationCount while having much lower lastMapIterationCountToPlot),
# since that is the whole purpose of using pypy rather than py. Using pypy while plotting all iterates is unlikely to speed up things much.

import sys

pass1 = sys.argv[1] # "pypy" or "py"

preset = 1

def iteration(initX, r, iterationCount, xList=None):
  x = initX
  if xList is not None:
    xList[0] = x
  for i in xrange(iterationCount - 1):
    x = r * x * (1 - x)
    if xList is not None:
      xList[i + 1] = x
  return x
  
saveFigure = True
markerSize = 0.001
markerStyle = "," # . - point; , - pixel
progressReportBlockLen = 1000
interfaceTxtFileName = "_orbitDiagTmp.txt"
if preset == 1: # Complete diagram on the conventional range
  rInterval = (-0.1, 4.1)
  rResolution = 0.0005
  mapIterationCount = 100000
  lastMapIterationCountToPlot = 1000
  startingIterationValue = 0.5
  yLim = (-0.1, 1.1)
  fileName = "LogisticMapOrbitDiag_PyPy.png"
if preset == 2: # Complete diagram, including the negative r
  rInterval = (-2.1, 4.1)
  rResolution = 0.0005
  mapIterationCount = 100000
  lastMapIterationCountToPlot = 1000
  startingIterationValue = 0.5
  yLim = (-0.6, 1.6)
  fileName = "LogisticMapOrbitDiag_InclNegR_PyPy.png"
if preset == 3: # Zoom on the leftmost bifurcation region
  rInterval = (3.54, 3.5705) #3.571
  rResolution = 0.000004
  mapIterationCount = 100000
  lastMapIterationCountToPlot = 500
  startingIterationValue = 0.5
  yLim = (0.33, 0.9)
  fileName = "LogisticMapOrbitDiag_BifurcZoom_PyPy.png"
if preset == 4: # Complete diagram on the conventional range, skip early iterates
  rInterval = (-0.1, 4.1)
  rResolution = 0.0005
  mapIterationCount = 100000
  lastMapIterationCountToPlot = 1000
  startingIterationValue = 0.5
  yLim = (-0.1, 1.1)
  fileName = "LogisticMapOrbitDiag_SkipEarlyIter_PyPy.png"
  progressReportBlockLen = 100
if preset == 5: # Zoom into the period-3 window
  rInterval = (3.825, 3.865)
  rResolution = 0.000005
  mapIterationCount = 100000
  lastMapIterationCountToPlot = 1000
  startingIterationValue = 0.5
  yLim = (0.1, 1.0)
  fileName = "LogisticMapOrbitDiag_P3Zoom_PyPy.png"

if pass1 != "py" and pass1 != "pypy":
  print("The only supported passes are pypy (to be used first) and py (to be used second).")
elif pass1 == "pypy":
  outf = open(interfaceTxtFileName, "w")
  sys.stdout = outf
  lims = [0.0] * mapIterationCount
  itemCount = 0
  idx = 0
  r = rInterval[0]
  while r <= rInterval[1]:
    r = rInterval[0] + idx * rResolution
    idx += 1
    itemCount += 1
    x = iteration(startingIterationValue, r, mapIterationCount - lastMapIterationCountToPlot, xList = None)
    print(r, x)
    if itemCount % progressReportBlockLen == 0:
      print("Parameter values processed:", itemCount, file=sys.stderr)
elif pass1 == "py":
  import matplotlib.pyplot as plt
  fig, biax = plt.subplots()
  fig.set_size_inches(16, 9)

  x = [0.0] * lastMapIterationCountToPlot
  itemCount = 0
  for line in open(interfaceTxtFileName):
    itemCount += 1
    rStr, xStr = line.split(" ")
    r, xVal = float(rStr), float(xStr)
    iteration(xVal, r, lastMapIterationCountToPlot, xList = x)#__
    biax.plot([r] * lastMapIterationCountToPlot,
              x, "b.",
              markersize=markerSize,
              marker=markerStyle)
    if itemCount % progressReportBlockLen == 0:
      print("Parameter values processed:", itemCount, file=sys.stderr)
            
  biax.set(xlabel="c", ylabel="x",
           title="Logistic map - r * x (1 - x) - bifurcation diagram/orbit diagram" + "\n"
                 "r resolution: " + str(rResolution) +
                 ", iteration count: " + str(mapIterationCount) +
                 ", last values to plot count: " + str(lastMapIterationCountToPlot))

  biax.set_ylim(yLim[0], yLim[1])
  if saveFigure:
    plt.savefig(fileName, dpi=600)
  else:
    plt.show()

Plotting in Matlab

The following code for Matlab for the Mandelbrot set quadratic map was in an old version of the Wikipedia article Bifurcation diagram:

c=0; 
y=0.0;  

hold on 
while c < 2
        for i=1:100; 
            y = y.^2 - c; %converge the iteration
            plot(c,y,'.'); % plot the converged points
        end 
        c=c+0.01;
end

Unlike the Mandelbrot set quadratic map, there is "-c" rather than "+c" above; thus, the diagram will look flipped compared to the Mandelbrot set quadratic map.

Plotting in Wolfram Alpha

Entering "bifurcation diagram of logistic map" into Wolfram Alpha plots the diagram. Wolfram Alpha plots not only the orbit diagram, but also an orbit time series, the cobweb diagram, etc.

Further reading:

bifurcation diagram of logistic map, wolframalpha.com

Plotting in Desmos

One can plot the orbit diagram/bifurcation diagram in Desmos, an online graphing calculator. A link is in the further reading.

Further reading:

Logistic Map 2, desmos.com

Plotting in the C language

One could plot the orbit diagram in the C language for maximum performance (unless one wants to go for assembly). The key idea is that calculating the initial iterates to be skipped (not plotted) requires no special facilities, unlike plotting the final iterates.

One could output the items to be plotted to a text file and then feed the text file to a plotting program such as Gnuplot or to a custom Python program that would use matplotlib. This would be useful for the manner of plotting in which one wishes to skip a large number of initial iterates and plot only a relatively small number of last iterates (it would keep the text file reasonably small). Thus, one could run e.g. 100,000 iterations in C without plotting and only output the subsequent 1000 values for each parameter c.

Alternatively, if one wants to have a high number of skipped iterates and low number of plotted iterates, one could use C to output the last non-plotted iterate to a text file and continue the iteration from Python. This would minimize the size of the text file used as an interface.

Alternatively, one could find a plotting library for C and directly call it from C.

Alternatively, one could create a Python module in C to be called from a Python script.

Alternatively, one could plot the results using a raster image file format that is easy to create in C, to be converted to PNG using an external tool.

Calculation precision

Since we are dealing with systems described in sources as showing chaos (a small difference in an initial condition or an intermediate state of the system gets greatly amplified), one may ask whether the Float32 or Float64 (default in Python) is sufficient for good accuracy. That is to say, one may wonder whether some plotted pixels are a result of accumulated numerical error, possibly very much so in the region very close to the boundary of the chaotic region (close to 4 for the logistic map, and close to -2 for the Mandelbrot set quadratic map). In Python, one may use the decimal module to obtain a higher precision. One would need to do more research to find out whether this problem is merely hypothetical.

History

An orbit diagram looking like one for the logistic map appears in the 1985 book Metamagical Themas by Douglas Hofstadter, where it is credited to Road to Chaos by Leo P. Kadanoff in Physics Today, 1983, and a reference is made to J. P. Crutchfield, D. Farmer and B. A. Huberman, Physics Reports, Vol. 92, pp. 45-82, December 1982. The history of the diagram can well reach farther. According to Nusse et al. 1994, bifurcation diagrams appeared in the seventies.^[5]

Term orbit diagram

This section is something of an appendix, to substantiate the use of the term "orbit diagram". The attesting quotations are as follows:

"The orbit diagram shows the long-term behaviour for all values of r at once: ..."^[6]
"Here is the orbit diagram for the quadratic family q_c(x) =x² + c for −2 ≤ c ≤0:"^[7]
"The algorithm for drawing the orbit diagram is described in many places, so we'll make it quick here, with the numbers actually used: at each of 4000 equally spaced values of c from left to right in [-2.25, 0.375], the first 2000 members of the orbit of zero for f_c(x) = x² + c were computed, and then the next 200 members actually plotted."^[8]
"The Orbit Diagram and the Mandelbrot Set"^[9]

Term bifurcation diagram

Chip Ross uses the term "bifurcation diagram" to refer to a somewhat different if related type of diagram.^[8] In that diagram, e.g. point [c=0, x=1] is plotted, which is a fixed point of x**2 + c for c=0 (since 1**2 --> 1) but is a repelling one, not attracting one, and the system starting from zero never visits that fixed point.

A page at vanderbilt.edu uses the term "bifurcation diagram" for what we here refer to as "orbit diagram"^[10]

Nusse et al. 1994 states: "The concept of bifurcation diagram includes a number of ways of plotting a phase variable on one axis and a parameter on another."^[5] This suggests a broader concept than just orbit diagram. (The source has the parameter of the map on the y-axis, while most diagrams have it on the x-axis.)

Term Feigenbaum diagram

Term "Feigenbaum diagram" is sometimes used to refer to this kind of diagram. It would need to be clarified whether the term only refers to an orbit diagram for the logistic map or whether it refers to any other discrete-time map orbit diagram as well.

Attesting quotations:

"The Feigenbaum diagram for the logistic map shows the periodic doubling regime"^[11]
"The bifurcation diagram (also known as the orbit diagram, the Feigenbaum diagram, or the road to chaos) is a plot of later iterates, showing the eventual state of the system."^[12]

References

↑ Logistic Map, mathworld.wolfram.com
↑ Quadratic Map, mathworld.wolfram.com
↑ Tent Map, mathworld.wolfram.com
↑ Ricker’s Population Model by Willie Bell, James Boffenmyer, Scott Dean, Andrew Stewart, 2010, math.lsu.edu
↑ ^5.0 ^5.1 Bifurcation Diagrams by Helena E. Nusse, James A. Yorke & Eric J. Kostelich , 1994
↑ Chaos from maps, physics.ox.ac.uk
↑ Dynamics, Chaos, and Fractals (part 2): Dynamics of One-Parameter Families by Evan Dummit, 2015, northeastern.edu
↑ ^8.0 ^8.1 Of Bifurcation and Orbit Diagrams by Chip Ross, abacus.bates.edu
↑ The Orbit Diagram and the Mandelbrot Set by Robert L. Devaney, 1991, tandfonline.com
↑ Bifurcation Diagram, vanderbilt.edu
↑ Studying the Cantor Dust at the Edge of Feigenbaum Diagrams, May 1998, maa.org
↑ The Road to Chaos is Filled with Polynomial Curves by Richard D. Neidinger and R. John Annen, The American Mathematical Monthly, 1996