Consider a one-parameter family of discrete-time dynamical systems (maps) \( x_{n+1} = f(x_n;\mu) \), where \( \mu \) is a real parameter. A bifurcation occurs when a small change in \( \mu \) causes a qualitative change in the long-term dynamics (e.g., a stable fixed point becomes unstable; a stable periodic orbit appears; chaos emerges).

A bifurcation diagram is a parametric portrait: for each \( \mu \), we plot asymptotic samples of the state (or some scalar observable). In practice, we iterate (or simulate), discard a transient, then record a set of iterates (or Poincaré samples) to approximate the invariant attracting set.

flowchart TD
  A["Choose parameter range (mu_min..mu_max)"] --> B["For each mu: initialize state x0"]
  B --> C["Evolve dynamics for N_transient steps (discard)"]
  C --> D["Collect N_keep samples (iterates or section hits)"]
  D --> E["Plot samples versus mu (scatter)"]
  E --> F["Interpret: fixed points, periods, bands, crises"]
  F --> G["Validate numerics: step size, transients, precision"]
        

In continuous time, \( \dot{\mathbf{x}} = \mathbf{F}(\mathbf{x};\mu) \), we cannot directly plot \( \mathbf{x}(t) \) for all \( t \) in a bifurcation diagram. We instead sample using: (i) maxima/minima of a coordinate, (ii) stroboscopic sampling \( t = kT \), or (iii) a Poincaré section producing a discrete return map (discussed in Section 4).

A fixed point \( x^\star(\mu) \) satisfies \( x^\star = f(x^\star;\mu) \). Define the perturbation \( e_n = x_n - x^\star \). A first-order Taylor expansion gives

\[ e_{n+1} = f(x^\star + e_n;\mu) - f(x^\star;\mu) = f_x(x^\star;\mu)\,e_n + \mathcal{O}(e_n^2), \]

where \( f_x \) denotes derivative with respect to \( x \). Hence for sufficiently small perturbations, the dynamics is approximately linear: \( e_{n+1} \approx a e_n \) with \( a = f_x(x^\star;\mu) \).

Local stability criterion (1D maps):

\[ \left| f_x(x^\star;\mu) \right| \; < \; 1 \;\;\Rightarrow\;\; x^\star \text{ is locally asymptotically stable}, \\ \left| f_x(x^\star;\mu) \right| \; > \; 1 \;\;\Rightarrow\;\; x^\star \text{ is unstable}. \]

Proof sketch (contraction argument): If \( |a| < 1 \), then the linear recursion gives \( |e_n| \approx |a|^n |e_0| \to 0 \). For the nonlinear system, pick a neighborhood where the remainder satisfies \( |\mathcal{O}(e_n^2)| \le c |e_n|^2 \); then for sufficiently small \( e_0 \) one can show \( |e_n| \) decreases geometrically, because the quadratic term becomes dominated by the linear contraction. (This is a standard local stability theorem.)

Bifurcation as loss of stability: A common route is when the multiplier \( f_x(x^\star;\mu) \) crosses the unit circle boundary. In 1D maps:

• \( f_x(x^\star;\mu_c) = +1 \) (typically saddle-node / fold or pitchfork-like behavior, depending on symmetry).
• \( f_x(x^\star;\mu_c) = -1 \) (flip / period-doubling bifurcation).

The flip case is the backbone of the period-doubling cascade used in bifurcation diagrams of many nonlinear systems.

The logistic map (introduced in Lesson 1 as a canonical nonlinear map) is \( x_{n+1} = f(x_n;r) = r x_n(1-x_n) \), typically studied on \( 0 \le x_n \le 1 \) with \( 0 \le r \le 4 \).

3.1 Fixed points and their stability

Fixed points satisfy \( x = r x(1-x) \), yielding:

\[ x_1^\star = 0,\qquad x_2^\star = 1 - \frac{1}{r}\quad (r \neq 0). \]

The derivative is \( f_x(x;r) = r(1-2x) \). Evaluate at the fixed points:

\[ f_x(x_1^\star;r) = r,\qquad f_x(x_2^\star;r) = r\left(1-2\left(1-\frac{1}{r}\right)\right)=2-r. \]

Stability conditions \( |f_x| < 1 \) give:

\[ |r| < 1 \Rightarrow x_1^\star \text{ stable (for the usual } r \ge 0 \text{, this is } 0 \le r < 1), \]

\[ |2-r| < 1 \;\;\Leftrightarrow\;\; 1 < r < 3 \Rightarrow x_2^\star \text{ stable}. \]

At \( r=3 \), the multiplier at \( x_2^\star \) is \( f_x = -1 \), the hallmark of a flip (period-doubling) bifurcation.

3.2 Birth of the period-2 orbit (explicit formula)

A period-2 orbit satisfies \( f(f(x;r);r)=x \) but is not a fixed point. Compute:

\[ f(f(x;r);r) - x = 0. \]

The polynomial factorization separates fixed points from true 2-cycles; the 2-cycle points are the roots of:

\[ r^2 x^2 - r^2 x - r x + r + 1 = 0, \]

yielding the two points of the 2-cycle:

\[ x_{\pm}^{(2)}(r) = \frac{r+1 \pm \sqrt{(r-3)(r+1)}}{2r}. \]

Real 2-cycle points require \( (r-3)(r+1) \ge 0 \), so for the standard regime \( r \ge 0 \) this starts at \( r \ge 3 \), consistent with the flip bifurcation at \( r=3 \).

3.3 Stability of the 2-cycle and the next doubling threshold

For a period-2 orbit \( \{x_1, x_2\} \), the linearized multiplier over one full period is the product: \( m_2 = f_x(x_1;r)f_x(x_2;r) \). Stability requires \( |m_2| < 1 \). For the logistic map, this simplifies to a closed form:

\[ m_2(r) = -r^2 + 2r + 4. \]

Therefore the 2-cycle is stable when:

\[ | -r^2 + 2r + 4 | \; < \; 1 \;\;\Rightarrow\;\; 3 \; < \; r \; < \; 1 + \sqrt{6}. \]

The next flip bifurcation (period-2 to period-4) occurs at \( r_2 = 1+\sqrt{6} \approx 3.44949 \). Repeating this process yields a cascade of doublings: \( 2 \to 4 \to 8 \to \cdots \), accumulating at a finite parameter value \( r_\infty \).

flowchart TD
  R1["r in (1,3): stable period-1 (fixed point)"] --> R2["r=3: flip bifurcation"]
  R2 --> R3["r in (3, r2): stable period-2"]
  R3 --> R4["r=r2: flip (2->4)"]
  R4 --> R5["period-4 then 8 then 16 ..."]
  R5 --> R6["accumulation at r_infty: onset of chaotic bands"]
        

3.4 Feigenbaum scaling (universality, statement)

Let \( r_n \) denote the parameter where the stable period \( 2^{n-1} \) orbit loses stability and a stable period \( 2^n \) orbit is born. For a large class of smooth unimodal maps (including the logistic map), the ratio of successive interval lengths converges to the Feigenbaum constant \( \delta \):

\[ \delta = \lim_{n\to\infty}\frac{r_n - r_{n-1}}{r_{n+1} - r_n} \approx 4.6692016\ldots \]

The key point is universality: the same \( \delta \) appears for many different nonlinear systems with the same local quadratic maximum structure in their return map, explaining why period-doubling cascades are so commonly observed in engineering models.

For continuous-time systems \( \dot{\mathbf{x}}=\mathbf{F}(\mathbf{x};\mu) \), an attractor is typically defined as a compact invariant set \( \mathcal{A} \subset \mathbb{R}^d \) with a basin of attraction such that trajectories starting in the basin approach \( \mathcal{A} \). Invariance means the flow map \( \varphi_t \) satisfies:

\[ \varphi_t(\mathcal{A}) = \mathcal{A}\quad \text{for all } t \ge 0. \]

An attractor is called strange when it combines attracting behavior with complicated geometry (often fractal-like) and strong orbit separation (sensitive dependence). In practice, a standard diagnostic is that nearby initial conditions diverge exponentially for some time, indicating at least one positive Lyapunov exponent (formally studied in Lesson 3).

Poincaré section and return map: Choose a hypersurface \( \Sigma \) transverse to the flow, and define the return map \( P:\Sigma \to \Sigma \) by mapping a point to its next intersection with \( \Sigma \). Then continuous-time dynamics reduces to a discrete map, and bifurcation diagrams can be built by iterating \( P \).

The Lorenz system is the canonical example of a continuous-time strange attractor:

\[ \dot{x} = \sigma(y-x),\qquad \dot{y} = x(\rho - z) - y,\qquad \dot{z} = xy - \beta z. \]

For classical parameters (e.g., \( \sigma=10, \rho=28, \beta=8/3 \)), trajectories approach a bounded set with repeated stretching and folding. Computationally, we visualize it by integrating the ODE and plotting \( (x(t),y(t),z(t)) \) after discarding transients.

How to create a bifurcation diagram for an ODE: pick a scalar observable \( s(t)=h(\mathbf{x}(t)) \) and record a discrete sample set, e.g. maxima of \( x(t) \) or Poincaré hits. Then plot these sample values versus a parameter (e.g., \( \rho \) in Lorenz).

Bifurcation diagrams are conceptually simple but numerically subtle. The standard pipeline is:

• Parameter sweep: choose \( \mu_i \) grid, often uniform, sometimes adaptive near transitions.
• Transient removal: discard the first \( N_{\text{trans}} \) iterates (or time units) so that the orbit approaches the attractor.
• Sampling: record \( N_{\text{keep}} \) asymptotic points (map) or event-based points (ODE).
• Plot: scatter plot \( (\mu_i, x_n) \) (or \( (\mu_i, h(\mathbf{x}(t_k))) \)).
• Plot: scatter plot \( (\mu_i, x_n) \) (or \( (\mu_i, h(\mathbf{x}(t_k))) \)).

Key pitfalls:

• Insufficient transient: may plot non-asymptotic behavior and fake branches.
• Step size / solver tolerance (ODE): coarse integration can destroy invariant geometry; always verify by refining step size/tolerances.
• Finite precision: chaotic systems amplify rounding; this is not a bug but must be managed (repeatability, seeding, reporting tolerances).
• Multistability: multiple coexisting attractors may exist; bifurcation diagram can depend on initial condition.

In the labs below, the logistic bifurcation diagram uses a vectorized sweep (Python) and simple CSV output (C++/Java). The Lorenz attractor is visualized by numerical integration (RK4 or ode45).

File: Chapter20_Lesson2.py
Libraries used: numpy, matplotlib. This script: (i) generates a logistic map bifurcation diagram (PNG + CSV) and (ii) simulates the Lorenz system via RK4 and saves a 3D plot (PNG + CSV).

# Chapter20_Lesson2.py
# System Dynamics — Chapter 20 (Chaos, Complex Dynamics, and Computational Tools)
# Lesson 2: Bifurcation Diagrams, Period Doubling, and Strange Attractors
#
# Outputs:
#   - logistic_bifurcation.png (bifurcation diagram)
#   - logistic_bifurcation.csv (points used in the diagram)
#   - lorenz_attractor.png (3D trajectory projection)
#   - lorenz_traj.csv (Lorenz trajectory samples)

import numpy as np
import matplotlib.pyplot as plt

def logistic_map(x, r):
    return r * x * (1.0 - x)

def bifurcation_logistic(r_min=2.5, r_max=4.0, n_r=6000,
                         n_transient=1200, n_keep=200, x0=0.123456):
    rs = np.linspace(r_min, r_max, n_r)
    x = np.full(n_r, x0, dtype=float)

    # Burn-in / transient
    for _ in range(n_transient):
        x = logistic_map(x, rs)

    # Keep last n_keep iterates for each r
    r_out = np.repeat(rs, n_keep)
    x_out = np.empty(n_r * n_keep, dtype=float)
    idx = 0
    for _ in range(n_keep):
        x = logistic_map(x, rs)
        x_out[idx:idx+n_r] = x
        idx += n_r

    # Save CSV (r, x)
    data = np.column_stack([r_out, x_out])
    np.savetxt("logistic_bifurcation.csv", data, delimiter=",", header="r,x", comments="")

    # Plot
    plt.figure(figsize=(10, 5))
    plt.plot(r_out, x_out, ",", markersize=1)
    plt.xlabel("r")
    plt.ylabel("x (asymptotic samples)")
    plt.title("Logistic map bifurcation diagram")
    plt.tight_layout()
    plt.savefig("logistic_bifurcation.png", dpi=200)
    plt.close()

def rk4_step(f, t, y, h):
    k1 = f(t, y)
    k2 = f(t + 0.5*h, y + 0.5*h*k1)
    k3 = f(t + 0.5*h, y + 0.5*h*k2)
    k4 = f(t + h, y + h*k3)
    return y + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

def lorenz_rhs(sigma=10.0, rho=28.0, beta=8.0/3.0):
    def f(t, y):
        x, yv, z = y
        return np.array([
            sigma*(yv - x),
            x*(rho - z) - yv,
            x*yv - beta*z
        ], dtype=float)
    return f

def simulate_lorenz(T=40.0, h=0.005, y0=(1.0, 1.0, 1.0),
                    sigma=10.0, rho=28.0, beta=8.0/3.0,
                    transient=5.0, sample_stride=4):
    f = lorenz_rhs(sigma, rho, beta)
    n = int(T / h)
    y = np.array(y0, dtype=float)
    t = 0.0

    # Run and store after transient
    traj = []
    for i in range(n):
        y = rk4_step(f, t, y, h)
        t += h
        if t >= transient and (i % sample_stride == 0):
            traj.append((t, y[0], y[1], y[2]))
    traj = np.array(traj, dtype=float)
    np.savetxt("lorenz_traj.csv", traj, delimiter=",", header="t,x,y,z", comments="")

    # Plot 3D trajectory (requires mpl_toolkits)
    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
    fig = plt.figure(figsize=(7, 6))
    ax = fig.add_subplot(111, projection="3d")
    ax.plot(traj[:,1], traj[:,2], traj[:,3], linewidth=0.4)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.set_title("Lorenz attractor trajectory (sampled)")
    plt.tight_layout()
    plt.savefig("lorenz_attractor.png", dpi=200)
    plt.close()

if __name__ == "__main__":
    bifurcation_logistic()
    simulate_lorenz()
    print("Done. Generated PNGs and CSVs in the current directory.")
      

File: Chapter20_Lesson2.cpp
This program outputs point clouds as CSV: logistic_bifurcation_cpp.csv and lorenz_traj_cpp.csv. Plot them using your preferred tool (Python/MATLAB/gnuplot).

// Chapter20_Lesson2.cpp
// System Dynamics — Chapter 20 (Chaos, Complex Dynamics, and Computational Tools)
// Lesson 2: Bifurcation Diagrams, Period Doubling, and Strange Attractors
//
// Build (example):
//   g++ -O2 -std=c++17 Chapter20_Lesson2.cpp -o ch20_l2
// Run:
//   ./ch20_l2
//
// Outputs:
//   logistic_bifurcation_cpp.csv  (r,x points)
//   lorenz_traj_cpp.csv           (t,x,y,z samples)

#include <cmath>
#include <fstream>
#include <iostream>
#include <vector>

static inline double logistic(double x, double r) {
    return r * x * (1.0 - x);
}

void bifurcation_logistic(double rmin=2.5, double rmax=4.0, int Nr=6000,
                          int nTransient=1200, int nKeep=200, double x0=0.123456) {
    std::vector<double> rs(Nr), x(Nr, x0);
    for (int i = 0; i < Nr; ++i) {
        rs[i] = rmin + (rmax - rmin) * (double(i) / double(Nr - 1));
    }

    // Transient
    for (int k = 0; k < nTransient; ++k) {
        for (int i = 0; i < Nr; ++i) x[i] = logistic(x[i], rs[i]);
    }

    std::ofstream out("logistic_bifurcation_cpp.csv");
    out << "r,x\n";
    for (int k = 0; k < nKeep; ++k) {
        for (int i = 0; i < Nr; ++i) x[i] = logistic(x[i], rs[i]);
        for (int i = 0; i < Nr; ++i) out << rs[i] << "," << x[i] << "\n";
    }
    out.close();
}

struct Vec3 {
    double x, y, z;
    Vec3 operator+(const Vec3& o) const { return {x+o.x, y+o.y, z+o.z}; }
    Vec3 operator-(const Vec3& o) const { return {x-o.x, y-o.y, z-o.z}; }
    Vec3 operator*(double a) const { return {a*x, a*y, a*z}; }
};

Vec3 lorenz_rhs(const Vec3& s, double sigma, double rho, double beta) {
    return {
        sigma*(s.y - s.x),
        s.x*(rho - s.z) - s.y,
        s.x*s.y - beta*s.z
    };
}

Vec3 rk4_step(Vec3 y, double h, double sigma, double rho, double beta) {
    Vec3 k1 = lorenz_rhs(y, sigma, rho, beta);
    Vec3 k2 = lorenz_rhs(y + k1*(0.5*h), sigma, rho, beta);
    Vec3 k3 = lorenz_rhs(y + k2*(0.5*h), sigma, rho, beta);
    Vec3 k4 = lorenz_rhs(y + k3*h, sigma, rho, beta);
    return y + (k1 + k2*2.0 + k3*2.0 + k4) * (h/6.0);
}

void simulate_lorenz(double T=40.0, double h=0.005,
                     Vec3 y0={1.0,1.0,1.0},
                     double sigma=10.0, double rho=28.0, double beta=8.0/3.0,
                     double transient=5.0, int stride=4) {
    int N = int(T / h);
    double t = 0.0;
    Vec3 y = y0;

    std::ofstream out("lorenz_traj_cpp.csv");
    out << "t,x,y,z\n";

    for (int i = 0; i < N; ++i) {
        y = rk4_step(y, h, sigma, rho, beta);
        t += h;
        if (t >= transient && (i % stride == 0)) {
            out << t << "," << y.x << "," << y.y << "," << y.z << "\n";
        }
    }
    out.close();
}

int main() {
    bifurcation_logistic();
    simulate_lorenz();
    std::cout << "Done. Wrote CSV files." << std::endl;
    return 0;
}
      

File: Chapter20_Lesson2.java
This Java implementation mirrors the C++ workflow and writes CSV files for bifurcation points and Lorenz trajectory samples.

// Chapter20_Lesson2.java
// System Dynamics — Chapter 20 (Chaos, Complex Dynamics, and Computational Tools)
// Lesson 2: Bifurcation Diagrams, Period Doubling, and Strange Attractors
//
// Compile:
//   javac Chapter20_Lesson2.java
// Run:
//   java Chapter20_Lesson2
//
// Outputs:
//   logistic_bifurcation_java.csv
//   lorenz_traj_java.csv

import java.io.FileWriter;
import java.io.IOException;

public class Chapter20_Lesson2 {

    static double logistic(double x, double r) {
        return r * x * (1.0 - x);
    }

    static void bifurcationLogistic(double rmin, double rmax, int Nr,
                                   int nTransient, int nKeep, double x0) throws IOException {
        double[] rs = new double[Nr];
        double[] x = new double[Nr];
        for (int i = 0; i < Nr; i++) {
            rs[i] = rmin + (rmax - rmin) * (double)i / (double)(Nr - 1);
            x[i] = x0;
        }

        for (int k = 0; k < nTransient; k++) {
            for (int i = 0; i < Nr; i++) x[i] = logistic(x[i], rs[i]);
        }

        FileWriter fw = new FileWriter("logistic_bifurcation_java.csv");
        fw.write("r,x\n");
        for (int k = 0; k < nKeep; k++) {
            for (int i = 0; i < Nr; i++) x[i] = logistic(x[i], rs[i]);
            for (int i = 0; i < Nr; i++) fw.write(rs[i] + "," + x[i] + "\n");
        }
        fw.close();
    }

    static class Vec3 {
        double x, y, z;
        Vec3(double x, double y, double z){ this.x=x; this.y=y; this.z=z; }
        Vec3 add(Vec3 o){ return new Vec3(x+o.x, y+o.y, z+o.z); }
        Vec3 mul(double a){ return new Vec3(a*x, a*y, a*z); }
    }

    static Vec3 lorenzRhs(Vec3 s, double sigma, double rho, double beta) {
        return new Vec3(
            sigma*(s.y - s.x),
            s.x*(rho - s.z) - s.y,
            s.x*s.y - beta*s.z
        );
    }

    static Vec3 rk4Step(Vec3 y, double h, double sigma, double rho, double beta) {
        Vec3 k1 = lorenzRhs(y, sigma, rho, beta);
        Vec3 k2 = lorenzRhs(y.add(k1.mul(0.5*h)), sigma, rho, beta);
        Vec3 k3 = lorenzRhs(y.add(k2.mul(0.5*h)), sigma, rho, beta);
        Vec3 k4 = lorenzRhs(y.add(k3.mul(h)), sigma, rho, beta);
        return y.add(k1.add(k2.mul(2.0)).add(k3.mul(2.0)).add(k4).mul(h/6.0));
    }

    static void simulateLorenz(double T, double h, Vec3 y0,
                              double sigma, double rho, double beta,
                              double transient, int stride) throws IOException {
        int N = (int)(T / h);
        double t = 0.0;
        Vec3 y = y0;

        FileWriter fw = new FileWriter("lorenz_traj_java.csv");
        fw.write("t,x,y,z\n");
        for (int i = 0; i < N; i++) {
            y = rk4Step(y, h, sigma, rho, beta);
            t += h;
            if (t >= transient && (i % stride == 0)) {
                fw.write(t + "," + y.x + "," + y.y + "," + y.z + "\n");
            }
        }
        fw.close();
    }

    public static void main(String[] args) throws Exception {
        bifurcationLogistic(2.5, 4.0, 6000, 1200, 200, 0.123456);
        simulateLorenz(40.0, 0.005, new Vec3(1.0,1.0,1.0),
                       10.0, 28.0, 8.0/3.0, 5.0, 4);
        System.out.println("Done. Wrote CSV files.");
    }
}
      

File: Chapter20_Lesson2.m
This script generates: (i) a logistic bifurcation diagram in MATLAB and (ii) a Lorenz attractor plot via ode45. It also includes an optional section that programmatically constructs a simple Simulink model of the logistic map (Unit Delay + MATLAB Function), if Simulink is installed.

% Chapter20_Lesson2.m
% System Dynamics — Chapter 20 (Chaos, Complex Dynamics, and Computational Tools)
% Lesson 2: Bifurcation Diagrams, Period Doubling, and Strange Attractors
%
% Outputs:
%   - logistic_bifurcation_matlab.png
%   - lorenz_attractor_matlab.png
%
% Notes:
%   (1) The bifurcation diagram is built by sweeping r, discarding transients,
%       and plotting the last iterates for each r.
%   (2) The Lorenz trajectory is integrated with ode45 (continuous-time ODE).
%   (3) A minimal Simulink model can be generated programmatically (optional).

clear; clc; close all;

%% Part A — Logistic map bifurcation diagram
rmin = 2.5; rmax = 4.0; Nr = 6000;
nTransient = 1200; nKeep = 200;
rs = linspace(rmin, rmax, Nr);
x = 0.123456 * ones(1, Nr);

for k = 1:nTransient
    x = rs .* x .* (1 - x);
end

Rplot = zeros(1, Nr*nKeep);
Xplot = zeros(1, Nr*nKeep);
idx = 1;
for k = 1:nKeep
    x = rs .* x .* (1 - x);
    Rplot(idx:idx+Nr-1) = rs;
    Xplot(idx:idx+Nr-1) = x;
    idx = idx + Nr;
end

figure('Position',[100 100 900 450]);
plot(Rplot, Xplot, '.', 'MarkerSize', 1);
xlabel('r'); ylabel('x (asymptotic samples)');
title('Logistic map bifurcation diagram');
grid on;
saveas(gcf, 'logistic_bifurcation_matlab.png');

%% Part B — Lorenz attractor (ODE) with ode45
sigma = 10; rho = 28; beta = 8/3;
f = @(t, y) [ sigma*(y(2)-y(1));
              y(1)*(rho - y(3)) - y(2);
              y(1)*y(2) - beta*y(3) ];

tspan = [0 40];
y0 = [1; 1; 1];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t, Y] = ode45(f, tspan, y0, opts);

% discard transient
mask = t >= 5;
Y2 = Y(mask, :);

figure('Position',[120 120 700 600]);
plot3(Y2(:,1), Y2(:,2), Y2(:,3), 'LineWidth', 0.5);
xlabel('x'); ylabel('y'); zlabel('z');
title('Lorenz attractor trajectory (ode45)');
grid on; view(45, 25);
saveas(gcf, 'lorenz_attractor_matlab.png');

%% Part C — Optional: programmatically build a simple Simulink map model
% This creates a discrete-time logistic map with a Unit Delay and MATLAB Function.
% Requires Simulink installed. If not available, this section will safely skip.
try
    hasSimulink = license('test','Simulink');
catch
    hasSimulink = false;
end

if hasSimulink
    mdl = 'Chapter20_Lesson2_LogisticMap_Simulink';
    new_system(mdl); open_system(mdl);

    add_block('simulink/Sources/Constant', [mdl '/r'], 'Value', '3.7', ...
              'Position', [50 80 90 110]);
    add_block('simulink/Sources/Constant', [mdl '/x0'], 'Value', '0.2', ...
              'Position', [50 160 90 190]);

    add_block('simulink/Discrete/Unit Delay', [mdl '/UnitDelay'], ...
              'InitialCondition', '0.2', 'SampleTime', '1', ...
              'Position', [170 150 220 200]);

    add_block('simulink/User-Defined Functions/MATLAB Function', [mdl '/LogisticFcn'], ...
              'Position', [280 130 420 220]);

    % Set MATLAB Function code
    fcn = [ ...
        'function x_next = fcn(x, r)\n' ...
        '%#codegen\n' ...
        'x_next = r * x * (1 - x);\n' ...
        'end\n' ];
    set_param([mdl '/LogisticFcn'], 'Script', fcn);

    add_block('simulink/Sinks/Scope', [mdl '/Scope'], ...
              'Position', [480 150 520 190]);

    % Wire: x0 -> UnitDelay input (sets IC only; still connect for convenience)
    add_line(mdl, 'x0/1', 'UnitDelay/1');
    add_line(mdl, 'UnitDelay/1', 'LogisticFcn/1');
    add_line(mdl, 'r/1', 'LogisticFcn/2');
    add_line(mdl, 'LogisticFcn/1', 'UnitDelay/1', 'autorouting', 'on');
    add_line(mdl, 'UnitDelay/1', 'Scope/1');

    set_param(mdl, 'StopTime', '200', 'Solver', 'FixedStepDiscrete');
    save_system(mdl);
    % run with: sim(mdl)
end

disp('Done.');
      

File: Chapter20_Lesson2.nb
This notebook (text form) computes a logistic bifurcation diagram and integrates the Lorenz system using NDSolveValue.

(* Chapter20_Lesson2.nb
   System Dynamics — Chapter 20 (Chaos, Complex Dynamics, and Computational Tools)
   Lesson 2: Bifurcation Diagrams, Period Doubling, and Strange Attractors *)

Notebook[{

Cell["Chapter 20 — Lesson 2: Bifurcation Diagrams, Period Doubling, and Strange Attractors", "Title"],

Cell["A) Logistic map bifurcation diagram", "Section"],

Cell[BoxData[
 ToBoxes[
  Module[{rmin=2.5, rmax=4.0, Nr=4000, nTrans=1200, nKeep=200, x0=0.123456,
          rs, x, pts},
   rs = Subdivide[rmin, rmax, Nr-1];
   x = ConstantArray[x0, Nr];
   Do[x = rs*x*(1 - x), {nTrans}];
   pts = Reap[
      Do[
        x = rs*x*(1 - x);
        Sow[Transpose[{rs, x}]];
      , {nKeep}]
   ][[2,1]];
   ListPlot[Flatten[pts, 1], PlotStyle -> PointSize[0.001],
            AxesLabel -> {"r","x"}, PlotRange -> All,
            PlotLabel -> "Logistic map bifurcation diagram"]
  ]
 ]
], "Input"],

Cell["B) Lorenz strange attractor (continuous-time ODE)", "Section"],

Cell[BoxData[
 ToBoxes[
  Module[{sigma=10, rho=28, beta=8/3, sol, tmax=40, y0={1,1,1}, traj},
   sol = NDSolveValue[
     {
      x'[t] == sigma (y[t]-x[t]),
      y'[t] == x[t] (rho - z[t]) - y[t],
      z'[t] == x[t] y[t] - beta z[t],
      x[0]==y0[[1]], y[0]==y0[[2]], z[0]==y0[[3]]
     },
     {x,y,z}, {t,0,tmax},
     Method -> {"StiffnessSwitching"}, AccuracyGoal -> 12, PrecisionGoal -> 12
   ];
   traj = Table[{sol[[1]][t], sol[[2]][t], sol[[3]][t]}, {t,5,tmax,0.01}];
   Graphics3D[{Thickness[0.001], Line[traj]},
     Axes -> True, AxesLabel -> {"x","y","z"},
     PlotLabel -> "Lorenz attractor trajectory (sampled)"]
  ]
 ]
], "Input"]

}]
      

Problem 1 (Fixed points and stability of the logistic map): For \( x_{n+1} = r x_n(1-x_n) \), find all fixed points and determine their stability intervals (for \( r \ge 0 \)).

Solution: Fixed points satisfy \( x = r x(1-x) \), so \( x_1^\star=0 \) and \( x_2^\star = 1-\frac{1}{r} \) for \( r\neq 0 \). Since \( f_x(x;r)=r(1-2x) \), we have \( f_x(x_1^\star;r)=r \), hence \( |r|<1 \Rightarrow 0 \le r < 1 \) stable. For \( x_2^\star \), \( f_x(x_2^\star;r)=2-r \), so \( |2-r|<1 \Rightarrow 1 < r < 3 \) stable.

Problem 2 (Explicit period-2 orbit): Show that the non-fixed solutions of \( f(f(x;r);r)=x \) are \( x_{\pm}^{(2)}(r) = \frac{r+1 \pm \sqrt{(r-3)(r+1)}}{2r} \), and determine when they are real.

Solution: Expanding and factoring \( f(f(x;r);r)-x \) separates fixed points from the 2-cycle polynomial \( r^2 x^2 - r^2 x - r x + r + 1 = 0 \). Solving the quadratic yields the stated formula. Reality requires \( (r-3)(r+1)\ge 0 \). For \( r\ge 0 \), this means \( r \ge 3 \).

Problem 3 (Stability interval of the 2-cycle): Let \( x_1, x_2 \) be the logistic map period-2 orbit. Prove that the period-2 multiplier satisfies \( m_2(r)=f_x(x_1;r)f_x(x_2;r) = -r^2+2r+4 \), and deduce the stability range for \( r\ge 0 \).

Solution: The linearized evolution over two steps is \( e_{n+2} \approx f_x(x_2;r)f_x(x_1;r)e_n \). For the logistic map, substituting the explicit 2-cycle points into the derivative product simplifies to \( -r^2+2r+4 \). Stability is \( |m_2(r)|<1 \) which gives \( 3 < r < 1+\sqrt{6} \) for \( r\ge 0 \). The upper endpoint is the next flip bifurcation (2→4).

Problem 4 (How a bifurcation diagram is computed): Suppose you want a bifurcation diagram for \( x_{n+1}=f(x_n;\mu) \). Specify a correct algorithm (including transient handling) and explain why the diagram may depend on the initial condition.

Solution: Choose a grid \( \{\mu_i\}_{i=1}^M \). For each \( \mu_i \), initialize \( x_0 \) and iterate \( N_{\text{trans}} \) steps to remove transients. Then record the next \( N_{\text{keep}} \) iterates \( x_{N_{\text{trans}}+1},\dots,x_{N_{\text{trans}}+N_{\text{keep}}} \) and plot pairs \( (\mu_i,x_n) \). Dependence on initial condition arises when multiple attractors coexist (multistability): different basins lead to different asymptotic sets, so a single diagram can show only the attractor reached by the chosen initialization rule.

Problem 5 (Building an ODE bifurcation diagram via sampling): For \( \dot{\mathbf{x}}=\mathbf{F}(\mathbf{x};\mu) \), propose a method to construct a bifurcation diagram versus \( \mu \) without explicitly computing equilibria. Use the Lorenz system as an example.

Solution: For each \( \mu \), numerically integrate long enough to discard transients. Then sample a scalar observable at discrete events: (i) maxima of \( x(t) \), (ii) intersections with a Poincaré section (e.g., \( z = z_0 \) with \( \dot{z} > 0 \)), or (iii) stroboscopic sampling \( t=kT \). Plot sampled values versus \( \mu \). For Lorenz, a common choice is to sweep \( \rho \), integrate, collect local maxima of \( x(t) \) after transient, and scatter-plot (rho, x_max). This reveals periodic windows, banded chaos, and crises.

We defined bifurcation diagrams as parameter-indexed portraits of asymptotic dynamics, derived the linear stability criterion \( |f_x| < 1 \) for 1D maps, and used it to rigorously explain the first steps of the period-doubling cascade in the logistic map: fixed-point stability \( 1 < r < 3 \), flip at \( r=3 \), explicit 2-cycle points, and the next doubling at \( r=1+\sqrt{6} \). We connected discrete bifurcation diagrams to continuous-time systems through sampling and Poincaré return maps, motivating why strange attractors (e.g., Lorenz) can be studied using discrete-time tools. Next lesson formalizes sensitive dependence and Lyapunov exponents as quantitative measures of chaos.