ellipj and elliptic12
routines compute both exactly.
The Dubins Car and the Group SE(2)
A car at position $(x, y)$ heading in direction $\theta$ satisfies:
$$\dot x = \cos\theta,\quad \dot y = \sin\theta,\quad \dot\theta = u,$$where the control $u$ is the steering rate (curvature). If we want to steer from state $(x_0, y_0, \theta_0)$ to $(x_1, y_1, \theta_1)$ while minimising total squared curvature
$$\int_0^L u^2\,ds = \int_0^L \kappa^2(s)\,ds,$$the Pontryagin maximum principle (PMP) tells us the optimal control satisfies an ODE equivalent to the nonlinear pendulum. The state space of this system is the Lie group $\mathrm{SE}(2) = \mathbb{R}^2 \rtimes S^1$ — the group of orientation-preserving isometries of the plane.
The PMP Hamiltonian along an extremal is
$$H = \frac{h_1^2 + h_2^2}{2},$$where $h_1 = p_x\cos\theta + p_y\sin\theta$ and $h_2 = p_\theta$ are the "body-frame" components of the costate. Hamilton's equations give, after the substitution $h_1 = r\cos\varphi$, $h_2 = r\sin\varphi$:
$$\dot\varphi = -r\cos\varphi \cdot \text{(bracket term)},$$which reduces to the nonlinear pendulum $\ddot\varphi = -c\sin\varphi$ for a constant $c > 0$. The costate component $h_2 = \kappa$ is the curvature of the optimal path.
Euler's Elastica and Jacobi Elliptic Functions
The curvature of the optimal path satisfies Euler's elastica equation (Euler 1744):
$$2\kappa'' + \kappa^3 - \lambda\kappa = 0,\qquad \lambda > 0,$$which is a conservative ODE with energy $E = (\kappa')^2 + V(\kappa)$, $V = -\lambda\kappa^2/2 + \kappa^4/4$. The potential $V$ has two wells at $\kappa = \pm\sqrt{\lambda}$ and a local maximum at $\kappa = 0$. Solutions split into three topological classes:
- Inflectional elastica ($k < 1$): $\kappa$ oscillates across zero (two-well oscillation through the barrier). $$\kappa(s) = 2k\,\mathrm{sn}(s\mid k^2).$$
- Non-inflectional elastica ($m < 1$): $\kappa$ stays positive (oscillation within one well). $$\kappa(s) = 2\,\mathrm{dn}(s\mid m).$$
- Euler spiral / Cornu spiral (separatrix, $k = 1$): $\kappa$ is a positive pulse that decays to zero at infinity. $$\kappa(s) = 2\,\mathrm{sech}(s) = 2\,\mathrm{dn}(s\mid 1).$$
The spatial period of $\kappa(s)$ — the arc length between successive inflections — is
$$T_\kappa = 4K(k^2)\quad\text{(inflectional)},\qquad T_\kappa = 2K(m)\quad\text{(non-inflectional)},$$diverging as $k \to 1^-$ (the separatrix has infinite period — it is the Euler spiral, which spirals into two point attractors at $\pm\infty$). The position is then recovered by integrating $\theta(s) = \int_0^s \kappa(t)\,dt$ and $(x,y)(s) = \int_0^s(\cos\theta,\sin\theta)\,dt$.
Interactive: The Elastica Family
Figure 1 shows $\kappa(s)$ (left) and the corresponding curve $(x(s), y(s))$ (right) as a function of the modulus $k$. Drag the slider to morph between curve types.
The Euler / Cornu Spiral (Separatrix)
At $k = 1$ the inflectional and non-inflectional families both converge to the Euler spiral (also called the Cornu spiral or clothoid). Its curvature $\kappa(s) = 2\,\mathrm{sech}(s)$ is a positive pulse centred at $s = 0$, and it is the only elastica that changes from being "more curved" to "less curved" (the curvature maximum is at $s=0$, and $\kappa \to 0$ as $s \to \pm\infty$).
In the context of the Moiseev–Sachkov paper, this curve corresponds to the heteroclinic orbit (separatrix) of the pendulum, connecting the unstable fixed points at $\varphi = \pm\pi$. In highway engineering and railway design the Euler spiral is used as a transition curve exactly because its curvature increases linearly with arc length — smoothly connecting a straight section to a circular arc.
Spatial Period Governed by $K(k^2)$
Figure 3 shows how the spatial period $T(k) = 4K(k^2)$ (inflectional) and $2K(k^2)$ (non-inflectional) depend on the modulus $k$. The key features:
- At $k = 0$: $T \to 2\pi$ (inflectional limit) or $\pi$ (non-inflectional) — the circular limit.
- At $k \to 1^-$: $T \to +\infty$ — the separatrix has infinite period. Formally $K(k^2) \sim \ln(4/\sqrt{1-k^2})$, a logarithmic divergence.
- Near $k = 0.9$: the period has already doubled relative to the circular limit — the curve is noticeably "slower" (longer wavelength) than a circle.
This is the same $K(k^2)$ divergence that controls the libration period of the nonlinear pendulum and the Maxwell cut time of the SE(2) geodesic problem.
Geodesic Family from a Fixed Base Point
Figure 4 shows the exponential map of SE(2): all inflectional elastica sharing the same initial point and heading, for $k = 0.1, 0.2, \ldots, 0.97$ (blue–green gradient) and the Euler spiral $k = 1$ (red). This is the structure that Sachkov and Moiseev analyse to locate the Maxwell strata — the loci where two geodesics of equal length first coincide, marking the end of global optimality.
Petitot's Model of Primary Visual Cortex V1
In 2003 Jean Petitot proposed that the functional architecture of primary visual cortex V1 implements a contact structure on SE(2). Each neuron in V1 is tuned to a local position and orientation: the cortex is not a 2D image map but a contact bundle over the visual field. When a contour is interrupted — an occlusion, or a gap between two collinear segments — V1 fills it in by activating neurons along the shortest path in this contact geometry.
Citti and Sarti (2006) made the connection precise: the horizontal cortico-cortical
connections in V1 implement exactly the sub-Riemannian metric on SE(2), and the
"association fields" of Field, Hayes & Hess (1993) — the empirical patterns of
which oriented edge pairs humans most readily group together — correspond to
sub-Riemannian geodesics. Those geodesics are Euler's elastica. Their curvature
is $\kappa(s) = 2k\,\mathrm{sn}(s\mid k^2)$, evaluated using the
ellipj function.
Moiseev & Sachkov (2010) — arXiv:0807.4731 — then completely characterised the Maxwell strata for this problem: the sets where the exponential map fails to be injective, i.e., where the visual-completion "decision" becomes ambiguous because two equally-short curves exist between the given boundary conditions. The first Maxwell time equals the period $4K(k^2)/\omega_0$, computed via the same integral.
The connection drove the original development of the elliptic package
(2003–2004): fast, vectorised evaluation of $\mathrm{sn}, \mathrm{cn}, \mathrm{dn}$
and $K(m), E(m)$ was needed to sweep large families of geodesics for the Maxwell
strata analysis.
| Modulus $k$ | Curve type | Period $T = 4K(k^2)$ | $\kappa(s)$ |
|---|---|---|---|
| $k \to 0$ | Nearly circular arc | $\to 2\pi$ | $\approx 0$ (straight) or const. (circle) |
| $k = 0.7$ | S-curve elastica | $\approx 8.7$ | $1.4\,\mathrm{sn}(s\mid 0.49)$ |
| $k = 0.9$ | Deep S-curve | $\approx 14.8$ | $1.8\,\mathrm{sn}(s\mid 0.81)$ |
| $k = 1$ | Euler spiral (sep.) | $\to\infty$ | $2\,\mathrm{sech}(s)$ |
Code
The Python snippet below sweeps the inflectional elastica family, computes the
spatial curve for each modulus using ellipj and elliptic12,
and evaluates the period via ellipticK — all vectorised over a batch of
$k$-values.
import numpy as np
from scipy.integrate import cumulative_trapezoid
from elliptic import ellipj, elliptic12 # moiseevigor/elliptic
def elastica_inflectional(k: float, n_periods: float = 1.0, N: int = 2000):
"""
Compute one inflectional elastica with modulus k ∈ (0, 1).
κ(s) = 2k · sn(s | k²), period T = 4K(k²)
Returns arrays s, kappa, theta, x, y.
"""
m = k * k
# Complete elliptic integral K(m) for the period
Km, _, _ = elliptic12(np.array([np.pi / 2]), np.array([m]))
T = 4 * Km[0]
s = np.linspace(-n_periods * T / 2, n_periods * T / 2, N)
# Jacobi sn evaluated at all arc-length points simultaneously (GPU-ready)
sn_vals, cn_vals, dn_vals = ellipj(s, m)
kappa = 2 * k * sn_vals
# Integrate θ(s) = ∫κ ds and (x,y) = ∫(cos θ, sin θ) ds
theta = cumulative_trapezoid(kappa, s, initial=0.0)
x = cumulative_trapezoid(np.cos(theta), s, initial=0.0)
y = cumulative_trapezoid(np.sin(theta), s, initial=0.0)
return s, kappa, theta, x, y
def elastica_noninflectional(m: float, N: int = 2000):
"""
Non-inflectional elastica: κ(s) = 2·dn(s | m), m ∈ (0, 1).
Period T = 2K(m).
"""
Km, _, _ = elliptic12(np.array([np.pi / 2]), np.array([m]))
T = 2 * Km[0]
s = np.linspace(-T / 2, T / 2, N)
_, _, dn_vals = ellipj(s, m)
kappa = 2 * dn_vals
theta = cumulative_trapezoid(kappa, s, initial=0.0)
x = cumulative_trapezoid(np.cos(theta), s, initial=0.0)
y = cumulative_trapezoid(np.sin(theta), s, initial=0.0)
return s, kappa, theta, x, y
# ── Sweep k-values and report period ──────────────────────────────────
k_vals = np.linspace(0.1, 0.95, 18)
m_vals = k_vals ** 2
# Vectorised K(m) for all k at once — the GPU payoff
Km_all, _, _ = elliptic12(np.full(len(m_vals), np.pi / 2), m_vals)
T_all = 4 * Km_all
print(f"{'k':>6} {'K(k²)':>9} {'Period T':>10}")
for k, Km, T in zip(k_vals, Km_all, T_all):
print(f"{k:6.3f} {Km:9.5f} {T:10.5f}")
% Euler's Elastica via Jacobi elliptic functions
% Requires the elliptic toolbox: https://github.com/moiseevigor/elliptic
function [s, kappa, theta, x, y] = elastica_inflectional(k, n_periods, N)
% Inflectional elastica: kappa(s) = 2k * sn(s | k^2)
% k in (0, 1), period T = 4*K(k^2)
if nargin < 2, n_periods = 1; end
if nargin < 3, N = 2000; end
m = k^2;
Km = ellipticK(m); % K(m) — vectorised, GPU-ready
T = 4 * Km;
s = linspace(-n_periods*T/2, n_periods*T/2, N)';
[sn, ~, ~] = ellipj(s, m); % Jacobi sn for all s at once
kappa = 2 * k * sn;
% Integrate theta and (x, y) via cumulative trapezoid
ds = s(2) - s(1);
theta = cumtrapz(s, kappa);
x = cumtrapz(s, cos(theta));
y = cumtrapz(s, sin(theta));
end
% ── Sweep moduli and print period table ───────────────────────────────
k_vals = linspace(0.1, 0.95, 18);
m_vals = k_vals .^ 2;
Km_all = ellipticK(m_vals); % vectorised call — one line for 18 values
T_all = 4 * Km_all;
fprintf('%6s %9s %10s\n', 'k', 'K(k²)', 'Period T');
for i = 1:numel(k_vals)
fprintf('%6.3f %9.5f %10.5f\n', k_vals(i), Km_all(i), T_all(i));
end
% ── Plot a family of elastica ─────────────────────────────────────────
figure; hold on; axis equal;
cmap = parula(numel(k_vals));
for i = 1:numel(k_vals)
[~, ~, ~, x, y] = elastica_inflectional(k_vals(i));
plot(x - x(1), y - y(1), 'Color', cmap(i,:), 'LineWidth', 1.5);
end
title('Inflectional elastica family, \kappa(s) = 2k\cdotsn(s | k^2)');
xlabel('x'); ylabel('y');
colorbar; clim([k_vals(1) k_vals(end)]);
References
- I. Moiseev & Yu. L. Sachkov (2010). "Maxwell strata in sub-Riemannian problem on the group of motions of a plane." ESAIM: Control, Optimisation and Calculus of Variations 16(2): 380–399. arXiv:0807.4731
- J. Petitot (2003). "The neurogeometry of pinwheels as a sub-Riemannian contact structure." Journal of Physiology–Paris 97(2–3): 265–309.
- G. Citti & A. Sarti (2006). "A cortical based model of perceptual completion in the roto-translation space." J. Math. Imaging and Vision 24(3): 307–326.
- Yu. L. Sachkov (2011). "Cut locus and optimal synthesis in the sub-Riemannian problem on the group of motions of a plane." ESAIM: COCV 17(4): 293–321. arXiv:0903.0727
- L. Euler (1744). Methodus Inveniendi Lineas Curvas Maximi Minimive Proprietate Gaudentes. Additamentum I: De curvis elasticis.
- D. A. Singer (1999). "Curves whose curvature depends on distance from the origin." The American Mathematical Monthly 106(9): 835–841.
- D. B. Henry (2004). "The Euler spiral: a mathematical history." UC Berkeley Technical Report EECS-2008-111.