Cobwebbing

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

A discrete one-dimensional recursive function is described by a mapping of the form

(Subscript[x, i + 1])  =  f[μ, (Subscript[x, i])]

where μ is a parameter that governs the dynamical properties of the system.  A typical example is the logistic equation, (Subscript[x, i + 1])  =  μ (Subscript[x, i])(
            1 - (Subscript[x, i])), often used to describe population dynamics for species with nonoverlapping generations.  See logistic map.nb for a more complete description of that model.

Points which are mapped upon themselves, such that (Subscript[x, i + 1]) = (Subscript[x, i]), are called  fixed points of the map and play a central role in the analysis of the behavior of iterative maps.  As the map is iterated, nearby points may remain at fixed distance, be driven toward, or be driven away from a fixed point.  Hence, fixed points may be classified as either neutral, stable, or unstable.

In our earlier numerical study of the logistic equation we discovered that the system can exhibit several qualitatively different types of behavior depending upon the value of the growth parameter, μ.  For small μ we found that there was a single stable fixed point at x = 0, while for somewhat larger μ the origin became unstable and a new stable fixed point emerged for finite x.  This fixed point remained stable until μ crossed a bifurcation value at which two new solutions emerged and the system settled upon a stable oscillation between two values of x, a 2-cycle.  As μ increased further, these solutions bifurcated again as a stable 4-cycle emerged.  As μ continued to increase the interval between successive period doublings decreased until a chaotic regime was reached.  For slightly larger μ a brief island of order interrupts the chaotic regime.

In this notebook we analyze these behaviors more rigorously, using graphical and analytical techniques to explain these behaviors and to calculate the values of the fixed points and the bifurcation parameters.  One of the basic tools is the cobweb plot, which illustrates the iteration sequence for a recursive function.  Using the cobweb plot, we will discover the conditions which govern the stability of a fixed point and will introduce the methodology of linear stability analysis.  Although we concentrate here upon the logistic equation, these techniques are much more general.

Much of the present discussion is based upon Chapter 2 of the book

Mathematical Biology, by J.D. Murray (Springer-Verlag, 1989)

but by using Mathematica we can develop versatile tools and empower the student to perform his/her own investigations.

Initialization

ClearAll["Global`*"]
Off[General::spell,  General::spell1,  Plot::plnr];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Cobweb plot

A manual calculation for a discrete one-dimensional recursive function proceeds by the following iteration process.

1) select (Subscript[x, 1]), i = 1, (Subscript[y, 1]) = (Subscript[x, 1]).

2) (Subscript[y, i + 1]) ↤ f[(Subscript[x, i])]

3) i ↤ i + 1

4) (Subscript[x, i]) ↤ (Subscript[y, i])

5) return to 2)

A cobweb plot tracks the sequence of functional evaluations using line segments to represent each assignment of either x or y during this iteration.  Hence, step 2) is represented by a vertical the line segment between the diagonal y ↤ x and the curve y = f[x] while step 4) is represented by a horizontal line segment from that curve back to the diagonal x ↤ y.

A function which produces a cobweb plot is defined below.  NestList is used to produce a sequence of iterates.  In order to focus on the long-term behavior, one can drop the first ndrop iterates to allow transients to decay.  When the entire sequence is displayed, the starting point is indicated by a prominent dot.  Line segments between the line (Subscript[x, i + 1]) = (Subscript[x, i]) and the continuous curve y = f[x] are constructed using Map and Line and drawn with Graphics.  Finally, Show is used to combine the line segments with plots of the diagonal line and the continuous function.

\(\(\(cobweb::usage\  = \ "cobweb[f,x0,nmax,ndrop] produces a cobweb plot for the recursive function f[x] with initial value x0, suppressing the first ndrop iterates.  Options may also be passed to Plot.";\)\(
\)
    \)		\)
\(cobweb[f_, x0_, nmax_, \ ndrop_:  0, \ opts___Rule] := 
      Module[{plot1, plot2, plot3}, 
plot1 = 
          Plot[{x, f[x]}, {x, 0, 1}, \ opts, Frame -> True, \ 
            FrameLabel -> {\*""(Subscript[x, i])"", \ \*""(Subscript[x, i + 1])""}, 
            PlotLabel -> 
              StringJoin["Cobweb Plot for: ", \ ToString[f[x]]\ ], 
            DisplayFunction -> Identity]; 
plot2 = 
          Graphics[\ 
            Map[\ \ Line[{{#, #}, {#, f[#]}, {f[#], f[#]}}]\  &, \ 
              Drop[NestList[f, x0, nmax], ndrop]], 
            DisplayFunction -> Identity]; 
plot3 = 
          ListPlot[{{x0, x0}}, PlotStyle -> PointSize[0.03], 
            DisplayFunction -> Identity]; 
If[
          ndrop == 0, 
          Return[Show[{plot1, plot2, plot3}, 
              DisplayFunction -> $DisplayFunction]], 
          Return[Show[{plot1, plot2}, 
              DisplayFunction -> $DisplayFunction]]]]\)

Example: logistic map

Definition of the mapping function

The logistic equation is represented by the pure function μ # (1-#)&.

logisticCobweb[μ_, x0_, nmax_, ndrop_:  0, opts___Rule]  := 
      cobweb[μ  # (1 - #) &, x0, nmax, ndrop, opts]

Extinction for μ<1

logisticCobweb[0.8, 0.5, 20, 0];

logisticCobweb[0.8, 0.8, 20, 0];

Approach to stable population for 1<μ<3

logisticCobweb[1.5, 0.1, 20, 0];

logisticCobweb[2.5, 0.7, 20, 0, 
        PlotRange -> {{0, 1}, {0, 1}}];

Emergence of stable cycles for 3<μ≲3.8

logisticCobweb[3.2, 0.2, 20, 0];

Clearly, with μ ≈ 3.2 this system exhibits a stable 2-cycle after the initial transients have decayed.

logisticCobweb[3.5, 0.2, 100, 50];

Thus, after the transients have died away, we find a stable 4-cycle for μ ≈ 3.5.

Chaos

logisticCobweb[3.8, 0.2, 200, 100];

In the chaotic regime the cobweb plot appears to visit practically every point within a large but finite range eventually.  Nevertheless, for a slightly larger value of μ order miraculously re-emerges!  However, the appearance of a 3-cycle is a new feature — for smaller μ we had observed (2^n) cycles arising from bifurcations.

logisticCobweb[3.84, 0.2, 200, 100];

Fixed points and cycles

Any point which is mapped upon itself, such that (Subscript[x, i + 1]) = (Subscript[x, i]), is called a fixed point of the map.  Similarly, points which repeat after n iterations, such that (Subscript[x, i + n]) = (Subscript[x, i]), are members of an n-cycle.  The dependence of these points upon the map parameter μ can be determined by examining the iterated function sequence (Subscript[f, 0])[x] = x, (Subscript[f, 1])[x] = f[x], (Subscript[f, 2])[x] = f[f[x]], ⋯, (Subscript[f, n])[x] = f[(Subscript[f, n - 1])[x]].

Clear[f];
(Subscript[f, n])_[μ_, x_] := Nest[μ # (1 - #) &, x, n]

Table[(Subscript[f, n])[μ, x], {n, 0, 4}] // Simplify // 
      TableForm

Thus, fixed points are located by solving the equation (Subscript[f, 1])[x] == x and n-cycles are constructed by solving (Subscript[f, n])[x] == x.

fixedPoints[n_] := Solve[(Subscript[f, n])[μ, x] == x, x]

It is instructive to display a graphical solution to this equation for a range of μ.

\(\(Plot[
        Evaluate[{x, Table[f\_1[μ, x], {μ, 0.5, 4, 0.5}]}], {x, 0, 
          1}, Frame -> True, \ 
        FrameLabel -> {\*""(Subscript[x, i])"", \ \*""(Subscript[x, i + 1])""}, AspectRatio -> Automatic];\)\)

For μ ≤ 1, the only solution is x = 0 corresponding to the situation in which a population that fails to replace itself becomes extinct.  For larger μ a new fixed point emerges.

sol1 = fixedPoints[1]

Linear stability analysis

A fixed point is described as stable if points in its neighborhood are attracted toward the fixed point or as unstable if small displacements from the fixed point are driven away as the map is iterated.  The property of the mapping function which distinguishes between these behaviors is the slope, λ, of the map function at the fixed point.  By visualizing the cobweb process, it becomes apparent that |λ| < 
        1 yields a stable approach toward the fixed point while |λ| > 
        1 results in repulsion from the fixed point.  More formally, suppose that (Subscript[x, 0]) is a fixed point of the mapping (Subscript[x, i + 1]) = f[μ, (Subscript[x, i])], such that (Subscript[x, 0]) = f[μ, (Subscript[x, 0])].  Consider a nearby point (Subscript[x, 0]) + δ x for which f[μ, 
          (Subscript[x, 0]) + δ x] ≈ 
        (Subscript[x, 0]) + (f^′)[μ, (Subscript[x, 0])] δ x where λ = (f^′)[μ, (Subscript[x, 0])] is the slope of f[μ, x] at the fixed point.  Hence, the iteration produces \(TraditionalForm\`x\_\(i + 1\) ≈ 
            x\_0 + λ\ δ x\_i\  ⇒ \ δ x\_\(i + 1\) = λ\ δ x\_i\) where δ (Subscript[x, i]) = (Subscript[x, i]) - (Subscript[x, 0])is the distance from the fixed point.  Thus, near a fixed point each iteration multiplies the distance from the fixed point by a factor of λ.  The slope at the fixed point is known as its  eigenvalue by analogy with discrete multidimensional mappings effected by matrix multiplication.  Clearly, if |λ| < 
        1 nearby points are driven toward the fixed point as the mapping is iterated,  if |λ| > 
        1 nearby points are driven away from the fixed point, while if |λ| = 
        1 nearby points orbit the fixed point at constant distance.  Therefore, the nature of a fixed point is determined by its eigenvalue and is classified as stable if the magnitude of the eigenvalue is less than unity, unstable if greater than unity, or neutrally stable if exactly equal to unity.  Furthermore, if the eigenvalue is positive the mapping approaches or recedes monotonically, while for negative eigenvalues it oscillates (orbits) around the fixed point.

This type of analysis is known as linear stability analysis and can be applied not only to arbitrary discrete mappings but to many continuous models of dynamical systems.  Many of the features characteristic of discrete maps, such as period doubling, chaos, order within chaos, and self-similarity are found in other dynamical systems governed by nonlinear differential equations also.

Thus, we define a function which evaluates the slope of the mapping function at the fixed points

(Subscript[λ, n])_[μ_] := 
      D[(Subscript[f, n])[μ, x], x] /. fixedPoints[n] // Simplify

and for the first iteration of the logistic map we find

(Subscript[λ, 1])[μ]

Therefore, stability for the first fixed point with x > 0 requires 1 < μ < 3, which is consistent with our earlier numerical investigations.  To illustrate this behavior, consider the following cobweb plots near the critical value μ=3.  We use PlotRange to zoom in upon the vicinity of a fixed point and use AspectRatio->Automatic to obtain a faithful picture of the slope of the mapping function at this point.

For μ just below the critical point we observe that nearby points are attracted toward the fixed point.

x /. 
        sol1[[2]] /. μ -> 
        2.8

logisticCobweb[2.8, 0.6, 10, 0, 
        PlotRange -> {{0.5, 0.7}, {0.5, 0.7}}, 
        AspectRatio -> Automatic];

Conversely, for μ just above the critical point we observe that nearby points are repelled by the fixed point.

x /. 
        sol1[[2]] /. μ -> 
        3.2

logisticCobweb[3.2, 0.7, 10, 0, 
        PlotRange -> {{0.6, 0.8}, {0.6, 0.8}}, 
        AspectRatio -> Automatic];

Pitchfork bifurcation mechanism

Our numerical investigations also suggested that for μ not too far above 3 the system settles upon a stable oscillation between two values, a 2-cycle.  These two values must be found among the solutions to the equation (Subscript[f, 2])[μ, x] == x,  which has the same structure as the original mapping except that the step size is 2.  Once again it is instructive to examine this equation graphically.

\(\(Plot[
        Evaluate[{x, Table[f\_2[μ, x], {μ, 2.8, 3.6, 0.2}]}], {x, 0, 
          1}, Frame -> True, 
        FrameLabel -> {\*""(Subscript[x, i])"", \*""(Subscript[x, i + 2])""}, AspectRatio -> Automatic];\)\)

We can apply the same analytical techniques to a cobweb plot for a two-step map as we did for a one-step map.  Thus, we identify the intersections of the quartic curve with the line as the fixed points of a 2-cycle and recognize that stability of such points can be classified according to their slopes, where again |λ| < 
        1 corresponds to stable and |λ| > 
        1 to unstable fixed points.

sol2 = fixedPoints[2] // Simplify

One observes that the fixed points of the first iterate remain fixed points of the second iterate — if a point is mapped into itself by one iteration it will surely be mapped into itself for two iterations also!  The point at the origin is the only stable fixed point for μ < 1, representing extinction for small growth parameters.  The second fixed point is stable only for growth parameters in the range 1 < μ < 3.  The third and fourth solutions of (Subscript[f, 2]) == x are complex for μ < 3, degenerate with the second solution for μ = 3, and real for μ > 3.  Examining the slopes,

slopes[2] = (Subscript[λ, 2])[μ] // Simplify

slopes[2] /. μ -> 3

we find that the three degenerate fixed points are neutrally stable, λ = 1, at the critical point μ = 3.  For larger μ, the previously stable parent becomes unstable while two stable offspring emerge and diverge from it as μ increases.  Thus, the previously stable fixed passes through a pitchfork bifurcation as it gives birth to a stable 2-cycle.

Show[{Plot[Evaluate[x /. sol1], {μ, 1, 4}, 
            DisplayFunction -> Identity], 
          Plot[Evaluate[x /. sol2], {μ, 3, 4}, 
            DisplayFunction -> Identity]}, 
        PlotRange -> {{2.5, 4.0}, {0.2, 1.0}}, 
        PlotLabel -> "First Pitchfork Bifurcation", 
        Frame -> True,  
        FrameLabel -> {"μ", "fixed point"}, 
        DisplayFunction -> $DisplayFunction];

However, their slopes decrease with μ until a new critical point is reached where the slopes cross through λ < -1.

Solve[slopes[
              2][[3]] == -1, μ][[2]]

% // N

At this point, each fixed point gives birth to two stable offspring as it crosses the threshold of instability.  Thus, this second period doubling produces a stable 4-cycle.  However, numerical evaluation of the fixed points is needed for cycles of order 3 or higher.  Nevertheless, the solutions can be displayed easily.

sol4  =  fixedPoints[4] // Simplify;

Plot[Evaluate[{x /. sol2, N[x /. sol4]}], {μ, 2.8, 3.84}, 
        PlotLabel -> "Second Pitchfork Bifurcation", 
        Frame -> True,  
        FrameLabel -> {"μ", "fixed point"}, 
        DisplayFunction -> $DisplayFunction];

The interval between successive bifurcations decreases as we approach the chaotic regime.

Tangent bifurcation and the onset of chaos

Sarkovskii proved a theorem for one-dimensional discrete maps that if there exists a stable 3-cycle for \(TraditionalForm\`μ\_c ≤ μ ≤ μ\_p\), then aperiodic or chaotic solutions exist in this range.  The bifurcation pattern appears to be random within the chaotic regime.  Also note that if (Subscript[μ, p]) is finite, order re-emerges for μ > (Subscript[μ, p]), but is usually confined to a narrow range before the another transition to chaos is encountered.  The order to chaos to order pattern then repeats with increasing frequency.

The condition that must be satisfied for the emergence of a stable 3-cycle is the existence of three fixed points in (Subscript[f, 3])[(Subscript[μ, c]), x] == x with λ = 1, such that (Subscript[f, 3]) is tangent to the diagonal at three places.  As illustrated below, this condition is satisfied for the logistic map when (Subscript[μ, c]) ≈ 3.828.

\(\(Plot[
        Evaluate[{x, Table[f\_3[μ, x], {μ, 3.728, 3.928, 0.1}]}], {x, 
          0, 1}, Frame -> True, \ 
        FrameLabel -> {\*""(Subscript[x, i])"", \ \*""(Subscript[x, i + 3])""}, PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, 
        PlotPoints -> 50];\)\)

For larger μ each of the tangent points bifurcates, but unlike the pitchfork bifurcation, fixed points that emerge via tangent bifurcation are not preceded by an earlier fixed point that loses its stability at the bifurcation point.  Furthermore, unlike λ = -1 bifurcations, bifurcations at λ = +1 typically produce unstable, aperiodic solutions.

sol3  =  fixedPoints[3] // Simplify;

Plot[Evaluate[N[x /. sol3]], {μ, 3.8, 4.0}, 
        PlotRange -> {Automatic, {0, 1}}, 
        PlotLabel -> "Tangent Bifurcation", Frame -> True,  
        FrameLabel -> {"μ", "fixed point"}, 
        DisplayFunction -> $DisplayFunction];