mathematicskit.fractals_chaos#

Lyapunov exponent estimation for 1D maps and flows; fractal (box-counting) dimension estimation; Mandelbrot and Julia set generation (Numba-accelerated escape-time iteration); iterated function systems (Barnsley fern, Sierpinski triangle/carpet) via the chaos game; elementary cellular automata (Wolfram rule numbering) and Conway’s Game of Life. The logistic map’s own Feigenbaum bifurcation cascade lives in mathematicskit.ode_dynamics.

mathematicskit.fractals_chaos: discrete dynamical systems and fractal geometry.

The logistic map’s Feigenbaum bifurcation cascade lives in mathematicskit.ode_dynamics (LogisticMap, bifurcation_diagram()); this domain covers what’s distinctly fractal: Lyapunov exponent estimation for 1D maps and flows (quantifying the chaos that cascade ends in); fractal (box-counting) dimension estimation; Mandelbrot and Julia set generation; iterated function systems (Barnsley fern, Sierpinski triangle/carpet); classical fractal curves (Weierstrass, Koch, Hilbert), Richardson’s divider length, and the similarity dimension; Lindenmayer systems with turtle graphics; the Hénon map; diffusion-limited aggregation; and elementary (Wolfram-numbered) cellular automata plus Conway’s Game of Life.

class mathematicskit.fractals_chaos.BarnsleyFern[source]#

Bases: IteratedFunctionSystem

Barnsley’s fern: four affine maps with Barnsley’s classic probabilities.

The four maps (stem, successively smaller leaflets) are Barnsley’s original 1988 coefficients. See Barnsley, Fractals Everywhere, 2nd ed., Ch. 3, Table III.

Examples

>>> import numpy as np
>>> fern = BarnsleyFern()
>>> points = fern.generate(2000, seed=0)
>>> points.shape
(2000, 2)
>>> # The fern fits within its known bounding box.
>>> bool(np.all(points[:, 0] >= -3.0) and np.all(points[:, 0] <= 3.0))
True
generate(n_points, seed=0)[source]#

Run the chaos game for n_points iterations.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n_points, 2)

class mathematicskit.fractals_chaos.BoxCountingResult(box_sizes, box_counts, dimension=0.0, extra=<factory>)[source]#

Bases: object

Container for a box-counting fractal-dimension estimate.

Parameters:
box_counts: ndarray#

Number of occupied boxes at each size.

Type:

ndarray, int

box_sizes: ndarray#

Box edge lengths used.

Type:

ndarray

dimension: float = 0.0#

Estimated box-counting dimension (slope of log(count) vs. log(1/size)).

Type:

float

extra: dict#

Free-form diagnostics slot (e.g. the linear fit’s residuals).

Type:

dict

class mathematicskit.fractals_chaos.CellularAutomaton[source]#

Bases: ABC

Common base for cellular automata (1D elementary or 2D Game of Life).

run(n_steps)[source]#

Run for n_steps generations, recording every state.

Parameters:

n_steps (int)

Return type:

ndarray

Returns:

ndarray, shape (n_steps + 1, …) – The initial state followed by each generation.

state: ndarray#
abstractmethod step()[source]#

Advance the automaton by one generation and return the new state.

Return type:

ndarray

class mathematicskit.fractals_chaos.ElementaryCA(rule, width, initial=None)[source]#

Bases: CellularAutomaton

1D elementary cellular automaton, Wolfram rule numbering (0-255).

Each cell’s next state is a function of its own and its two neighbors’ current states – 8 possible 3-cell neighborhoods, so \(2^8 = 256\) possible rules, each identified by the 8-bit number whose bits give the output for each neighborhood (MSB-first: 111, 110, 101, 100, 011, 010, 001, 000), Wolfram’s numbering convention. See Wolfram, A New Kind of Science, 2002, Ch. 2.

Parameters:
  • rule (int) – Wolfram rule number, 0 <= rule <= 255.

  • width (int) – Number of cells (periodic boundary).

  • initial (ndarray, shape (width,), optional) – Initial state (0s and 1s); defaults to a single 1 at the center cell, the standard way to display a rule’s characteristic pattern.

Examples

>>> import numpy as np
>>> # Rule 90 is the XOR of a cell's two neighbors, producing a
>>> # discrete Sierpinski triangle from a single seed cell.
>>> ca = ElementaryCA(rule=90, width=7)
>>> ca.state.tolist()
[0, 0, 0, 1, 0, 0, 0]
>>> ca.step().tolist()
[0, 0, 1, 0, 1, 0, 0]
step()[source]#

Advance the automaton by one generation and return the new state.

Return type:

ndarray

class mathematicskit.fractals_chaos.EscapeTimeResult(iterations, extent=(-2.0, 1.0, -1.5, 1.5), max_iter=100)[source]#

Bases: object

Container for a Mandelbrot/Julia escape-time grid.

Parameters:
extent: tuple = (-2.0, 1.0, -1.5, 1.5)#

(re_min, re_max, im_min, im_max) plotted region.

Type:

tuple

iterations: ndarray#

Escape iteration count at each pixel (equal to max_iter for points that never escaped).

Type:

ndarray, shape (ny, nx), int

max_iter: int = 100#

Iteration cap used.

Type:

int

class mathematicskit.fractals_chaos.GameOfLife(initial)[source]#

Bases: CellularAutomaton

Conway’s Game of Life: a 2D outer-totalistic cellular automaton (rule B3/S23).

A dead cell with exactly 3 live neighbors becomes alive (“birth”); a live cell with 2 or 3 live neighbors survives, otherwise it dies (“death” by isolation or overcrowding). Periodic (toroidal) boundary. The per-cell neighbor count uses a Numba-compiled kernel (_life_step()), following physicskit’s factory/dispatcher pattern for performance-critical stepping. See Gardner (1970), Scientific American 223.

Parameters:

initial (ndarray, shape (ny, nx)) – Initial grid (0s and 1s).

Examples

>>> import numpy as np
>>> # A "blinker": a 3-cell line, period-2 oscillator.
>>> grid = np.zeros((5, 5), dtype=np.int64)
>>> grid[2, 1:4] = 1
>>> life = GameOfLife(grid)
>>> life.step().tolist() == [[0]*5, [0,0,1,0,0], [0,0,1,0,0], [0,0,1,0,0], [0]*5]
True
>>> np.array_equal(life.step(), grid)
True
step()[source]#

Advance the automaton by one generation and return the new state.

Return type:

ndarray

class mathematicskit.fractals_chaos.IteratedFunctionSystem[source]#

Bases: ABC

Common base for random (chaos-game) iterated function systems.

Concrete subclasses set self.transforms (a list of (A, b, p) affine-map / probability triples, x -> A x + b chosen with probability p) in __init__.

abstractmethod generate(n_points, seed=0)[source]#

Run the chaos game for n_points iterations.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n_points, 2)

transforms: list = []#
class mathematicskit.fractals_chaos.SierpinskiCarpet[source]#

Bases: IteratedFunctionSystem

Sierpinski carpet via the chaos game: eight third-scale maps (no center).

The unit square is divided into a 3x3 grid of sub-squares, and every map but the center one is kept, each chosen with equal probability \(1/8\). The resulting point cloud has box-counting dimension \(\log 8/\log 3 \approx 1.893\). See Barnsley, Fractals Everywhere, 2nd ed., Ch. 3.

Examples

>>> carpet = SierpinskiCarpet()
>>> points = carpet.generate(2000, seed=0)
>>> points.shape
(2000, 2)
generate(n_points, seed=0)[source]#

Run the chaos game for n_points iterations.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n_points, 2)

class mathematicskit.fractals_chaos.SierpinskiTriangle(vertices=None)[source]#

Bases: IteratedFunctionSystem

Sierpinski triangle via the chaos game: three vertex-halving maps.

Each map is \(x \mapsto \tfrac12(x + v_k)\) for a triangle vertex \(v_k\), chosen with equal probability \(1/3\). The resulting point cloud has box-counting dimension \(\log 3/\log 2 \approx 1.585\) (see mathematicskit.fractals_chaos.systems.box_counting.box_counting_dimension()). See Barnsley, Fractals Everywhere, 2nd ed., Ch. 3.

Parameters:

vertices (ndarray, shape (3, 2), optional) – Triangle vertices; defaults to a unit equilateral triangle.

Examples

>>> tri = SierpinskiTriangle()
>>> points = tri.generate(2000, seed=0)
>>> points.shape
(2000, 2)
generate(n_points, seed=0)[source]#

Run the chaos game for n_points iterations.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n_points, 2)

mathematicskit.fractals_chaos.box_counting_dimension(points, box_sizes=None, n_sizes=12)[source]#

Estimate a 2D point set’s box-counting (Minkowski-Bouligand) dimension.

Covers the bounding box of points with a grid of side size, at each of several sizes, and counts \(N(\text{size})\) = the number of grid cells containing at least one point. \(N(\text{size}) \sim \text{size}^{-D}\) defines the box-counting dimension \(D\). Taking logs gives \(\log N = D\log(1/\text{size}) + c\), so \(D\) is the slope of \(\log N\) against \(\log(1/\text{size})\) (least-squares fit via numpy.polyfit()) – inverting the size axis is what turns the minus sign in the exponent into a positive slope. See Falconer, Fractal Geometry, 3rd ed., Ch. 3, equation (3.1).

Parameters:
  • points (ndarray) – A point cloud sampled from (or lying on) the set whose dimension is being estimated, e.g. from generate().

  • box_sizes (ndarray | None) – Box edge lengths to use; defaults to n_sizes sizes geometrically spaced between half the point set’s extent and 3% of it. The largest scale is excluded (a box near the full extent is prone to an off-by-one edge artifact from points sitting exactly at the boundary) and the smallest is capped at 3% (below which a finite point sample saturates – every box holds at most one point regardless of the set’s true dimension – biasing the estimate).

  • n_sizes (int) – Number of box sizes to use when box_sizes is not given.

Return type:

BoxCountingResult

Returns:

BoxCountingResult

Examples

>>> import numpy as np
>>> # A dense fill of the unit square has dimension ~2.
>>> rng = np.random.default_rng(0)
>>> square = rng.uniform(0.0, 1.0, size=(20000, 2))
>>> result = box_counting_dimension(square)
>>> 1.8 < result.dimension < 2.1
True
mathematicskit.fractals_chaos.chaos_game(transforms, n_points, seed=0, x0=None)[source]#

Run the “chaos game”: repeatedly apply a randomly chosen affine map.

At each step, one of the (A_k, b_k, p_k) triples is chosen with probability \(p_k\) and the current point \(x\) is replaced by \(A_k x + b_k\). For a contractive iterated function system (every \(A_k\) a contraction), the orbit converges (in distribution) onto the IFS’s unique attractor regardless of the starting point – the standard randomized algorithm for rendering self-similar fractals such as the Sierpinski triangle and the Barnsley fern. See Barnsley, Fractals Everywhere, 2nd ed., Ch. 3.

Parameters:
  • transforms (list[tuple[ndarray, ndarray, float]]) – (A, b, p) triples: a 2x2 matrix, a length-2 translation, and a selection probability (must sum to 1).

  • n_points (int) – Number of points to generate.

  • seed (int) – Random seed.

  • x0 (ndarray | None) – Starting point; defaults to the origin.

Return type:

ndarray

Returns:

ndarray, shape (n_points, 2)

Examples

>>> import numpy as np
>>> # A single constant map (A=0) with probability 1 sends every point to b.
>>> transforms = [(np.zeros((2, 2)), np.array([1.0, 2.0]), 1.0)]
>>> pts = chaos_game(transforms, n_points=5, x0=np.zeros(2))
>>> np.allclose(pts, [[1.0, 2.0]] * 5)
True
mathematicskit.fractals_chaos.divider_length(curve, ruler)[source]#

Length of a polygonal curve measured by walking a pair of dividers of opening ruler along it.

Starting at the first vertex, each step jumps to the first later point of the curve at distance ruler. Lewis Fry Richardson found empirically that measured coastline lengths grow as \(L(\varepsilon) \propto \varepsilon^{1-D}\) as the ruler \(\varepsilon\) shrinks, and Benoit Mandelbrot identified \(D\) as a fractal dimension.

Parameters:
  • curve (ndarray) – Vertices of a polyline, finely sampled compared with ruler.

  • ruler (float)

Return type:

float

Returns:

float – Number of full steps times ruler, plus the final partial step.

Examples

>>> line = np.column_stack([np.linspace(0, 1, 1001), np.zeros(1001)])
>>> round(divider_length(line, 0.1), 10)
1.0
mathematicskit.fractals_chaos.dla_cluster(n_particles, seed=0, size=801)[source]#

Grow a two-dimensional DLA cluster on the square lattice.

Particles are released one at a time on a circle just outside the cluster and random-walk until they touch it, then stick. Walkers that stray too far are relaunched. The branched clusters that result have fractal dimension about 1.71.

Parameters:
  • n_particles (int) – Target cluster size, including the seed particle.

  • seed (int) – Random seed.

  • size (int) – Lattice width; growth stops early if the cluster nears the edge.

Return type:

ndarray

Returns:

ndarray, shape (m, 2), int – Lattice coordinates of the cluster’s particles relative to the seed, in the order they stuck (m <= n_particles).

Examples

>>> cluster = dla_cluster(200, seed=1)
>>> cluster.shape
(200, 2)
>>> cluster[0].tolist()
[0, 0]
mathematicskit.fractals_chaos.henon_map(n_points, a=1.4, b=0.3, x0=(0.0, 0.0), discard=100)[source]#

An orbit of the Hénon map \((x, y) \mapsto (1 - a x^2 + y,\; b x)\).

With Michel Hénon’s parameters \(a = 1.4\), \(b = 0.3\) the orbit settles onto a strange attractor with box-counting dimension about 1.26. The map contracts areas by the factor \(|b|\) at every step.

Parameters:
  • n_points (int) – Number of orbit points returned.

  • a (float)

  • b (float)

  • x0 (tuple of float)

  • discard (int) – Initial transient iterations dropped.

Return type:

ndarray

Returns:

ndarray, shape (n_points, 2)

Examples

>>> orbit = henon_map(1000)
>>> bool(np.all(np.abs(orbit[:, 0]) < 1.5))
True
mathematicskit.fractals_chaos.hilbert_curve(order)[source]#

Vertices of the order-order Hilbert curve, visiting every cell of a \(2^n \times 2^n\) grid.

David Hilbert’s 1891 space-filling curve is the limit of these polygons. Consecutive vertices are grid neighbours, so points close along the curve stay close in the plane. Uses the standard index-to-coordinate conversion (Hacker’s Delight, Sec. 16-2).

Parameters:

order (int)

Return type:

ndarray

Returns:

ndarray, shape (4**order, 2), int – Integer grid coordinates in [0, 2**order).

Examples

>>> hilbert_curve(1).tolist()
[[0, 0], [0, 1], [1, 1], [1, 0]]
mathematicskit.fractals_chaos.julia_set(c, extent=(-2.0, 2.0, -2.0, 2.0), resolution=400, max_iter=100)[source]#

Escape-time grid for the filled Julia set of \(z \mapsto z^2 + c\), c fixed.

Unlike the Mandelbrot set (which sweeps c with a fixed \(z_0=0\)), the Julia set fixes c and sweeps the initial condition \(z_0\) over the grid. See Devaney, A First Course in Chaotic Dynamical Systems, 2nd ed., Ch. 15.

Parameters:
  • c (complex) – The fixed parameter.

  • extent (tuple of float) – (re_min, re_max, im_min, im_max) region of the complex plane.

  • resolution (int) – Number of pixels along the real axis.

  • max_iter (int) – Iteration cap.

Return type:

EscapeTimeResult

Returns:

EscapeTimeResult

Examples

>>> result = julia_set(c=-0.4 + 0.6j, resolution=10, max_iter=50)
>>> result.iterations.shape[1]
10
>>> # For c=0, z_{n+1}=z_n^2: |z0|<1 always stays bounded.
>>> from mathematicskit.fractals_chaos.systems.mandelbrot_julia import _escape_time_julia
>>> import numpy as np
>>> int(_escape_time_julia(np.array([[0.5]]), np.array([[0.0]]), 0.0, 0.0, 200)[0, 0])
200
mathematicskit.fractals_chaos.koch_curve(order, start=(0.0, 0.0), end=(1.0, 0.0))[source]#

Vertices of the Koch curve after order refinements of a segment.

Each step replaces every segment by four segments one third as long, with the middle third raised into an equilateral bump. Helge von Koch’s 1904 construction is continuous, nowhere differentiable, and has similarity dimension \(\log 4/\log 3 \approx 1.2619\).

Parameters:
  • order (int)

  • start (tuple of float) – Endpoints of the initial segment.

  • end (tuple of float) – Endpoints of the initial segment.

Return type:

ndarray

Returns:

ndarray, shape (4**order + 1, 2)

Examples

>>> koch_curve(2).shape
(17, 2)
mathematicskit.fractals_chaos.koch_snowflake(order)[source]#

Vertices of the closed Koch snowflake: three Koch curves on an equilateral triangle.

Parameters:

order (int)

Return type:

ndarray

Returns:

ndarray, shape (3 * 4**order + 1, 2) – The first vertex is repeated at the end to close the curve.

Examples

>>> koch_snowflake(1).shape
(13, 2)
mathematicskit.fractals_chaos.lsystem(axiom, rules, iterations)[source]#

Rewrite every symbol of axiom in parallel, iterations times.

Symbols without a rule are copied unchanged.

Parameters:
  • axiom (str)

  • rules (dict) – {symbol: replacement}.

  • iterations (int)

Return type:

str

Returns:

str

Examples

>>> lsystem("A", {"A": "AB", "B": "A"}, 5)  # Lindenmayer's algae: lengths are Fibonacci numbers
'ABAABABAABAAB'
mathematicskit.fractals_chaos.lyapunov_exponent_1d_map(f, fprime, x0, n_transient=500, n_iterations=2000)[source]#

Lyapunov exponent of a 1D map, \(\lambda = \lim_{n\to\infty} \tfrac1n \sum_{k=0}^{n-1} \ln|f'(x_k)|\).

Measures the average exponential rate of separation of nearby orbits: \(\lambda > 0\) signals sensitive dependence on initial conditions (chaos), \(\lambda < 0\) a stable periodic orbit or fixed point, \(\lambda = 0\) marks a bifurcation point. See Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 10.5, eq. (10.5.1)-(10.5.2).

Parameters:
  • f (Callable[[float], float]) – The map and its derivative.

  • fprime (Callable[[float], float]) – The map and its derivative.

  • x0 (float) – Initial condition.

  • n_transient (int) – Iterations discarded before averaging.

  • n_iterations (int) – Iterations averaged over.

Return type:

float

Returns:

float

Examples

>>> # Logistic map at r=4 is exactly conjugate to a tent map with
>>> # known Lyapunov exponent ln(2) ~ 0.693.
>>> f = lambda x: 4.0 * x * (1.0 - x)
>>> fprime = lambda x: 4.0 - 8.0 * x
>>> lam = lyapunov_exponent_1d_map(f, fprime, x0=0.4, n_iterations=100000)
>>> round(float(lam), 2)
0.69
mathematicskit.fractals_chaos.lyapunov_exponent_flow(rhs, x0, dt=0.01, n_steps=20000, renorm_every=10)[source]#

Largest Lyapunov exponent of a flow, via a shadow-trajectory method.

Integrates the system alongside a nearby “shadow” trajectory, periodically measuring and renormalizing their separation (Benettin et al.’s algorithm) – this avoids the naive approach of just watching two nearby trajectories diverge, which overflows once they separate enough to leave the linear regime. Uses simple forward-Euler stepping (adequate for an exponent estimate, where the log-average washes out an O(dt) local error) rather than mathematicskit.integrators, since renormalization must happen between steps at a resolution finer than any fixed-step RK4/Dormand-Prince call would expose. See Wolf et al. (1985), Physica D 16, and Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., Ch. 9.3.

Parameters:
  • rhs (Callable[[ndarray], ndarray]) – rhs(state) -> dstate/dt (autonomous).

  • x0 (ndarray) – Initial condition.

  • dt (float) – Integration step.

  • n_steps (int) – Total steps.

  • renorm_every (int) – Steps between renormalizations of the shadow separation.

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> # A simple contracting linear flow has a negative exponent.
>>> rhs = lambda x: -2.0 * x
>>> lam = lyapunov_exponent_flow(rhs, np.array([1.0]), n_steps=5000)
>>> bool(lam < 0)
True
mathematicskit.fractals_chaos.mandelbrot_set(extent=(-2.0, 1.0, -1.5, 1.5), resolution=400, max_iter=100)[source]#

Escape-time grid for the Mandelbrot set: \(z_0 = 0\), c sweeps the grid.

The Mandelbrot set is \(\{c \in \mathbb{C} : z_{n+1}=z_n^2+c,\ z_0=0 \text{ stays bounded}\}\). See Devaney, A First Course in Chaotic Dynamical Systems, 2nd ed., Ch. 15.

Parameters:
  • extent (tuple of float) – (re_min, re_max, im_min, im_max) region of the complex plane.

  • resolution (int) – Number of pixels along the real axis (the imaginary-axis count is scaled to keep pixels roughly square).

  • max_iter (int) – Iteration cap.

Return type:

EscapeTimeResult

Returns:

EscapeTimeResult

Examples

>>> result = mandelbrot_set(resolution=10, max_iter=50)
>>> result.iterations.shape[1]
10
>>> # c = -1 is in the Mandelbrot set (a period-2 cycle -1, 0, -1, ...).
>>> from mathematicskit.fractals_chaos.systems.mandelbrot_julia import _escape_time_mandelbrot
>>> import numpy as np
>>> int(_escape_time_mandelbrot(np.array([[-1.0]]), np.array([[0.0]]), 200)[0, 0])
200
mathematicskit.fractals_chaos.similarity_dimension(ratios)[source]#

The similarity dimension \(D\) solving Moran’s equation \(\sum_i r_i^D = 1\).

For a self-similar set built from maps with contraction ratios \(r_i\) that do not overlap too much (the open set condition), Patrick Moran proved in 1946 that \(D\) equals the Hausdorff dimension. With \(m\) equal ratios \(r\) it reduces to \(\log m / \log(1/r)\).

Parameters:

ratios (sequence of float) – Contraction ratios, each in \((0, 1)\).

Return type:

float

Returns:

float

Examples

>>> round(similarity_dimension([1 / 3] * 4), 6)  # Koch curve: log 4 / log 3
1.26186
>>> round(similarity_dimension([0.5, 0.25, 0.25]), 6)  # 0.5^D + 2 * 0.25^D = 1
1.0
mathematicskit.fractals_chaos.turtle_path(commands, angle, step=1.0, heading=90.0, draw='FG')[source]#

Interpret an L-system word as turtle-graphics moves.

F and G (or the symbols in draw) move forward drawing a line, f moves without drawing, +/- turn left/right by angle degrees, and [/] push and pop the turtle’s state (for branching plants). Other symbols are ignored.

Parameters:
  • commands (str)

  • angle (float) – Turning angle in degrees.

  • step (float)

  • heading (float) – Initial heading in degrees (90 = up).

  • draw (str) – Symbols that draw a forward line.

Return type:

list

Returns:

list of ndarray – Polylines, each of shape (k, 2); a new polyline starts after each pen-up move or state pop.

Examples

>>> [line.round(6).tolist() for line in turtle_path("F+F", 90, heading=0)]
[[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]]
mathematicskit.fractals_chaos.weierstrass_function(x, a=0.5, b=7.0, n_terms=30)[source]#

Partial sum of the Weierstrass function \(W(x) = \sum_{n \ge 0} a^n \cos(b^n \pi x)\).

Karl Weierstrass showed in 1872 that for \(0 < a < 1\), odd integer \(b\), and \(ab > 1 + 3\pi/2\) the sum is continuous everywhere but differentiable nowhere. G. H. Hardy (1916) extended this to all \(ab \ge 1\). Its graph has box-counting dimension \(2 + \log a / \log b\).

Parameters:
  • x (array_like)

  • a (float) – Amplitude ratio, \(0 < a < 1\).

  • b (float) – Frequency ratio.

  • n_terms (int) – Number of terms kept.

Return type:

ndarray

Returns:

ndarray

Examples

>>> round(float(weierstrass_function(0.0, a=0.5, n_terms=40)), 6)  # sum of 0.5^n
2.0