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:
IteratedFunctionSystemBarnsley’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
- class mathematicskit.fractals_chaos.BoxCountingResult(box_sizes, box_counts, dimension=0.0, extra=<factory>)[source]#
Bases:
objectContainer for a box-counting fractal-dimension estimate.
- class mathematicskit.fractals_chaos.CellularAutomaton[source]#
Bases:
ABCCommon base for cellular automata (1D elementary or 2D Game of Life).
- class mathematicskit.fractals_chaos.ElementaryCA(rule, width, initial=None)[source]#
Bases:
CellularAutomaton1D 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:
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]
- class mathematicskit.fractals_chaos.EscapeTimeResult(iterations, extent=(-2.0, 1.0, -1.5, 1.5), max_iter=100)[source]#
Bases:
objectContainer for a Mandelbrot/Julia escape-time grid.
- class mathematicskit.fractals_chaos.GameOfLife(initial)[source]#
Bases:
CellularAutomatonConway’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
- class mathematicskit.fractals_chaos.IteratedFunctionSystem[source]#
Bases:
ABCCommon 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 + bchosen with probabilityp) in__init__.
- class mathematicskit.fractals_chaos.SierpinskiCarpet[source]#
Bases:
IteratedFunctionSystemSierpinski 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)
- class mathematicskit.fractals_chaos.SierpinskiTriangle(vertices=None)[source]#
Bases:
IteratedFunctionSystemSierpinski 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)
- 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. fromgenerate().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:
- 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:
- Return type:
- 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
ruleralong 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:
- Return type:
- 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:
- Return type:
- 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:
- Return type:
- 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-
orderHilbert 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:
- 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\),
cfixed.Unlike the Mandelbrot set (which sweeps
cwith a fixed \(z_0=0\)), the Julia set fixescand sweeps the initial condition \(z_0\) over the grid. See Devaney, A First Course in Chaotic Dynamical Systems, 2nd ed., Ch. 15.- Parameters:
- Return type:
- 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
orderrefinements 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:
- Return type:
- 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:
- 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
axiomin parallel,iterationstimes.Symbols without a rule are copied unchanged.
- Parameters:
- Return type:
- 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:
- Return type:
- 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:
- Return type:
- 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\),
csweeps 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:
- Return type:
- 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:
- 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.
FandG(or the symbols indraw) move forward drawing a line,fmoves without drawing,+/-turn left/right byangledegrees, and[/]push and pop the turtle’s state (for branching plants). Other symbols are ignored.- Parameters:
- Return type:
- 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:
- Return type:
- Returns:
ndarray
Examples
>>> round(float(weierstrass_function(0.0, a=0.5, n_terms=40)), 6) # sum of 0.5^n 2.0