Which Shape Encloses the Most Area for a Given Perimeter?
Imagine you have a fixed length of fencing and want to enclose the largest possible field. Should you build a square, a rectangle, a triangle — or something else entirely? This is the classical isoperimetric problem, one of the oldest optimization problems in mathematics, dating back to the legend of Queen Dido and the founding of Carthage.
The answer, proven rigorously in the 19th century, is the circle. Formally, the isoperimetric inequality states that for any simple closed curve with perimeter $L$ enclosing an area $A$:
$$
4\pi A \le L^2
$$
with equality if and only if the curve is a circle. Equivalently, we can define the isoperimetric ratio
$$
Q = \frac{4\pi A}{L^2}, \qquad 0 < Q \le 1
$$
where $Q = 1$ only for the circle, and $Q < 1$ for every other shape.
In this article, we verify this inequality computationally in three ways:
- Direct comparison of simple shapes (circle, square, triangle, rectangle) that all share the same perimeter.
- Numerical optimization: starting from an arbitrary “wavy” closed curve, we numerically deform it to maximize area while keeping perimeter fixed, and watch it converge to a circle.
- A 3D landscape of the isoperimetric ratio $Q$ as a function of shape perturbations, visually showing the circle sitting at the global maximum.
Mathematical Setup
We represent a closed curve in polar-like form:
$$
r(\theta) = 1 + \sum_{k=1}^{n} \left[ a_k \cos(k\theta) + b_k \sin(k\theta) \right], \qquad \theta \in [0, 2\pi)
$$
For such a curve, two classical formulas give us the quantities we need:
$$
A = \frac{1}{2}\int_0^{2\pi} r(\theta)^2 , d\theta
\qquad\text{(polar area formula)}
$$
$$
L = \int_0^{2\pi} \sqrt{r(\theta)^2 + r’(\theta)^2} ; d\theta
\qquad\text{(arc-length formula)}
$$
When all Fourier coefficients $a_k, b_k$ are zero, $r(\theta) \equiv 1$ and the curve is exactly the unit circle. The optimization problem becomes:
$$
\max_{a_k, b_k} ; A \quad \text{subject to} \quad L = L_0
$$
and the isoperimetric theorem predicts that the optimal solution is $a_k = b_k = 0$ for all $k$.
Source Code
1 | # ============================================================ |
Code Walkthrough
Part 1 — Comparing shapes by formula
We fix a perimeter $L = 4\pi$ and derive the side lengths of a square, an equilateral triangle, and a 2:1 rectangle that all have exactly this perimeter, using elementary formulas (e.g. for a square, $s = L/4$, $A = s^2$). We then compute each shape’s isoperimetric ratio $Q = 4\pi A / L^2$. Since these are closed-form formulas, this part is instantaneous — there is no computational bottleneck here at all.
Part 2 — Numerical optimization with SLSQP
This is the computational heart of the article. We describe an arbitrary closed curve using a truncated Fourier series in the radius $r(\theta)$. The function radius_and_derivative is fully vectorized: instead of looping over each $\theta$ value in Python (which would be slow), it uses np.outer to build a matrix of $\cos(k\theta)$ and $\sin(k\theta)$ values for all modes $k$ and all angles $\theta$ simultaneously, then collapses it with a single matrix-vector product (a @ cosk). This means evaluating the curve at hundreds of angles costs only a couple of NumPy matrix multiplications rather than a Python-level loop — the code stays fast even if you increase the number of Fourier modes or angular resolution.
curve_area and curve_perimeter implement the two closed-curve formulas from the math section using np.trapz for the numerical integration.
We then hand the problem to scipy.optimize.minimize with the SLSQP (Sequential Least Squares Programming) method, which is well suited to constrained nonlinear optimization with a small number of variables. The objective is -curve_area (since minimize always minimizes), and the perimeter constraint is expressed as an equality constraint perimeter(c) - L_TARGET = 0. A callback records the area at every iteration so we can plot the convergence path afterward.
Starting from a distinctly non-circular, asymmetric “wavy” shape, the optimizer should drive nearly all Fourier coefficients toward zero, and the isoperimetric ratio $Q$ should climb toward $1$.
Part 3 — 3D landscape of the isoperimetric ratio
Rather than optimize, this part directly maps out how $Q$ behaves for a simple one-parameter family of shapes, $r(\theta) = 1 + a\cos(k\theta)$, where $a$ is the perturbation amplitude and $k$ is the mode number (i.e., how many “lobes” the perturbation has). For every combination of $a$ and $k$ on a grid, we compute $Q$ using the same vectorized area/perimeter formulas — but this time processing an entire array of amplitudes at once via NumPy broadcasting (amps[:, None] * np.cos(...)), so the only actual Python-level loop is over the 8 mode numbers. This keeps the whole 3D scan running in a fraction of a second even with 1,500 integration points per curve.
The resulting surface should look like a ridge: $Q = 1$ exactly along $a = 0$ (the circle, regardless of $k$), and $Q$ decreases smoothly as $|a|$ grows in either direction — a direct visual proof that any deviation from a circle strictly decreases the area-to-perimeter efficiency.
Execution Results
Run the cell above in Google Colaboratory. Paste your outputs into the marked areas below.
PART 1 — console output
======================================================= PART 1: Areas for a fixed perimeter L = 12.5664 ======================================================= Circle Area = 12.5664 Q = 4*pi*A/L^2 = 1.0000 Square Area = 9.8696 Q = 4*pi*A/L^2 = 0.7854 Equilateral Triangle Area = 7.5976 Q = 4*pi*A/L^2 = 0.6046 Rectangle (2:1) Area = 8.7730 Q = 4*pi*A/L^2 = 0.6981
PART 1 — figure (shape comparison overlay)
PART 2 — console output
======================================================= PART 2: Numerical optimization result ======================================================= Optimizer success : True Iterations : 40 Initial shape Area=3.39905 Perimeter=7.39715 Q=0.78062 Final shape Area=3.14159 Perimeter=6.28319 Q=1.00000 Final Fourier coefficients (should be near 0): [-0. 0. -0. -0. -0. -0. -0. -0.]
PART 2 — figure (initial vs. optimized shape + convergence curve)
PART 3 — console output
======================================================= PART 3: 3D surface summary ======================================================= Maximum Q on the grid : 1.000000 (theoretical max = 1.0) Location of max Q : amplitude=0.000, mode=1
PART 3 — figure (3D isoperimetric ratio surface)
Why This Matters
The isoperimetric problem isn’t just a mathematical curiosity — it explains why soap bubbles are spherical (they minimize surface area for a given enclosed volume), why cross-sections of blood vessels tend toward circularity to minimize the energy needed to maintain their boundary, and why many engineered enclosures (pipes, tanks, cross-sections of beams) default to circular or near-circular shapes for material efficiency. The same underlying principle — the circle as the unique optimizer of the area-to-perimeter trade-off — resurfaces across physics, biology, and engineering. What we’ve done here numerically is exactly what physical systems do continuously: relax toward the shape that minimizes “boundary cost” for a given “enclosed quantity.”
