In this chapter we study deterministic systems whose trajectories can look “random-like.” We keep three notions distinct:

• Nonlinear map (discrete-time): a recurrence \( x_{n+1}=f(x_n) \) on a state space (often an interval or manifold).
• Nonlinear flow (continuous-time): an ODE \( \dot{\mathbf{x}}=\mathbf{F}(\mathbf{x}) \) generating a trajectory \( \mathbf{x}(t) \).
• Chaos (informal here): bounded deterministic motion with strong sensitivity to initial conditions and complicated aperiodic structure. Formal quantification (Lyapunov exponents) is deferred to Lesson 3.

For maps, “sensitivity” can be hinted by linearization: for small perturbations \( \delta_n \), one often has \( \delta_{n+1}\approx f'(x_n)\,\delta_n \), so that

\[ |\delta_n| \approx |\delta_0| \prod_{k=0}^{n-1} |f'(x_k)|. \]

For flows, local behavior near equilibria uses the Jacobian \( \mathbf{J}(\mathbf{x}^\star)=\frac{\partial \mathbf{F}}{\partial \mathbf{x}}(\mathbf{x}^\star) \) and eigenvalues of \( \mathbf{J} \).

flowchart TD
  A["Choose model type"] --> B["Discrete map: x[n+1]=f(x[n])"]
  A --> C["Continuous flow: xdot = F(x)"]
  B --> D["Iterate: x[0] -> x[1] -> ..."]
  C --> E["Integrate numerically: x(t0) -> x(t)"]
  D --> F["Analyze: fixed points, stability via |f'(x*)|"]
  E --> G["Analyze: equilibria, Jacobian eigenvalues"]
  F --> H["Observe complexity: aperiodic bounded orbits"]
  G --> H
  H --> I["Compare nearby initial states (sensitivity)"]
        

The logistic map is the scalar discrete-time system \( x_{n+1}=f(x_n)=r\,x_n(1-x_n) \), typically studied on \( x_n\in[0,1] \) with parameter \( r\in(0,4] \). As a stylized population model, \( x_n \) is a normalized population, \( r x_n \) is reproduction, and \( (1-x_n) \) encodes limited resources.

A key structural fact is that for \( 0 < r \le 4 \), the interval \( [0,1] \) is forward invariant (with a mild caveat at boundaries). This bounds the dynamics and allows nontrivial long-term behavior.

2.1 Invariance of the interval \( [0,1] \)

Suppose \( x\in[0,1] \) and \( 0 < r \le 4 \). Since \( x(1-x)\ge 0 \), we have \( f(x)\ge 0 \). Also, the concave parabola \( x(1-x) \) satisfies \( 0 \le x(1-x)\le \tfrac{1}{4} \) with maximum at \( x=\tfrac{1}{2} \). Hence

\[ 0 \le f(x)=r\,x(1-x) \le r\cdot \tfrac{1}{4} \le 1, \]

which proves \( f([0,1])\subseteq[0,1] \). Therefore, trajectories starting in \( [0,1] \) remain in \( [0,1] \).

A fixed point \( x^\star \) satisfies \( x^\star=f(x^\star) \). For the logistic map:

\[ x^\star = r x^\star(1-x^\star) \;\Longrightarrow\; r(x^\star)^2 - (r-1)x^\star = 0 \;\Longrightarrow\; x^\star\in\left\{0,\;1-\frac{1}{r}\right\}. \]

Local stability of a fixed point for a 1D map is governed by the derivative: if \( |f'(x^\star)| < 1 \) then the fixed point is locally asymptotically stable; if \( |f'(x^\star)| > 1 \) it is unstable.

Here \( f'(x)=r(1-2x) \). Evaluate:

3.1 Stability of \( x^\star=0 \)

\[ f'(0)=r,\quad |f'(0)| < 1 \iff 0 < r < 1. \]

Thus \( x^\star=0 \) is stable for \( 0 < r < 1 \) and unstable for \( r > 1 \).

3.2 Stability of \( x^\star=1-\tfrac{1}{r} \)

\[ f'\!\left(1-\frac{1}{r}\right)=r\left(1-2+ \frac{2}{r}\right)= -r+2, \quad | -r+2 | < 1 \iff 1 < r < 3. \]

Hence the nonzero equilibrium is stable for \( 1 < r < 3 \), and loses stability at \( r=3 \). What happens afterward (period-doubling cascades and bifurcations) is the focus of Lesson 2.

Key takeaway: already at the level of a single derivative test, a nonlinear map can undergo qualitative changes as a parameter crosses thresholds (a first glimpse of bifurcation behavior).

On \( [0,1] \), the logistic map is a unimodal parabola. For larger \( r \), the slope magnitude \( |f'(x)| \) can exceed 1 over substantial portions of the interval, meaning nearby states can separate (stretching). However, because \( [0,1] \) is invariant, trajectories must also “fold” back into the interval. The interplay of stretching + folding under repeated iteration is the core geometric intuition behind chaotic sets in 1D maps.

A standard visualization is the cobweb plot, showing repeated application of \( f \) and the diagonal \( y=x \). This is computationally light and reveals convergence, periodicity, or irregular wandering depending on \( r \).

The Lorenz system is the 3D nonlinear ODE with parameters \( \sigma, \rho, \beta \):

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

For the classical chaotic regime one often uses \( \sigma=10 \), \( \rho=28 \), \( \beta=\tfrac{8}{3} \). This model emerged from simplified convection equations and became a paradigm for deterministic chaos in flows.

5.1 Equilibria (fixed points of the flow)

Equilibria satisfy \( \dot{x}=\dot{y}=\dot{z}=0 \). From \( \dot{x}=0 \Rightarrow y=x \). Then \( \dot{y}=0 \Rightarrow x(\rho-z)-x=0 \), so either \( x=0 \) or \( z=\rho-1 \). Finally \( \dot{z}=0 \Rightarrow xy=\beta z \). Hence:

• \( \mathbf{x}^\star_0=(0,0,0) \)
• If \( \rho > 1 \), two additional equilibria: \( \mathbf{x}^\star_{\pm}=\left(\pm\sqrt{\beta(\rho-1)},\;\pm\sqrt{\beta(\rho-1)},\;\rho-1\right) \)

5.2 Dissipativity via volume contraction (a rigorous quick proof)

Consider the divergence of the vector field \( \mathbf{F}=(\dot{x},\dot{y},\dot{z}) \):

\[ \nabla\cdot \mathbf{F} =\frac{\partial \dot{x}}{\partial x}+\frac{\partial \dot{y}}{\partial y}+\frac{\partial \dot{z}}{\partial z} =(-\sigma)+(-1)+(-\beta)=-(\sigma+1+\beta). \]

Since \( \sigma+1+\beta > 0 \), the divergence is strictly negative everywhere: volumes in phase space contract exponentially under the flow. This implies the existence of an attracting set of zero Lebesgue volume (an “attractor”), which in the classical parameter regime is a strange attractor.

Let \( \mathbf{x}=[x\;y\;z]^\top \). The Jacobian is \( \mathbf{J}(\mathbf{x})=\frac{\partial \mathbf{F}}{\partial \mathbf{x}} \):

\[ \mathbf{J}(x,y,z)= \begin{bmatrix} -\sigma & \sigma & 0 \\ \rho-z & -1 & -x \\ y & x & -\beta \end{bmatrix}. \]

6.1 Stability of the origin

At \( (0,0,0) \):

\[ \mathbf{J}(0,0,0)= \begin{bmatrix} -\sigma & \sigma & 0 \\ \rho & -1 & 0 \\ 0 & 0 & -\beta \end{bmatrix}. \]

One eigenvalue is \( -\beta \) (stable in the \( z \)-direction). The remaining 2D block has characteristic polynomial:

\[ \lambda^2 + (\sigma+1)\lambda + \sigma(1-\rho)=0. \]

If \( \rho < 1 \), all coefficients are positive, hence both roots have negative real parts, so the origin is locally asymptotically stable. If \( \rho > 1 \), the constant term \( \sigma(1-\rho) \) is negative, so one eigenvalue is positive: the origin becomes unstable.

(The local stability of \( \mathbf{x}^\star_\pm \) depends on parameters and involves a Hopf bifurcation, which we will not develop here; it belongs naturally with the bifurcation discussion in Lesson 2.)

Discrete maps and ODEs require different numerical workflows:

• Map: iterate exactly (up to floating-point arithmetic) \( x_{n+1}=f(x_n) \).
• ODE: approximate the flow using a time-stepper (e.g., RK4, ode45, solve_ivp), controlling local error and step size.

flowchart TD
  M1["Logistic map"] --> M2["Pick r in (0,4], choose x0 in [0,1]"]
  M2 --> M3["Iterate: x[n+1]=r*x[n]*(1-x[n])"]
  M3 --> M4["Diagnostics: time series, cobweb, return map"]
  O1["Lorenz ODE"] --> O2["Pick parameters sigma,rho,beta and x(0)"]
  O2 --> O3["Integrate: RK4 / ode45 / solve_ivp"]
  O3 --> O4["Diagnostics: time series, 3D trajectory, sensitivity"]
        

In later lessons we will add systematic tools: bifurcation diagrams (Lesson 2) and Lyapunov exponents (Lesson 3).

We implement (i) logistic-map iteration and cobweb plotting, (ii) Lorenz integration via a from-scratch RK4, and (iii) a sensitivity demonstration by integrating two nearby initial conditions.

File: Chapter20_Lesson1.py

# Chapter20_Lesson1.py
"""
System Dynamics — Chapter 20, Lesson 1
Nonlinear Maps and Continuous-Time Chaotic Systems
(Logistic Map, Lorenz System)
"""

import numpy as np
import matplotlib.pyplot as plt

def logistic_map(r: float, x0: float, n: int) -> np.ndarray:
    x = np.empty(n + 1, dtype=float)
    x[0] = x0
    for k in range(n):
        x[k + 1] = r * x[k] * (1.0 - x[k])
    return x

def cobweb_data(f, x0: float, n: int):
    xs, ys = [], []
    x = x0
    for _ in range(n):
        y = f(x)
        xs.extend([x, x]); ys.extend([x, y])   # vertical
        xs.extend([x, y]); ys.extend([y, y])   # horizontal
        x = y
    return np.array(xs), np.array(ys)

def plot_logistic_demo(r=3.8, x0=0.2, n=80):
    f = lambda x: r * x * (1.0 - x)
    x = logistic_map(r, x0, n)

    plt.figure()
    plt.plot(np.arange(n + 1), x)
    plt.xlabel("n"); plt.ylabel("x_n")
    plt.title(f"Logistic map time series (r={r}, x0={x0})")
    plt.grid(True)

    plt.figure()
    grid = np.linspace(0, 1, 600)
    plt.plot(grid, f(grid), label="f(x)")
    plt.plot(grid, grid, label="y=x")
    xs, ys = cobweb_data(f, x0, n=40)
    plt.plot(xs, ys, linewidth=1.0)
    plt.xlim(0, 1); plt.ylim(0, 1)
    plt.xlabel("x_n"); plt.ylabel("x_{n+1}")
    plt.title(f"Cobweb plot (r={r}, x0={x0})")
    plt.legend(); plt.grid(True)

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

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

def integrate_lorenz_rk4(x0=(1.0,1.0,1.0), t0=0.0, tf=40.0, h=0.005,
                         sigma=10.0, rho=28.0, beta=8.0/3.0):
    n = int(np.ceil((tf - t0) / h))
    t = np.linspace(t0, t0 + n*h, n + 1)
    X = np.empty((n + 1, 3))
    X[0, :] = np.array(x0, dtype=float)
    for k in range(n):
        X[k + 1, :] = rk4_step(lorenz_rhs, t[k], X[k, :], h,
                               sigma=sigma, rho=rho, beta=beta)
    return t, X

def plot_lorenz_demo():
    t, X = integrate_lorenz_rk4()
    x, y, z = X[:,0], X[:,1], X[:,2]

    plt.figure()
    plt.plot(t, x, label="x(t)")
    plt.plot(t, y, label="y(t)")
    plt.plot(t, z, label="z(t)")
    plt.xlabel("t"); plt.ylabel("states")
    plt.title("Lorenz system states vs time (RK4)")
    plt.legend(); plt.grid(True)

    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    ax.plot(x, y, z, linewidth=0.7)
    ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z")
    ax.set_title("Lorenz attractor (RK4)")

def sensitivity_demo_lorenz(eps=1e-9):
    x0 = np.array([1.0, 1.0, 1.0])
    x0b = x0 + np.array([eps, 0.0, 0.0])

    t, Xa = integrate_lorenz_rk4(x0=tuple(x0))
    _, Xb = integrate_lorenz_rk4(x0=tuple(x0b))
    d = np.linalg.norm(Xa - Xb, axis=1)

    plt.figure()
    plt.semilogy(t, d + 1e-30)
    plt.xlabel("t"); plt.ylabel("||delta(t)||")
    plt.title(f"Sensitivity in Lorenz system (eps={eps})")
    plt.grid(True)

if __name__ == "__main__":
    plot_logistic_demo()
    plot_lorenz_demo()
    sensitivity_demo_lorenz()
    plt.show()
      

Exercise file: Chapter20_Lesson1_Ex1.py

# Chapter20_Lesson1_Ex1.py
import numpy as np
import matplotlib.pyplot as plt

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

def iterate(r, x0, n):
    x = np.empty(n + 1)
    x[0] = x0
    for k in range(n):
        x[k + 1] = logistic(r, x[k])
    return x

if __name__ == "__main__":
    r = 3.9
    n = 200
    x0 = 0.2000000000
    y0 = 0.2000000001

    x = iterate(r, x0, n)
    y = iterate(r, y0, n)
    d = np.abs(x - y)

    plt.figure()
    plt.semilogy(np.arange(n + 1), d + 1e-30)
    plt.xlabel("n"); plt.ylabel("|x_n - y_n|")
    plt.title("Log-distance between nearby logistic-map trajectories")
    plt.grid(True)
    plt.show()
      

The C++ implementation writes logistic.csv and lorenz.csv so you can plot in any tool.

File: Chapter20_Lesson1.cpp

// Chapter20_Lesson1.cpp
#include <cmath>
#include <fstream>
#include <iostream>
#include <vector>

static std::vector<double> logistic_map(double r, double x0, int n) {
    std::vector<double> x(n + 1);
    x[0] = x0;
    for (int k = 0; k < n; ++k) x[k + 1] = r * x[k] * (1.0 - x[k]);
    return x;
}

struct Vec3 { double x, y, z; };

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

static Vec3 add(const Vec3& a, const Vec3& b) { return Vec3{a.x+b.x,a.y+b.y,a.z+b.z}; }
static Vec3 mul(double c, const Vec3& a) { return Vec3{c*a.x,c*a.y,c*a.z}; }

static Vec3 rk4_step(double h, const Vec3& x, double sigma, double rho, double beta) {
    Vec3 k1 = lorenz_rhs(x, sigma, rho, beta);
    Vec3 k2 = lorenz_rhs(add(x, mul(0.5*h, k1)), sigma, rho, beta);
    Vec3 k3 = lorenz_rhs(add(x, mul(0.5*h, k2)), sigma, rho, beta);
    Vec3 k4 = lorenz_rhs(add(x, mul(h, k3)), sigma, rho, beta);
    return add(x, mul(h/6.0, add(add(k1, mul(2.0,k2)), add(mul(2.0,k3), k4))));
}

int main() {
    // Logistic map
    double r = 3.8, x0 = 0.2; int n = 2000;
    auto x = logistic_map(r, x0, n);
    std::ofstream flog("logistic.csv");
    flog << "n,x\n";
    for (int k = 0; k <= n; ++k) flog << k << "," << x[k] << "\n";
    flog.close();

    // Lorenz
    double sigma = 10.0, rho = 28.0, beta = 8.0/3.0;
    double t = 0.0, tf = 40.0, h = 0.005;
    int steps = (int)std::ceil((tf - t)/h);

    std::ofstream flor("lorenz.csv");
    flor << "t,x,y,z\n";
    Vec3 s{1.0,1.0,1.0};
    flor << t << "," << s.x << "," << s.y << "," << s.z << "\n";
    for (int k = 0; k < steps; ++k) {
        s = rk4_step(h, s, sigma, rho, beta);
        t += h;
        flor << t << "," << s.x << "," << s.y << "," << s.z << "\n";
    }
    flor.close();

    std::cout << "Wrote logistic.csv and lorenz.csv\n";
    return 0;
}
      

File: Chapter20_Lesson1.java

// Chapter20_Lesson1.java
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;

public class Chapter20_Lesson1 {

    static double[] logisticMap(double r, double x0, int n) {
        double[] x = new double[n + 1];
        x[0] = x0;
        for (int k = 0; k < n; k++) x[k + 1] = r * x[k] * (1.0 - x[k]);
        return x;
    }

    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 b) { return new Vec3(this.x + b.x, this.y + b.y, this.z + b.z); }
        Vec3 mul(double s) { return new Vec3(s * this.x, s * this.y, s * this.z); }
    }

    static Vec3 lorenzRHS(Vec3 s, double sigma, double rho, double beta) {
        double dx = sigma * (s.y - s.x);
        double dy = s.x * (rho - s.z) - s.y;
        double dz = s.x * s.y - beta * s.z;
        return new Vec3(dx, dy, dz);
    }

    static Vec3 rk4Step(Vec3 x, double h, double sigma, double rho, double beta) {
        Vec3 k1 = lorenzRHS(x, sigma, rho, beta);
        Vec3 k2 = lorenzRHS(x.add(k1.mul(0.5*h)), sigma, rho, beta);
        Vec3 k3 = lorenzRHS(x.add(k2.mul(0.5*h)), sigma, rho, beta);
        Vec3 k4 = lorenzRHS(x.add(k3.mul(h)), sigma, rho, beta);

        return x.add(k1.mul(h/6.0))
                .add(k2.mul(h/3.0))
                .add(k3.mul(h/3.0))
                .add(k4.mul(h/6.0));
    }

    public static void main(String[] args) throws IOException {
        // Logistic
        double r = 3.8, x0 = 0.2; int n = 2000;
        double[] x = logisticMap(r, x0, n);
        try (BufferedWriter w = new BufferedWriter(new FileWriter("logistic.csv"))) {
            w.write("n,x\n");
            for (int k = 0; k <= n; k++) w.write(k + "," + x[k] + "\n");
        }

        // Lorenz
        double sigma = 10.0, rho = 28.0, beta = 8.0/3.0;
        double t = 0.0, tf = 40.0, h = 0.005;
        int steps = (int)Math.ceil((tf - t)/h);

        Vec3 s = new Vec3(1.0, 1.0, 1.0);
        try (BufferedWriter w = new BufferedWriter(new FileWriter("lorenz.csv"))) {
            w.write("t,x,y,z\n");
            w.write(t + "," + s.x + "," + s.y + "," + s.z + "\n");
            for (int k = 0; k < steps; k++) {
                s = rk4Step(s, h, sigma, rho, beta);
                t += h;
                w.write(t + "," + s.x + "," + s.y + "," + s.z + "\n");
            }
        }

        System.out.println("Wrote logistic.csv and lorenz.csv");
    }
}
      

File: Chapter20_Lesson1.m

% Chapter20_Lesson1.m
clear; close all; clc;

%% Logistic map
r = 3.8; x0 = 0.2; N = 80;
x = zeros(N+1,1); x(1) = x0;
for k = 1:N
    x(k+1) = r*x(k)*(1 - x(k));
end
figure; plot(0:N, x, 'LineWidth', 1.2);
grid on; xlabel('n'); ylabel('x_n');
title(sprintf('Logistic map time series (r=%.3g, x0=%.3g)', r, x0));

% Cobweb
f = @(u) r*u.*(1-u);
gridx = linspace(0,1,600);
figure; plot(gridx, f(gridx), 'LineWidth', 1.2); hold on;
plot(gridx, gridx, 'LineWidth', 1.2);
xk = x0;
for k = 1:40
    yk = f(xk);
    plot([xk xk], [xk yk], 'k-');
    plot([xk yk], [yk yk], 'k-');
    xk = yk;
end
xlim([0 1]); ylim([0 1]); grid on;
xlabel('x_n'); ylabel('x_{n+1}');
title(sprintf('Cobweb plot (r=%.3g, x0=%.3g)', r, x0));

%% Lorenz via ode45
sigma = 10; rho = 28; beta = 8/3;
lorenz = @(t,s) [ sigma*(s(2)-s(1));
                  s(1)*(rho - s(3)) - s(2);
                  s(1)*s(2) - beta*s(3) ];
tspan = [0 40];
s0 = [1; 1; 1];
opts = odeset('RelTol',1e-9,'AbsTol',1e-12,'MaxStep',0.01);
[t,S] = ode45(lorenz, tspan, s0, opts);

figure; plot(t, S(:,1), t, S(:,2), t, S(:,3), 'LineWidth', 1.0);
grid on; xlabel('t'); ylabel('states');
title('Lorenz states vs time (ode45)'); legend('x','y','z');

figure; plot3(S(:,1), S(:,2), S(:,3), 'LineWidth', 0.7);
grid on; xlabel('x'); ylabel('y'); zlabel('z');
title('Lorenz attractor (ode45)');

%% Sensitivity demo
eps = 1e-9;
s0b = s0 + [eps; 0; 0];
[t2,S2] = ode45(lorenz, tspan, s0b, opts);
d = vecnorm((S - S2).', 2).';
figure; semilogy(t, d + 1e-30, 'LineWidth', 1.0);
grid on; xlabel('t'); ylabel('||delta(t)||');
title(sprintf('Sensitivity in Lorenz system (eps=%.1e)', eps));

%% Optional Simulink model builder is included in the downloadable file (commented).
      

File: Chapter20_Lesson1.nb

(* Chapter20_Lesson1.nb (Wolfram Language content) *)

logisticIter[r_, x0_, n_] := NestList[(r # (1 - #)) &, x0, n];

r = 3.8;
x = logisticIter[r, 0.2, 200];
ListLinePlot[x]

sigma = 10; rho = 28; beta = 8/3;

lorenz = {
  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] == 1, y[0] == 1, z[0] == 1
};

sol = NDSolve[lorenz, {x, y, z}, {t, 0, 40}, MaxStepSize -> 0.01];

ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 40}]
      

Problem 1 (Invariant interval for the logistic map): Show that for \( 0 < r \le 4 \), if \( x_0\in[0,1] \) then \( x_n\in[0,1] \) for all \( n\ge 0 \).

Solution: For \( x\in[0,1] \), \( 0\le x(1-x)\le \tfrac{1}{4} \). Thus \( 0\le f(x)=r x(1-x)\le r/4\le 1 \). Hence \( f([0,1])\subseteq[0,1] \) and induction yields the claim.

Problem 2 (Fixed points and stability thresholds): For \( x_{n+1}=r x_n(1-x_n) \), find fixed points and determine for which \( r \) each fixed point is locally stable.

Solution: Fixed points satisfy \( x=r x(1-x) \), so \( x^\star\in\{0,\;1-\tfrac{1}{r}\} \). With \( f'(x)=r(1-2x) \), \( |f'(0)|=|r| < 1 \) gives stability for \( 0 < r < 1 \). For \( 1-\tfrac{1}{r} \), \( f'(1-\tfrac{1}{r})=2-r \); stability requires \( |2-r| < 1 \), i.e. \( 1 < r < 3 \).

Problem 3 (Derivative-based perturbation growth in maps): Let \( x_{n+1}=f(x_n) \) where \( f \) is differentiable. Show that for a small perturbation \( \delta_0 \), the linearized perturbation evolves as \( \delta_n \approx \delta_0 \prod_{k=0}^{n-1} f'(x_k) \).

Solution: Linearize: \( x_{n+1}+\delta_{n+1}=f(x_n+\delta_n)\approx f(x_n)+f'(x_n)\delta_n \). Subtract \( x_{n+1}=f(x_n) \) to get \( \delta_{n+1}\approx f'(x_n)\delta_n \). Iterating yields \( \delta_n \approx \delta_0 \prod_{k=0}^{n-1} f'(x_k) \).

Problem 4 (Lorenz equilibria): Find all equilibria of the Lorenz system and state the condition for the nonzero equilibria to exist.

Solution: From \( \dot{x}=0 \Rightarrow y=x \). Then \( \dot{y}=0 \Rightarrow x(\rho-z)-x=0 \Rightarrow x=0 \) or \( z=\rho-1 \). If \( x=0 \), then \( y=0 \) and \( \dot{z}=0 \Rightarrow z=0 \), giving \( (0,0,0) \). If \( z=\rho-1 \), then \( \dot{z}=0 \Rightarrow x^2=\beta(\rho-1) \), so equilibria \( \left(\pm\sqrt{\beta(\rho-1)},\;\pm\sqrt{\beta(\rho-1)},\;\rho-1\right) \) exist iff \( \rho > 1 \).

Problem 5 (Dissipativity by divergence): Compute \( \nabla\cdot\mathbf{F} \) for the Lorenz vector field and conclude that phase-space volumes contract.

Solution: With \( \dot{x}=\sigma(y-x),\;\dot{y}=x(\rho-z)-y,\;\dot{z}=xy-\beta z \), we have \( \partial \dot{x}/\partial x=-\sigma \), \( \partial \dot{y}/\partial y=-1 \), \( \partial \dot{z}/\partial z=-\beta \). Therefore

\[ \nabla\cdot\mathbf{F}=-(\sigma+1+\beta) < 0. \]

By Liouville’s formula, an infinitesimal volume element \( V(t) \) evolves as \( \dot{V}=(\nabla\cdot\mathbf{F})V \), so \( V(t)=V(0)\exp\!\left(-(\sigma+1+\beta)t\right) \), proving exponential volume contraction.

We introduced two foundational chaotic-system models: the logistic map (discrete-time) and the Lorenz system (continuous-time). For the logistic map we proved invariance of \( [0,1] \), computed fixed points, and derived stability ranges using \( |f'(x^\star)| \). For the Lorenz system we derived equilibria, computed the Jacobian, showed when the origin becomes unstable, and proved global dissipativity via negative divergence. We also implemented iteration and RK4/ode45-style integration pipelines across multiple languages.