DEV Community: Rat Wolf (ratwolf)

Hopalong Map – fast, parallel, low-memory – two-pass checkpointed orbit decomposition

Rat Wolf (ratwolf) — Wed, 28 Oct 2026 12:32:25 +0000

Draws the Hopalong map (Barry Martin's map) from hundreds of millions or even billions of points, using all CPU cores, with memory that does not grow with the number of points.

Exact: the picture is bit-for-bit identical to a plain, single-threaded implementation, and identical for any number of threads (checked with 1 and 8).
Fast: all cores are used, and each core is kept busier than in a naive loop.
Small: memory stays fixed however many points you compute (about 76 MB on an 8-thread machine, on top of Python and its libraries).
Two output flavours: counts (hits per pixel) or a proper probability density on the colour bar.
Three colour scales: linear, square root, logarithmic.
Interactive window (matplotlib) plus PNG export.

Github Repository

Hopalong Map – fast, parallel, low-memory – two-pass checkpointed orbit decomposition
- Contents
- Quick start
- The Hopalong map
- The general approach (checkpoints and segments)
- The problem
- The idea in one picture
- Checkpoints (flags) and segments
- Why the result is exactly the same
- What is stored, and how much
- Lanes: keeping each core busy
- Where the time goes
- Counts vs. density
- Counts (default)
- Density
- What the density option is not
- Colour scales
- Using it from Python
- Configuration
- Code map
- Verification and benchmarks
- Limitations and notes

Quick start

Requirements: Python 3 (developed with 3.15), numpy, numba, matplotlib. Works on macOS and Linux (the memory read-out uses the Unix-only resource module; on Windows remove display_resource_usage).

pip install numpy numba matplotlib
python hopalong_fast.py

The script asks for six things:

Prompt	Meaning
`a`, `b`, `c`	The three parameters of the map (any real numbers; `a = 0` together with `c = 0` is rejected because the orbit would never leave the origin).
`n`	Number of points to compute, at least 1000. Scientific notation such as `1e8` is accepted.
Colour scale	`linear` (default), `sqrt` or `log` – see Colour scales.
Colour bar shows	`counts` (default) or `density` – see Counts vs. density.

Example session (times and memory will differ on your machine; this one was run on an Apple MacMini M2):

Enter a float value for "a": -3.14
Enter a float value for "b": 0.3
Enter a float value for "c": 0.3
Enter a positive integer value > 1000 for "n": 1e8
Colour scale (linear / sqrt / log) [linear]: sqrt
Colour bar shows (counts / density) [counts]: density
Pass 1 (extents + checkpoints, serial): 0.81 s
Pass 2 (parallel plot):                 0.08 s
Wall-clock time (compute): 0.89 s
Memory (RAM): 221.25 MByte used (peak)
Check: sum(density) * cell area = 1.000000
Image saved to: <OUTPUT_DIR>/py_num=100000000_a=-3.14_b=0.3_c=0.3_sqrt_density.png

Then a window opens with the picture, and a PNG is written to OUTPUT_DIR. Change OUTPUT_DIR at the top of the script to a folder that exists on your machine. The PNG is the raw 1000 × 1000 picture without axes; the file name encodes the parameters, the colour scale (_sqrt, _log) and the mode (_density).

The first start takes a few seconds because numba compiles the kernels; they are cached on disk afterwards.

The "Memory" line is the peak of the whole Python process, including Python, numba and matplotlib themselves (roughly 150 MB before anything is computed).

The Hopalong map

Every point is computed from the previous one. Starting at the origin, (x, y) = (0, 0):

where sign(x) is +1 for positive and -1 for negative x. The three numbers a, b, c decide what the picture looks like; for example a = -3.14, b = 0.3, c = 0.3 gives a ring-and-petal pattern.

Two facts about this map drive the whole design:

It is strictly serial. Point number 1,000,000 cannot be computed without first computing points 1 to 999,999. You cannot simply hand different parts of the orbit to different CPU cores.
It is deterministic. The same (x, y) always leads to the same next point. There is no randomness. If you know where you are, you know everything that follows.

A note on "map" versus "attractor". The Hopalong map is often called an attractor, but it does not pull orbits in from elsewhere. Away from the lines x = 0 and b * x = c, its Jacobian determinant (the factor by which the map stretches areas) is exactly 1, so the map preserves area. You can check this by hand: the new x changes with y at slope 1, and the new y changes with x at slope -1 and does not depend on y, so the determinant is 1. Orbits therefore stay in regions that are determined by where they start, and different start points generally give different pictures. The picture is the orbit from (0, 0), not a property of a, b and c alone. That is why the exact orbit is replayed from checkpoints, instead of starting the segments at independent random points.

The general approach (checkpoints and segments)

The problem

To draw the picture you lay a grid over the plane (here 1000 × 1000 cells) and, for every point of the orbit, add one tally mark to the cell it lands in. Cells that were visited often end up bright. Such a grid of tally counts is called a histogram.

Two things make this harder than it sounds:

You must know the size of the picture before you can start tallying. Where is the leftmost point, the topmost? That is only known after the whole orbit has been walked. This is why every implementation needs two passes: one to measure, one to draw.
Storing the orbit is out of the question. One point is two decimal numbers, 16 bytes. A billion points would be 16 GB.

A plain implementation therefore walks the orbit twice, one after the other, on a single core.

The idea in one picture

Think of the orbit as a very long hiking trail where each step depends on the previous one, so one hiker has to walk it in order. The trick:

Pass 1: one hiker walks the whole trail once. On the way they note the extreme positions (how big must the map be?) and plant a flag every so often. A flag records just the exact position at that spot.
Pass 2: now many hikers start at the same time, each at a different flag, and each walks only the stretch up to the next flag. Because the trail is deterministic, they all take exactly the steps the single hiker would have taken. Each hiker makes tally marks on their own private tally sheet, so nobody gets in anybody's way. At the end, the tally sheets are added up.

 PASS 1  -  one worker walks the whole orbit once (serial)
   |         remembers the smallest and largest x and y   -> how big is the picture?
   |         plants a flag at the start of every segment  -> a few bytes each
   v
 PASS 2  -  many workers, all at the same time (parallel)

   worker 1:  walks segments 0-3    -> marks its own tally sheet 1 --+
   worker 2:  walks segments 4-7    -> marks its own tally sheet 2 --+
   worker 3:  walks segments 8-11   -> marks its own tally sheet 3 --+--> add up all
   ...                                                               |    tally sheets
   worker k:  walks segments ...    -> marks its own tally sheet k --+        |
                                                                              v
                                                                       final picture

Checkpoints (flags) and segments

The orbit is cut into consecutive pieces called segments. A checkpoint is the exact position (x, y) at the start of a segment.

   flag 0        flag 1        flag 2        flag 3        flag 4
     |             |             |             |             |
     v             v             v             v             v
     +-------------+-------------+-------------+-------------+--------- ...
     |  segment 0  |  segment 1  |  segment 2  |  segment 3  |  segment 4
     +-------------+-------------+-------------+-------------+--------- ...
   start (0,0)                                                    the orbit goes on

A worker that starts at flag k just keeps computing until it reaches the position stored in flag k+1, adding a tally mark for every point on the way.
If the number of points does not divide evenly into segments, the leftover points are simply added to the last segment, so every point is visited exactly once.

Why the result is exactly the same

Starting from an exact checkpoint, a worker repeats exactly the same arithmetic as the single hiker did, so it produces exactly the same numbers. The set of tally marks is the same as in a single sequential run, only the order in which they are made differs, and adding up whole numbers does not depend on the order. So the picture is identical, no matter how many cores are used.

This matters more than it may seem: the map is chaotic, so even a difference in the 16th decimal place would grow until the replayed orbit looked completely different. Both passes therefore use the very same step function and plain IEEE floating-point arithmetic (do not switch on numba's fastmath).

What is stored, and how much

What	Size	Depends on `n`?
The orbit itself	not stored at all	–
Checkpoints	16 bytes per segment, about 1 KB in total	no
One tally sheet per task	4 MB each (1000 × 1000 counters of 4 bytes)	no
Final picture	8 MB	no

With T CPU threads the script uses 2 × T tasks and 4 segments per task, so 8 × T segments. On an 8-thread machine that is 64 segments, i.e. 1 KB of checkpoints, instead of 16 GB for a billion stored points. The only memory that grows is the number of tally sheets (one per task), and it grows with the number of cores and with the picture resolution, never with n.

The counters of a tally sheet are 4-byte numbers, which overflow at about 4.3 billion. To be safe, no single task is allowed to plot more than 2 billion points (MAX_POINTS_PER_TASK); for very large n the script simply uses more tasks. The sheets are then added up in 8-byte numbers.

Lanes: keeping each core busy

Each step needs the result of the previous step, including a square root, which takes the CPU a moment. A core working on one chain of steps spends part of its time waiting. So each task works on four consecutive segments in lock-step ("lanes"): while the CPU waits for the square root of lane 1, it already works on lanes 2, 3 and 4. This uses no extra cores, it just makes each core busier. In the benchmark below this makes pass 2 about 3.7 times faster per core (0.82 s down to 0.22 s on an Apple M2; the gain depends on the CPU).

Where the time goes

Pass 1 is serial: it costs roughly as much as the walk itself.
Pass 2 is parallel: its time shrinks with the number of cores.

So the total time is roughly pass 1 + pass 2 / number of cores. Beyond a handful of cores, pass 1 dominates, and extra cores no longer help much. The script prints both times so you can see this. On an 8-thread Apple M2, pass 1 takes 0.81 s of the 0.87 s total (see the benchmark below).

Counts vs. density

Both options run the same two passes and show practically the same picture. They differ in how points are assigned to cells and in what the colour bar says.

Counts (default)

The colour bar shows how many points landed in each pixel.

A point is assigned to the nearest grid point, where the grid points lie exactly on the minimum and maximum of the orbit. With W columns:

(rows work the same way with y and H). This reproduces the output of the plain single-threaded two-pass implementation exactly. It has one property that makes it unsuitable for a density: the first and last row and column only collect points from half a pixel, so they are only half as wide as all other cells.

Density

The colour bar shows the probability density: how concentrated the orbit's points are per unit area of the x-y plane. For this the grid must consist of equal-width cells, exactly like numpy.histogram2d:

(a point exactly on the right edge goes into the last cell, as in numpy). With the cell sizes

the density in cell (i, j), with n_ij hits and N points in total, is

Added up over the whole area, the density gives exactly 1:

The script prints this as Check: sum(density) * cell area = 1.000000.

Two points are worth stressing:

The density is just the counts times one constant, 1 / (N * dx * dy). The counts stay exact whole numbers; nothing is rounded, smoothed or re-binned when converting. So the picture is the same and only the numbers on the colour bar change. (N equals the n you typed, because the extents were measured on the very same orbit and therefore every point falls inside the grid.)
It is exactly what numpy.histogram2d(x, y, bins=..., density=True) would give for the same points and the same range, only without ever holding the points in memory. Counting into bins is also all histogram2d does; the binning loses exactly the same information in both cases.

What the density option is not

Not the exact density of the map. It is a histogram estimate on a 1000 × 1000 grid. For many parameter sets the points pile up on thin curves and rings; for such patterns the density value in a cell depends on how large the cell is, and a finer grid gives higher values on the curves. There may be no resolution-independent "true" value to converge to. A larger n reduces the noise (cells with few hits are noisy), but it does not remove the dependence on the cell size.
Not a probability per pixel. It is a probability per unit area. The probability that a point falls into one cell is n_ij / N, which is the density multiplied by the cell area:

Not independent random samples. The points come from one deterministic orbit and are strongly correlated. The numbers describe how often this orbit, started at (0, 0), visits each cell; since the map preserves area (see The Hopalong map), other start points give different pictures. Textbook error estimates for independent samples (such as the square root of the count) are only a rough guide here.
Not smoothed. There is no kernel density estimate and no blurring; every cell holds the raw tally.
Not normalised to a maximum of 1. The numbers on the colour bar are real density values. Only the colours are scaled relative to the busiest cell.
Not affected by the colour scale. sqrt and log change which colour a value gets, never the numbers on the colour bar.
Not stored in the PNG. The PNG contains colours, not data. To get the density values themselves, use the Python interface below.

Colour scales

A pixel's relative brightness (0 to 1) is derived from its hit count n and the highest hit count n_max in the picture. The scale only decides how counts are mapped to colours; the colour bar ticks always show the real counts (or densities).

linear – proportional to the count. Rarely visited regions end up (almost) black:

sqrt – lifts sparse regions:

log – strongest lift, even single hits become visible:

The tick labels are spaced evenly along the colour bar, not evenly in count value, so with log you see labels like 0, 1, 3, 8, 17, 36, ...

Using it from Python

The file can be imported; nothing runs on import.

import numpy as np
import hopalong_fast as hf

hf.warm_up()  # compile, or load the cached kernels, once

a, b, c, n = -3.14, 0.3, 0.3, 100_000_000

# uniform_bins=True -> equal-width cells, as needed for a density
counts, extents, (t_pass1, t_pass2) = hf.compute_trajectory_image(
    a, b, c, n, (1000, 1000), uniform_bins=True
)

density = counts * hf.density_factor(counts, extents)   # integrates to 1 over the area
np.save("density.npy", density)

counts is an unsigned 64-bit array of shape (rows, columns) = (y, x). Row 0 is the lowest y value. numpy.histogram2d returns shape (x, y), so compare with density.T.
extents is (x_min, x_max, y_min, y_max). The cell edges are np.linspace(x_min, x_max, columns + 1) and np.linspace(y_min, y_max, rows + 1).
With uniform_bins=False (the default) you get the pixel-mode counts of the plain single-threaded two-pass implementation.

Configuration

Constants at the top of hopalong_fast.py:

Constant	Default	Meaning
`LANES`	`4`	Segments advanced in lock-step per task. The kernel is written out for exactly 4 lanes; do not change this number without rewriting `second_pass`.
`PARALLEL_THRESHOLD`	`1_000_000`	Below this many points a single task is used; the parallel set-up would cost more than it saves.
`MAX_POINTS_PER_TASK`	`1 << 31`	Upper bound on the points one task plots, so that its 4-byte counters cannot overflow.
`OUTPUT_DIR`	`/Users/ralf/Projects/hopalong_pictures`	Where PNG files are written. Change this.
`SCALES`	`linear`, `sqrt`, `log`	Allowed colour scales.
`OUTPUTS`	`counts`, `density`	Allowed colour-bar modes.

Other knobs:

Number of threads: set the environment variable NUMBA_NUM_THREADS, for example NUMBA_NUM_THREADS=4 python hopalong_fast.py.
Resolution and colour map: main(image_size=(1000, 1000), color_map='hot'). Every task needs rows × columns × 4 bytes, so memory grows with the number of tasks and with the resolution: at 4000 × 4000 each tally sheet is 64 MB.

Code map

Function	Role
`make_plan(n)`	Chooses the number of tasks and cuts `n` points into segments.
`_step`	One Hopalong iteration. The same function is used in both passes.
`first_pass`	Pass 1: walks the orbit serially; returns the extents and the checkpoints.
`_hit`	Maps a point to its cell (pixel mode or equal-width mode) and counts it. Points outside the grid, or NaN/infinite points, are discarded.
`second_pass`	Pass 2: parallel over tasks (`prange`), four lanes per task, one private tally sheet per task.
`compute_trajectory_image`	Runs both passes, sums the tally sheets, returns counts, extents and timings.
`density_factor`	The constant `1 / (N * dx * dy)` that turns counts into a density.
`make_norm`, `colorbar`, `render_trajectory_image`, `save_png`	Colour scales, colour bar, window and PNG output.

Verification and benchmarks

Measured on: Apple M2, 8 threads, 8 GB RAM, macOS (Darwin 27.2.0), Python 3.15.1, numpy 2.5.4, numba 0.68.0.

All numbers in this section are produced by bench_readme.py; run python bench_readme.py to reproduce them (add --big for the 100-million-point memory test).

Checks:

Counts mode vs. the plain single-threaded two-pass implementation: the extents and the complete 1000 × 1000 image are bit-for-bit identical for 4 different parameter sets, including values of n that do not divide evenly into segments and runs with several tasks.
Density mode vs. numpy.histogram2d(..., density=True): with 2 to 3 million points on 3 parameter sets, no cell differs in its count, the densities agree to floating-point rounding (relative difference about 2e-13), and the integral is within 5e-14 of 1.
Colour scales: the saved PNGs were compared with the three formulas above (using a grey colour map) and agree to within 2 grey levels out of 255, which is the rounding of a 256-entry colour map and 8-bit PNG; the orientation (y pointing up) was checked at the same time.

Speed, 100,000,000 points, parameters a = 0.4, b = 1.0, c = 0.0, 1000 × 1000 grid, median of 3 runs. Both versions were timed on this machine in the same session; the plain loop is single-threaded:

	Pass 1	Pass 2	Total
Plain two-pass numba loop (single-threaded -->"hopalong_basic.py")	0.80 s	0.82 s	1.62 s
`hopalong_fast.py`, 1 thread	0.81 s	0.22 s	1.02 s
`hopalong_fast.py`, 8 threads	0.81 s	0.06 s	0.87 s

Memory, density histogram of the same orbit (extra peak RAM on top of what the interpreter and the libraries already use; hopalong_fast.py with 8 threads = 16 tasks, so at most 16 × 4 MB of tally sheets plus the 8 MB result; the operating system only counts the parts that were actually touched):

	Extra peak RAM, 20,000,000 points	100,000,000 points
Generate the points into arrays, then `numpy.histogram2d`	1.1 GB	not run (estimated at about 6 GB)
`hopalong_fast.py`	76 MB	76 MB

Limitations and notes

Pass 1 is serial. That is the floor for the total time; see Where the time goes. It cannot simply be skipped by starting segments at independent random points, because the picture depends on the start point (see The Hopalong map).
First start compiles. numba compiles the kernels on first use (a few seconds) and caches them in __pycache__.
Not independent samples, and not a continuous density. See What the density option is not.
Extreme parameters. If the orbit overflows to infinity or NaN, the extents and the picture are meaningless.
Do not enable fastmath. Pass 1 and pass 2 must compute bit-identical orbits for the checkpoints to be valid (the map is chaotic, see Why the result is exactly the same).
Platform. The memory read-out uses the Unix-only resource module (macOS and Linux). Peak memory is reported in the unit appropriate for each platform.

One-pass Hopalong orbit density plot

Rat Wolf (ratwolf) — Sun, 25 Oct 2026 07:05:26 +0000

This script renders the classic Hopalong map as a visit-density image. It
streams the orbit into a fixed-size histogram and adapts the represented
coordinate range when necessary, rather than storing every orbit point or
making a full preliminary pass to discover the bounds.

Target

The goal is to plot a long orbit with:

Memory bounded by the histogram dimensions, rather than the number of iterations.
Exact integer visit counts at the histogram's final cell resolution.
Automatic coverage when the orbit moves beyond the current grid, at the finest cell size that still fits the whole orbit.
A useful density visualization, optionally Gaussian-smoothed.

The implementation uses NumPy arrays, Matplotlib, and Numba. Save the script
as hopalong_orbit_onepass.py or substitute its actual filename in the
commands below.

Implementation details

Orbit recurrence

Starting at (x, y) = (0, 0), each step uses the previous x:

old_x = x
x = y - sign(old_x) * sqrt(abs(b * old_x - c))
y = a - old_x

The default parameters are a=0.7, b=0.4, and c=0. The first --skip
steps are treated as a transient and are not added to the histogram.

Starting grid

The first 4,096 plotted points are held in a small buffer (64 KB). Their
bounding box selects the starting cell size: the smallest power of two times
H_MIN = 2**-32 at which the box fits the grid. The grid is centred on the
box, the buffered points are binned, and the orbit simply continues from where
it was. No recurrence step is computed twice, and there is no separate pilot
orbit.

This is a sizing heuristic only. Later points are still checked, and the grid
is reframed whenever a point falls outside it. The attractor's extent usually
keeps growing with the number of points, so the buffer cannot predict the final
size; it only avoids a burst of cheap-but-pointless regridding at the start.

Streaming histogram and dyadic growth

The histogram is a square size x size array of signed 64-bit integer
counts. A point inside the current bounds is mapped to a cell by scaling its
coordinates by the inverse cell width and converting to integer indices. The
hot loop does nothing else: it carries no bookkeeping of its own.

When a plotted point falls outside the grid, the grid is reframed:

The occupied block of cells is located with one pass over the array.
If the occupied block and the new point still fit at the current cell size, the counts are only translated, so the block is centred. Nothing is merged and no resolution is lost.
Otherwise the cell width is doubled: each 2x2 block of counts is merged into one cell by integer addition, and the merged block is centred in the grid. This repeats until everything fits.
The triggering point is binned and orbit processing resumes.

Centring leaves room on both sides, so a reframe is rarely needed again soon.
Only the occupied block is read, written and cleared, using two preallocated
arrays that swap roles, so a reframe costs a few milliseconds at 2048 x 2048.

All cell widths are H_MIN * 2**k, so old and new cell boundaries stay
aligned: translating and 2x2 merging are exact aggregations onto the new
lattice. Earlier points are not lost or approximated in count. However, their
original fine-cell locations cannot be recovered after merging: all points are
ultimately represented at the final, coarser resolution.

Because the grid doubles only when the occupied region no longer fits, the
final occupied region covers between about 50% and 100% of the grid side.

Output image

Empty outer rows and columns are cropped before display. By default the
script displays log1p(count) with a magma color map. With --smooth, it
applies a separable Gaussian filter with sigma 1.5 cells before taking the
logarithm. --save writes hopalong_orbit.png to the selected output
directory instead of opening a window.

Tuning resolution with `--bins`

Use --bins to set the side length of the square count grid. For example:

python hopalong_orbit_onepass.py --bins 4096 --save

Larger values generally produce a smaller final cell width in the status
output, meaning finer spatial resolution. The relationship is not perfectly
linear: the grid only doubles its cell width when the orbit outgrows it, so
the final cell width also depends on the orbit's coordinate range. Compare
cell values for the same parameters and iteration count; a smaller value
means finer bins.

The memory cost grows quadratically with --bins. The two 64-bit count grids
alone use about 16 * bins * bins bytes: approximately 64 MiB at 2048 and
256 MiB at 4096, before temporary arrays used for plotting or smoothing. The
option is rounded up to an even grid side internally because the grid growth
merges cells in 2x2 blocks.

Result

The result is a two-dimensional image of how often the orbit visits each
cell. Dense regions appear brighter, making the orbit's spatial structure
visible without keeping the full point sequence in memory.

On normal completion, the script reports the number of binned points, grid
size, final cell width, number of grid doublings, elapsed time, and occupied
image dimensions and coordinate extent. If the orbit reaches the escape
threshold, processing stops and the report indicates that it escaped; fewer
than the requested number of points may then be plotted. If it escapes before
any point is plotted, the script exits with a message.

Caveats

Resolution decreases on growth. Each doubling of the coordinate span also doubles the cell width. This preserves a fixed-size grid, but an orbit whose extent keeps growing with the number of points ends up with a coarser final image the longer it runs. The cell width is the smallest power of two at which the whole occupied region fits.
Counts are exact only for the final bins. Integer merges preserve visit totals, but do not preserve sub-cell positions or allow a later higher-resolution reconstruction.
Floating-point bin edges. Orbit coordinates and grid-to-cell mapping use floating-point arithmetic. Points extremely close to a cell boundary can be assigned to either neighboring cell due to rounding.
Reframing has a small fixed cost. Each reframe scans the grid once and moves only the occupied block. Typical runs need zero to a few dozen reframes, and heavy-tailed orbits somewhat more; the total is on the order of 0.1 s at 2048 x 2048 and does not scale with the number of iterations.
Finer images touch more memory. A sharper final grid has more occupied cells, so the scattered histogram updates miss the cache a little more often. Measured on one machine, this was a few percent, and it is the price of the better resolution rather than of the reframing code.
Memory is bounded by grid size, not necessarily small. Two count arrays are kept (about 16 * bins * bins bytes in total), plus temporary memory for plotting/smoothing.
Finest cell size. Cells never get finer than 2**-32. This only matters for orbits whose entire extent is below about 5e-7.
Escape guard. The orbit stops if either coordinate's absolute value exceeds 1e6. This is a safety cutoff, not a general mathematical test that every parameter choice produces a bounded orbit.
Smoothing changes the displayed density. Gaussian smoothing is a visualization step, not part of the visit counts. Its sigma is measured in cells, so the apparent blur depends on the final cell size. Zero padding can reduce smoothed values near the image boundary.

Use cases

Exploring the spatial density and structure of Hopalong orbits.
Quickly comparing plots produced with different a, b, and c parameters.
Producing a bounded-memory visualization for very long runs when a fixed-resolution image is acceptable.
Creating a smoothed, presentation-friendly density image with --smooth.

This method is intended for visualization, not for retaining a high
precision point cloud or supporting later zooming and re-binning. For those
needs, store orbit points or use a sparse/finer-resolution representation,
accepting the associated memory or processing cost.

Requirements

Install the dependencies in the Python environment used to run the script:

python -m pip install numpy matplotlib numba

Usage

python hopalong_orbit_onepass.py
python hopalong_orbit_onepass.py --smooth
python hopalong_orbit_onepass.py --save
python hopalong_orbit_onepass.py --a 1.1 --b 0.5 --c 0.3 --iters 2000000 --save

--bins sets the square grid side length. If omitted, the script derives it
from --iters and caps it at 2,048; the resulting side length is made even
for the 2x2 merge operation. Increase it for finer resolution at greater
memory cost; check the reported cell value to compare the actual final
resolution. --outdir selects where the PNG is saved.

Code and README: GitHub repository

The Hopalong "Attractor" That Isn't One

Rat Wolf (ratwolf) — Sun, 25 Oct 2026 06:53:31 +0000

A classic computer-art map has no attractor. Here is how to see its regular and chaotic regions with a few lines of Python and numba

Orbits of the Hopalong map. Cool colours are regular orbits circling islands, warm colours are chaotic orbits.

Search for "Hopalong attractor" and you get glowing, symmetric pictures made from nothing but two lines of arithmetic. They are fun to produce and easy to share. But the word "attractor" hides what the picture actually is.

This article walks through what the Hopalong map does, why the name is misleading, and how to see all of its phase space at once instead of one orbit at a time. None of this is new to dynamical-systems people. The point is that a short numba script makes it easy to see for yourself.

The map

Barry Martin's Hopalong map [2] turns each point of an orbit into the next one, using three parameters a, b and c:

Start at (0,0), repeat millions of times, plot every point. With a = 0.7, b = 0.4 and c = 0 you get a picture like the one below. (In German the map appeared as the program HÜPFER in A.K. Dewdney's column in Spektrum der Wissenschaft [3].)

One orbit, one picture

The core loop is short. With numba it compiles to machine code and runs ten million iterations in about a second:

import numpy as np
from numba import njit

@njit
def orbit(a, b, c, n):
    xs, ys = np.empty(n), np.empty(n)
    x, y = 0.0, 0.0
    for i in range(n):
        x, y = y - np.sign(x) * np.sqrt(abs(b * x - c)), a - x
        xs[i], ys[i] = x, y
    return xs, ys

(This is a simplified version of the script's loop, which also discards a short transient and stops if an orbit escapes to infinity.)

A scatter plot of ten million points is slow and saturates into a solid blob. A log-scaled 2D histogram is faster and shows structure:

Notice the black regions. The orbit never enters them, and they are the first hint that something is going on.

The picture also depends on how long you run the orbit. In one run it stayed inside the same box from 2 million up to 100 million points, then leaked into a larger region somewhere between 100 and 300 million points, and the box grew from about x −0.76 to 1.02 to about x −1.38 to 1.92. The exact point at which this happens depends on floating-point details, so treat it as an observation, not a constant. The script measures the extent over the full run so the image is framed correctly either way.

For long runs the script fills the 2D histogram while the orbit runs instead of storing the points, which keeps memory constant (in one test, 300 million points took about 11 seconds and about 200 MB) and gave counts identical to NumPy's histogram2d in the checks.

Why it is not an attractor

An attractor is a set that nearby starting points are pulled toward. That requires the dynamics to shrink areas in phase space. The Hopalong map does not. Wherever it is differentiable, its Jacobian is

and its determinant is

The map is also invertible: equation (2) gives the old x from the new y, and then equation (1) gives the old y. A map that is invertible and has determinant 1 preserves area. Nothing contracts, so there is nothing to attract.

What you see in the classic picture is the footprint of one orbit, and a different starting point can give a different picture. The name is a loose label, and it is common to call any shape drawn by an iterated map an attractor. It is fine to keep using it, as long as you know what the image is.

Two more facts about the map will matter later. It keeps its form under the scaling

for any scale factor s, so a and b behave like lengths and c like a length squared. And its fixed points satisfy

with y = a − x.

Seeing all of phase space

If one orbit only tells part of the story, we want a way to classify every starting point. The standard tool is the Lyapunov exponent: how fast two nearby orbits separate.

A regular orbit stays on a closed curve or a chain of islands. Nearby orbits separate slowly and the exponent is about zero.
A chaotic orbit fills an area. Nearby orbits separate exponentially and the exponent is positive.

To estimate it, a tangent vector is pushed forward by the Jacobian (3) and renormalised at every step; the hat marks a vector rescaled to length 1:

The largest finite-time exponent over N steps is the mean logarithmic stretching:

In code, it is the same loop with a few extra lines:

@njit
def lyapunov(a, b, c, x, y, n):
    vx, vy = 1.0, 0.0
    total = 0.0
    for i in range(n):
        d = b * x - c
        r = max(np.sqrt(abs(d)), 1e-12)       # guard the sqrt singularity
        dfdx = -np.sign(x) * b * np.sign(d) / (2.0 * r)
        x, y = y - np.sign(x) * np.sqrt(abs(d)), a - x
        vx, vy = dfdx * vx + vy, -vx
        norm = np.sqrt(vx * vx + vy * vy)
        vx, vy = vx / norm, vy / norm
        total += np.log(norm)
    return total / n

(Also simplified: the script additionally discards the first 500 steps so the tangent vector can align, and stops if the orbit escapes.)

Run it over a 400 × 400 grid of starting points (numba's prange spreads the grid over your cores) and colour each point by its exponent:

In the plotted window (x from −2.5 to 2.9, y from −2.2 to 3.2), about 70% of the starting points are regular and about 30% are chaotic, with 8,000 iterations per point. The dark blobs are islands of regular motion, and the bright bands between them are the chaotic region. The cross marks (0,0). Its orbit has an exponent of about 0.018 and lies in the chaotic region, which is why the classic picture has those empty holes: it is one chaotic orbit, and most of the empty holes are regular islands it does not enter (in a pixel-by-pixel check, about 81% of the empty pixels start regular orbits).

Starting points far from the centre, such as (50,50) or (1000,−1000), turned out to be regular, but a coarse survey still found chaotic orbits at every distance we probed (up to 64 units from the origin), so the structure has no sharp outer edge.

A prettier version

The classification map is informative but pixelated. For a more attractive picture, start a few hundred orbits along a line through phase space and draw each one as its own set of points. Regular orbits become nested rings around the islands, and chaotic ones fill the bands in between:

The warm colours (yellow to red) are chaotic orbits, brighter for a larger exponent. The cool colours (cyan to magenta) are regular orbits, but there they only encode the position of the starting point along the line, so they are decoration, not data.

The picture shows a ring of seven large islands around the centre and an inner ring of about ten smaller ones. The line of starting points only grazes the rims of the large islands, so their interiors stay empty. We did not analyse these island chains further.

What this is, and what it isn't

Mixed regular and chaotic motion is the standard behaviour of area-preserving maps, and this article reproduces it for one parameter set. The references at the end give the sources of the map; the results here have not been compared with them or with other literature on its dynamics. A few limits are worth knowing:

The percentages depend on the window. The chaotic share is about 31% in a window twice as wide as the one above, but about 20% at four times as wide. The numbers describe the plotted region, not the whole plane.
The threshold matters a little. The chaotic share stays at 30–31% for thresholds from 0.002 to 0.004 and for 4,000 to 20,000 iterations. At a threshold of 0.008 it drops to about 25–27%, because weakly chaotic orbits are no longer counted.
The exponent is finite-time. Orbits that stick near an island can look regular over a short run, so the chaotic share can be under-counted near island boundaries. Even the exponent of the orbit from (0,0) varies by about 10% between runs of different length.
The square-root term is singular at b·x = c. The derivative is unbounded there, so classical results for smooth maps are an analogy, not a guarantee.
One parameter set. Nothing here says how the picture changes as a, b and c vary.

Try it yourself

The full script produces all three views and the header image, with a, b, c and the iteration counts as command-line options:

pip install numpy matplotlib numba
python hopalong_phase_space.py --save
python hopalong_phase_space.py --a 1.1 --b 0.5 --c 0.3 --iters 2000000

The script runs in two passes. A quick first pass finds the length scale of the map from a, b and c, the centre of its structure (the map's fixed point) and the exact extent of the orbit from (0,0). The second pass draws the three views with that framing, so you can change a, b, c and the number of points and still get a framed picture. How much of the structure fits the window depends on the parameters, though: the window is 3.8 length scales around the fixed point (or around the middle of the orbit if there is none), which is a zoom choice, and for maps with much larger orbits you will want --window-scale.

Code and README: GitHub repository

A pretty picture is a good way into a real question. Here the question is how a map that preserves area still produces such rich structure, and the answer is that its phase space is not one thing but a patchwork of regular and chaotic regions. The natural next step is to ask how that patchwork changes as the parameters move.

References

[1] J. Lansdown and R. A. Earnshaw (eds.), Computers in Art, Design and Animation. New York: Springer-Verlag, 1989.

[2] Barry Martin, "Graphic Potential of Recursive Functions," in Computers in Art, Design and Animation [1], pp. 109–129.

[3] A.K. Dewdney, Psychotapeten, algorithm HÜPFER, in Spektrum der Wissenschaft: Computer Kurzweil. Heidelberg: Spektrum der Wissenschaft Verlagsgesellschaft mbH & Co., 1988. (German edition of Scientific American.)

Discreteness as the Late-Time Limit of a Smooth Flow

Rat Wolf (ratwolf) — Fri, 09 Oct 2026 03:55:32 +0000

An exactly solvable gradient flow whose time-evolution maps are smooth bijections of the real line at every finite time, and whose infinite-time limit is rounding to the nearest integer.

Abstract

We study the gradient flow of the periodic potential

Integers are its stable fixed points and half-integers its unstable ones. The flow can be solved in closed form. At every finite time it is a smooth, strictly increasing bijection of the real line; the family of maps forms a one-parameter group; and as time tends to infinity every open cell between consecutive half-integers collapses onto the integer it contains. "Integer-valued" therefore appears as the endpoint of a continuous process with a single, exactly composable resolution parameter, the time. The integers are not generated by the flow: they enter through the period of the potential. What the flow provides is an exactly solvable, exactly composable route to that prescribed discrete target. We state the closed form, derive the main properties, give a geometric interpretation as a hyperbolic Möbius dilation of the circle, and give the exact time it takes to resolve a point to a given accuracy.

The exact solution

The full derivations and proofs, the verification script and the visualization code are in the GitHub repo: Github

I Built a Finite-Element Solver With Nothing to Hide

Rat Wolf (ratwolf) — Thu, 16 Jul 2026 01:13:07 +0000

A transparent, dependency-free capacitor simulator — and an honest account of exactly where a simple mesh lets it down.

Every physics course teaches the parallel-plate capacitor formula in the first week:

Clean. Memorable. Wrong the moment your capacitor isn't two infinite plates in a vacuum. Real capacitors have edges. Real dielectrics are layered. Real geometries are curved. The formula survives exactly as long as the idealization does.

So you reach for simulation instead — and run into two unsatisfying options.

Option one: a commercial FEM package. Powerful, polished, and a complete black box. You get a number. You don't get to see the mesh, the matrix, or the fifty implicit decisions that produced it.

Option two: a serious open-source stack — FEniCS, deal.II, scikit-fem with gmsh underneath. Genuinely open, genuinely capable, and also genuinely heavy: compiled mesh generators, several abstraction layers, a real learning curve before you can even see where the physics happens.

I wanted a third thing. So I set myself three rules:

Every formula in the code has to be traceable to something you could write on a whiteboard.
Every accuracy claim has to be measured, not assumed.
Every known limitation has to be printed by the program itself, not buried in a comment nobody reads.

Here's what building to those rules actually looked like.

The Physics, Compressed

A capacitor problem is one equation. Given some conductors at fixed voltages and a dielectric — possibly several, possibly varying in space — filling the gaps between them, the electric potential $V$ satisfies:

That's the whole physics. Everything else — field, displacement, stored energy, capacitance — falls out of V:

That last formula is worth a pause. The obvious way to get capacitance out of a simulation is to integrate surface charge on a conductor. This solver doesn't do that — it computes the total energy stored in the field instead and rearranges that same relationship to solve for C. Charge integration means differentiating a numerically noisy field exactly at a boundary, the single worst place to differentiate anything. Energy integration only needs a field the solver has already computed everywhere. Small choice, but the kind that separates code that runs from code that's numerically honest.

What "No Black Box" Actually Means

Every finite element method runs on the same trick: chop the domain into small pieces, assume the solution is simple on each piece, turn a differential equation into a large sparse system of linear equations.

For a triangular piece, the entire physics collapses into one formula for how strongly two corners of that triangle interact:

𝐴 is the triangle’s area. 𝑏𝑖 and 𝑐𝑖 aren’t mysterious — they’re just numbers describing the edge directly across from corner 𝑖: bigger for a longer opposite edge, and they flip sign depending on which way that edge points. Feed in three corner coordinates and a material, and out comes a small table of these numbers, one for every pair of corners. That’s the whole trick. There’s no hidden layer between this formula and the code that computes it — about twenty lines of NumPy, every one of them readable.

Assemble that over every triangle, apply boundary conditions — a conductor is just a region held at a fixed voltage — solve the sparse system, and you have 𝑉 everywhere. The same twenty lines handle a flat plate, a circular cable cross-section, or any shape you can build from simple ones combined with union, intersection, and difference. The geometry never touches the physics.

The Deliberate Trade-Off

Here's the part most "look what I built" posts skip: an honest accounting of what it costs.

Real finite element tools mesh the actual geometry — triangles that hug a curved boundary exactly, dense where the field is doing something interesting, sparse where it isn't. Building that kind of mesh is genuinely hard, and it's exactly what pulls in the heavy dependencies — a geometry kernel, mesh-quality algorithms, compiled binaries — that make FEM tooling opaque in the first place.

I made a different choice: a plain Cartesian grid, split into triangles. No mesh generator. No native dependencies — just numpy.linspace and index arithmetic. A circle on a grid like this isn't a circle. It's a staircase of grid cells that looks roughly circular from far enough away.

That's not an oversight. It's the trade, made on purpose, to keep the whole tool inside three lines of pip install. The question is whether the trade is honest about its cost — and that's where most projects stop.

Proving It's Not a Toy

Before trusting any result, I checked the solver against a case with a known exact answer: plates spanning the entire simulation width, so no fringing is even geometrically possible. Just,

exactly, with nothing to approximate on the physics side.

The solver matched it to 0.0000% error, at every mesh resolution tested. That's the baseline: the finite element machinery has no bugs hiding in it. Whatever error shows up elsewhere comes from the mesh, not the math.

That check used a deliberately artificial setup — plates spanning the whole domain so fringing couldn’t happen — specifically so the formula would be exact and the test would be clean.
Real capacitors don’t look like that. So next I ran the two real examples, with real, finite edges — a coax cable filled with polyethylene, and a parallel plate with a glass slab across half the gap — through an actual convergence study instead of reporting one number and hoping.
The two reference values below aren’t the same kind of thing, on purpose: the coax formula is exact, so its gap is genuine mesh error. The parallel-plate formula, by contrast, is deliberately fringing-free — the classic no-fringing approximation, now extended to a glass-and-air sandwich instead of a single material — so its gap here isn’t a solver problem, it’s the real physics that formula was never built to capture.

The coax case increases at every step, though not at a perfectly steady rate — there’s a plateau partway through, still climbing but less steeply.

The parallel-plate case looks different again — across the four resolutions actually tested, it runs 102.0 → 101.9 → 98.8 → 101.8 pF/m, dipping in the middle before recovering, rather than settling smoothly toward the table’s two endpoints.

Two things are tangled together in that dip.
One is mundane: a fixed grid can only land exactly on some target sizes, and this particular resolution happens to round the capacitor’s own gap and dielectric thickness to slightly different values than its neighbors do — a small, real change in what’s actually being simulated, not an error.

The other is more interesting: a conductor’s sharp corner carries a genuine field singularity ( (𝐸∼𝑟^−1/3, for anyone curious), and independent meshes built at different resolutions aren’t strict refinements of each other, so the usual “finer mesh, better answer” guarantee doesn’t hold between them.

I could have shown just the two endpoints and let it look like ordinary convergence. Instead the program prints the whole sweep, and flags which of the two effects is at play at each point, rather than leaving the reader to guess.

A tool that only shows you its clean case isn’t more trustworthy than one that shows you both. It’s just quieter about its limits.

What This Is — and What It Isn't

This is	This is not
A from-scratch, mathematically correct FEM electrostatics solver	A replacement for COMSOL, Ansys Maxwell, or a production EDA tool
Small enough to read in one sitting — about 1,000 lines, three dependencies	An adaptive, conforming-mesh solver — curved boundaries are approximated, and quantifiably so
Validated against exact solutions, not just "it runs and the number looks plausible"	A 3D solver, or one that ingests arbitrary CAD geometry
Honest, in its own printed output, about exactly where and why it's approximate	Trying to be any of the above

If you need production-grade accuracy on complex geometry, this isn't your tool. If you want to actually see what a finite element solver is doing — every matrix entry, every approximation, every place it's wrong and by how much — this was built for exactly that.

Where It Goes From Here

The single biggest lever left on the table is the mesh itself. A few concrete next steps, in roughly the order they'd pay off:

Swap the Cartesian grid for a real unstructured, conforming mesh, fixing the staircased boundaries and under-resolved corners in one move — at the cost of the dependency-free simplicity that's the whole point of this version.
Short of that: grade the existing grid finer near sharp corners while staying dependency-free, which would tighten the parallel-plate numbers specifically without touching the coax case.
Add the classic extras once the mesh is solid: floating conductors, nonlinear and anisotropic dielectrics, 3D.

None of that happens silently. That's kind of the point.

The full code, the complete derivation, and both worked examples with figures are on: GitHub.

A Measured Formulation of the Shadow View of Collatz

Rat Wolf (ratwolf) — Fri, 03 Jul 2026 06:53:02 +0000

A Measured Formulation of the Shadow View of Collatz

A more careful description

A useful way to think about the Collatz map is to separate the orbit into two layers:

the parity pattern of the terms;
the arithmetic values that carry that pattern.

If one records the orbit as a sequence of even/odd decisions, the resulting binary word can be viewed as a symbolic encoding of the dynamics. In the accelerated version of the map, where one removes all factors of 2 immediately after each odd step, the evolution of the values becomes closely tied to this parity word. At that symbolic level, the dynamics look simple: the orbit is encoded by a binary sequence, and the next state is determined by the next bit of that sequence.

This is a genuine simplification of the picture, but it is a simplification of the representation, not a solution of the problem.

The Accelerated Collatz Formulation

This perspective uses the accelerated (or shortcut) Collatz map rather than the classical 3n+1 rule. In the accelerated version:

If n is even: next = n / 2
If n is odd: next = (3n + 1) / 2

By folding the immediate division by 2 after each 3n+1 step into a single operation, the accelerated map produces a shorter sequence in which every step corresponds directly to one bit in the parity word (0 = divide by 2, 1 = apply (3n+1)/2). For example, starting from 27 the classical orbit requires 111 steps, while the accelerated orbit reaches 1 in 71 steps. This formulation makes the parity shadow cleaner and more tightly coupled to the symbolic dynamics, while preserving the essential behavior of the original conjecture.

The Collatz Shadow Explorer

To actively investigate these dynamics, this repository includes a Python-based interactive desktop application (shadow_explorer.py). GitHub

Rather than viewing static plots, the application allows you to input starting values, adjust sequence parameters, and dynamically generate up to four concurrent analytical views of the parity shadow.

Analytical Views Included

The application provides several lenses through which to study the sequence:

Orbit + Shadow Track: The fundamental view. It overlays the numeric orbit (the arithmetic layer) with the binary parity sequence (the symbolic layer) on the same step axis.
Rolling Density Analysis: Displays the ratio of odd to even steps over a sliding window, highlighting localized phases of growth or decay.
: A histogram showing the frequency of $k$-length binary blocks (e.g., how often 011 or 101 occurs), revealing combinatorial constraints.
2D Turtle Walk (Fractal): Treats the parity word as directional instructions (0 = Left, 1 = Right) on a 2D plane. This maps the sequence to spatial geometry, often revealing structural symmetries hidden in the linear data.
Shannon Block Entropy: Measures the informational unpredictability of the sequence over time. It highlights chaotic phases versus highly ordered phases as the orbit collapses toward 1.
Fast Fourier Transform (FFT): Calculates the power spectrum of the $0-1$ signal, allowing you to check for dominant frequencies or periodic harmonics in the shadow word.
De Bruijn Transition Network: A directed network graph of parity block transitions. It visualizes the chronological "flow" of the sequence, essentially mapping the pathways that drain into the standard 4-2-1 attractor.
more see development folder on GitHub:

What this perspective is

This viewpoint is best understood as a reformulation. It highlights that:

the Collatz map can be studied through its parity pattern;
the parity pattern gives a symbolic description of the orbit;
the arithmetic part of the map is pushed into the encoding of the symbolic sequence.

In this sense, the problem can be seen as a question about which binary sequences arise from ordinary integers under this encoding, and how those sequences are constrained by arithmetic.

What this perspective is not

It is important not to overstate what this does. It does not:

solve the Collatz conjecture;
prove convergence for all starting values;
replace the classical numerical and probabilistic analysis;
introduce a fundamentally new dynamical system.

The main point is not that the problem becomes easy in this language. The point is that the structure becomes more legible: the symbolic part is simple, while the arithmetic constraints are the part that remains difficult.

Why it can still be useful

Even if it does not solve the conjecture, this perspective can be valuable for at least three reasons:

it gives a compact symbolic representation of the orbit;
it separates the simple combinatorial aspect from the arithmetic constraint;
it connects the problem with broader ideas from symbolic dynamics, signal processing, and 2-adic analysis.

So the idea is not that one has found a new mechanism for proving convergence. It is rather that one has found a different coordinate system—aided by computational visualization—in which the problem looks more structured and more interpretable.

Installation and Usage

To run the interactive explorer, you will need Python 3 installed along with a few standard scientific and GUI libraries.


bash

# Install required dependencies
pip install PyQt6 matplotlib numpy networkx

# Launch the explorer
python shadow_explorer.py

A Triplet-Based Parameterization for the Local Asymptotic Characterization of Polynomial Roots

Rat Wolf (ratwolf) — Tue, 31 Mar 2026 08:24:01 +0000

Abstract

This paper introduces a compact three-parameter framework for characterizing the local geometry of polynomial roots. For each root, the framework records its position, its algebraic multiplicity, and a newly defined quantity called the characteristic deflection distance 𝛿. This third parameter acts as a natural geometric scale: it measures how sharply or gradually the polynomial departs from zero in the immediate vicinity of the root, and it encodes the collective influence of all other roots through their distances from the one being analyzed.

The characteristic deflection distance generalizes the classical condition number of a simple root to roots of arbitrary multiplicity, and it allows direct geometric comparison across roots of different degrees. A key finding is that multiplicity alone does not determine geometric dominance — a lower-multiplicity root can have a larger spatial footprint than a higher-multiplicity one, depending on the global root configuration.

Symbolic and numerical Python implementations are provided, along with a worked example. The framework extends naturally to polynomials over the complex numbers.

DOI (Zenodo)

To visualize the spatial competition between these triplets, we map the 𝛿-Normalized Distance Field and the Newton Flow. This provides a "topographical" view of the polynomial's geometry where every root is evaluated against its own intrinsic scale.

GitHub

High-Precision Companion Matrix Root Finder

Rat Wolf (ratwolf) — Fri, 13 Mar 2026 06:53:50 +0000

High-Precision Companion Matrix Root Finder

GitHub
Medium

High-Precision Companion Matrix Root Finder

A robust polynomial solver leveraging arbitrary-precision floating-point arithmetic to compute roots of ill-conditioned polynomials. This tool primarily addresses rounding and finite-precision issues present in standard 64-bit floating-point implementations by constructing and solving a companion matrix within the mpmath environment.

The solver cannot eliminate the intrinsic mathematical ill-conditioning of certain polynomials, but it can significantly reduce the numerical errors caused by rounding, making the computed roots far more reliable.

It features modular code structure, residual checks, and robust input validation to improve usability and reliability.

Key Idea

The solver does not attempt to solve the fundamental conditioning problem of polynomial root finding.

Instead, it focuses on reducing numerical rounding errors using high‑precision arithmetic, which in practice leads to far more reliable root approximations for difficult polynomials.

Fast Chladni Pattern Simulation with a Simplified Membrane Model Using Vectorization and Stochastic Particles View

Rat Wolf (ratwolf) — Thu, 05 Mar 2026 06:25:48 +0000

Chladni Figures Simulation

This Python code simulates the nodal line patterns of Chladni figures, approximating Ernst Chladni’s 1787 experiments with vibrating plates by modeling the resonant modes of a square membrane. It captures the visual essence of these wave patterns without simulating the complex physics of elastic bending and the dynamics of particle motion.

Python Code, detailed simulation description, physics and mathematics see:

GitHub

Magnitude View: A 2D static plot showing the time-independent displacement magnitude ∣Z(x,y,f)∣ of the steady-state response at a driving frequency f matching the degenerate modes (3,5) and (5,3). A small damping factor (γ = 0.01) ensures these modes dominate. The magnitude highlights nodal lines and makes positive and negative antinodes indistinguishable, mimicking experimental Chladni patterns.

Phase View: This view displays the signed spatial amplitude of the steady-state response at a given frequency. The sign indicates whether a region oscillates in phase or π out of phase with a reference. Nodal lines appear where the amplitude vanishes. The membrane motion remains oscillatory in time; this view represents the spatial structure and phase relationships of that motion.

Sand View: The Sand View acts as a visually intuitive “nodal detector” derived from the same resonance-weighted field used in the magnitude and phase views.

Relation to Physical Experiments

In real Chladni experiments:

Grains bounce dynamically due to plate acceleration. Migration occurs through complex frictional and collisional processes.
Accumulation is influenced by time evolution and nonlinear interactions.

The Sand View abstracts these dynamics into a stationary probability model. It captures the final accumulation pattern without modeling transient motion.

Limitations

Field-based simulation; does not model particle dynamics.

Note: This simulation visualizes the nodal patterns of a unit square membrane L_x = L_y with fixed edges, governed by the wave equation, capturing the essence of Chladni figures. In real experiments, elastic plates follow complex bending physics, but the membrane model approximates their nodal patterns using simpler wave mechanics. Focusing on nodal patterns avoids the unresolved complexity of grain motion, highlighting the universal physics of resonance, independent of material properties.

Only finite max_mode included.
Damping gamma is uniform; real plates have non-uniform damping.
Geometry is idealized unit square membrane.
Uses a membrane model for computational efficiency, which simplifies the physics of a true plate with bending stiffness.

Simulation–Experiment Match Fidelity

Mode Category	Typical Example(s)	Match Fidelity	Notes
Exact Symmetric Modes	(2,2), (3,3), (4,4)	100%	Perfect match: Clean axial nodal lines align exactly with experiments.
Low-order Degenerate Pairs	(1,2)+(2,1), (1,3)+(3,1)	~50–60%	Partial match: Shows tilted/blended nodals instead of sharp crosses seen in some low-order plate patterns.
High-order Degenerate Pairs	(3,5)+(5,3), (1,6)+(6,1)	~90–95%	Excellent match: Complex, interior nodals closely mimic historical figures; diagonals may flip orientation.

Notes:

Symmetric modes m,m are exact (100%) and visually identical to experiments.

Asymmetric modes m ≠ n always form degenerate pairs in this square model, blending into symmetric patterns—both contribute fully e.g., (1,2) and (2,1) dominate together. No isolated "single" asymmetric modes appear due to symmetry.

Low-order pairs can look less crisp; high-order excel at reproducing rich, experimental-like topology, including accidental degeneracies where symmetric modes blend with pairs (e.g., (10,10) + (2,14)/(14,2)). Larger gamma adds asymmetry from nearby modes, enhancing realism.

Ratings prioritize nodal patterns (positions and shapes); the sim's blending captures the visual spirit of Chladni figures effectively, even if idealized.

Historical Context

Ernst Chladni (1756–1827), known as the father of acoustics, pioneered the study of vibrational modes by observing how elastic plates, sprinkled with sand or powder, formed Chladni figures—geometric patterns along nodal lines where the surface remains stationary. These patterns reveal standing waves and remain a subject of research today.

In his 1787 experiments:

An elastic plate (thin brass plate) was fixed at its center or edges.
It was vibrated with a violin bow at various frequencies.
Particles migrated from high-vibration areas to nodal lines, creating striking patterns.

Chladni figures laid foundations for acoustics, wave physics, and applications like mechanical engineering, with their nodal patterns inspiring both historical and modern scientific inquiry.

References

Chladni, Ernst Florens Friedrich. Entdeckungen über die Theorie des Klanges. Leipzig: Weidmanns Erben und Reich, 1787. ETH-Bibliothek Zürich,
Rar 5284. https://doi.org/10.3931/e-rara-4235. Public Domain Mark.
Herman, Russell. Vibrations of Rectangular Membranes. (2024, September 4) University of North Carolina Wilmington. CC BY‑NC‑SA. https://math.libretexts.org/@go/page/90264
Wikipedia contributors. “Ernst Chladni.” Wikipedia. https://en.wikipedia.org/wiki/Ernst_Chladni

An old Classic with new Perspective in 2D and 3D

Rat Wolf (ratwolf) — Wed, 04 Dec 2024 14:31:34 +0000

Historical Context

The Hopalong* attractor, authored by Barry Martin of Aston University in Birmingham, England [2], was popularized by A.K. Dewdney in the September 1986 issue of Scientific American. In Germany, the Hopalong attractor gained further popularity through an algorithm called HÜPFER in the article Psychotapeten published in Spektrum der Wissenschaft. [3].

*Nicknamed by A.K. Dewdney.

Python Code, detailed description and mathematics see:

GitHub

The Hopalong Attractor Functions

The mathematical definition of the Hopalong attractor is given by the following system of recursive functions, defined in [2]:

Where:

The sequence starts from the initial point (x0 , y0) = (0 , 0).
xn and yn represent the coordinates at the n-th iteration of the attractor.
a, b, and c are parameters influencing the attractor's dynamics.
sgn is the signum function. However, the programs use math.copysign(), which is defined as follows:

Computational Approach

The available Python programs calculate and visualize the attractor by iterating the described function system.

Goal:

Visual Representation:
The attractor is rendered as a density heatmap, with color intensity reflecting areas of high or low point concentration.
Efficiency:
- Large numbers of iterations are computed with low memory usage and high processing speed.
- Just-in-time (JIT) compilation is applied and supported by a low-complexity code structure.

Core Algorithm:

The algorithm consists of two main passes:

First Pass: Determines the spatial extent of the attractor trajectory.
Second Pass: Maps the sequentially generated trajectory points in continuous space directly to a discrete pixel grid, keeping track of the number of pixel hits used to generate the density heatmap matrix.

References

[1] J. Lansdown and R. A. Earnshaw (eds.), Computers in Art, Design and Animation. New York: Springer-Verlag, 1989.
e-ISBN-13: 978-1-4612-4538-4.

[2] Barry Martin, "Graphic Potential of Recursive Functions," in Computers in Art, Design and Animation [1], pp. 109–129.

[3] A.K. Dewdney, Psychotapeten, algorithm HÜPFER, in Spektrum der Wissenschaft: Computer Kurzweil.
Spektrum der Wissenschaft Verlagsgesellschaft mbH & Co., Heidelberg, 1988.
(German version of Scientific American).
ISBN-10: 3922508502, ISBN-13: 978-3922508502.

DEV Community: Rat Wolf (ratwolf)

Hopalong Map – fast, parallel, low-memory – two-pass checkpointed orbit decomposition

Contents

Quick start

The Hopalong map

The general approach (checkpoints and segments)

The problem

The idea in one picture

Checkpoints (flags) and segments

Why the result is exactly the same

What is stored, and how much

Lanes: keeping each core busy

Where the time goes

Counts vs. density

Counts (default)

Density

What the density option is not

Colour scales

Using it from Python

Configuration

Code map

Verification and benchmarks

Limitations and notes

One-pass Hopalong orbit density plot

Target

Implementation details

Orbit recurrence

Starting grid

Streaming histogram and dyadic growth

Output image

Tuning resolution with --bins

Result

Caveats

Use cases

Requirements

Usage

The Hopalong "Attractor" That Isn't One

A classic computer-art map has no attractor. Here is how to see its regular and chaotic regions with a few lines of Python and numba

The map

One orbit, one picture

Why it is not an attractor

Seeing all of phase space

A prettier version

What this is, and what it isn't

Try it yourself

References

Discreteness as the Late-Time Limit of a Smooth Flow

Abstract

The exact solution

I Built a Finite-Element Solver With Nothing to Hide

The Physics, Compressed

What "No Black Box" Actually Means

The Deliberate Trade-Off

Proving It's Not a Toy

What This Is — and What It Isn't

Where It Goes From Here

A Measured Formulation of the Shadow View of Collatz

A Measured Formulation of the Shadow View of Collatz

A more careful description

The Accelerated Collatz Formulation

The Collatz Shadow Explorer

Analytical Views Included

What this perspective is

What this perspective is not

Why it can still be useful

Installation and Usage

A Triplet-Based Parameterization for the Local Asymptotic Characterization of Polynomial Roots

High-Precision Companion Matrix Root Finder

High-Precision Companion Matrix Root Finder

Key Idea

Fast Chladni Pattern Simulation with a Simplified Membrane Model Using Vectorization and Stochastic Particles View

Chladni Figures Simulation

Limitations

Simulation–Experiment Match Fidelity

Historical Context

References

An old Classic with new Perspective in 2D and 3D

Historical Context

The Hopalong Attractor Functions

Computational Approach

Tuning resolution with `--bins`