diff --git a/_posts/2026-04-25-geometry-of-seeing-visual-cortex-se2.md b/_posts/2026-04-25-geometry-of-seeing-visual-cortex-se2.md index d0d63f1f4f4..f9017a298cb 100644 --- a/_posts/2026-04-25-geometry-of-seeing-visual-cortex-se2.md +++ b/_posts/2026-04-25-geometry-of-seeing-visual-cortex-se2.md @@ -119,7 +119,8 @@ Mapped over a patch of cortex, the preferred orientations rotate continuously completing a full $\pi$-rotation over a distance of about 1 mm. The resulting pattern of orientation preferences is called an *orientation map*, and it exhibits a characteristic structure of *pinwheels*: point singularities -around which the preferred orientation rotates by $\pm\pi/2$. +of topological charge $\pm\tfrac12$ — encircle one and the preferred orientation +sweeps through a full half-turn, $\pm\pi$. satisfies the Hörmander (bracket-generating) condition. By the -**Chow–Rashevskii theorem** (1938, 1938) any two configurations in +**Chow–Rashevskii theorem** (Rashevskii 1938, Chow 1939) any two configurations in $\mathrm{SE}(2)$ can therefore be joined by a *horizontal* path — a curve whose velocity lies in $\mathcal{H}$ at every point. Geometrically, $\mathcal{H}$ is a **contact structure** on $\mathrm{SE}(2)$, and Petitot's @@ -496,11 +501,12 @@ The spatial projection $(x(t), y(t))$ of the solution is the **perceptually comp contour** that the visual system infers between two oriented line elements $(x_0,\theta_0)$ and $(x_1,\theta_1)$. -The simplification $u_1 = 1$ (unit forward speed) reduces the problem to minimising -$$\int_0^L \kappa^2(s)\,ds$$ -where $L$ is arc length and $\kappa = u_2/u_1 = \dot\theta$ is the signed curvature. -This is the **Euler elastica functional**: the total squared bending energy of the -projected curve. +With the normalisation $u_1 = 1$ (unit forward speed), $s$ becomes arc length and +$\kappa = u_2/u_1 = \dot\theta$ is the signed curvature. The sub-Riemannian extremals +then project onto the critical curves of the **Euler elastica functional** +$$\int_0^L \kappa^2(s)\,ds,$$ +the total squared bending energy of the projected curve — the classical bridge, +made precise in Part 2, between geodesics on SE(2) and Euler's elastica. @@ -527,9 +533,9 @@ projected curve. $k \in (0, 1)$ — inflectional, $\kappa(s) = 2k\,\mathrm{cn}(s\mid k^{2})$; $k = 1$ — the Euler / Cornu spiral, the separatrix $\kappa(s) = 2\,\mathrm{sech}\,s$; - $k > 1$ — non-inflectional, parametrised by - $m = 2 - k \in (0, 1)$ with $\kappa(s) = 2\,\mathrm{dn}(s\mid m)$ - (closed-loop curves with one-signed curvature). + $k > 1$ — non-inflectional, $\kappa(s) = 2\,\mathrm{dn}(s\mid m)$ + with its own modulus $m \in (0, 1)$ (closed-loop curves with one-signed + curvature); the slider just sweeps $m$ downward as $k$ runs past 1. Drag the slider to set the maximum $k$ rendered; the vertical bar on the right is the colour scale, with the red tick marking the separatrix $k = 1$. As $k \to 1$ the period $4K(k^{2})$ diverges and the inflectional @@ -583,8 +589,10 @@ $L_{\mathrm{SR}}$: equal or smaller $L_{\mathrm{SR}}$. At a Maxwell point, two distinct globally-optimal geodesics meet with *exactly* the same SR length. -For a sub-Riemannian manifold as symmetric as SE(2), the cut and Maxwell loci coincide -(this is part of what Sachkov and I proved in arXiv:0807.4731). +For a sub-Riemannian manifold as symmetric as SE(2), the cut and Maxwell loci are +tightly linked: characterising the Maxwell strata is the main result of +arXiv:0807.4731 (Sachkov and I), and pinning the cut locus down exactly is the +subject of arXiv:0903.0727. The first point on the cut locus along a given geodesic is the **cut time** $t_\mathrm{cut}$, and it equals the first time the exponential map is no longer injective. diff --git a/_posts/2026-04-28-geometry-of-seeing-elastica-jacobi.md b/_posts/2026-04-28-geometry-of-seeing-elastica-jacobi.md index 97c7869e09c..01a389537e1 100644 --- a/_posts/2026-04-28-geometry-of-seeing-elastica-jacobi.md +++ b/_posts/2026-04-28-geometry-of-seeing-elastica-jacobi.md @@ -85,9 +85,13 @@ The **Hamiltonian** for the PMP maximisation condition is $$H = h_1 u_1 + h_2 u_2 - \nu \sqrt{u_1^2 + u_2^2},$$ -where $\nu \in \{0, \tfrac{1}{2}\}$ is the abnormality constant. -For **normal extremals** ($\nu = \tfrac{1}{2}$), maximising over $u_1, u_2$ -gives the **normal Hamiltonian** +where $\nu \in \{0, 1\}$ is the abnormality constant ($\nu = 1$ for +**normal** extremals, $\nu = 0$ for abnormal; Appendix A3 carries out the +maximisation in full). +Minimising length is equivalent to minimising the energy +$\tfrac12\!\int(u_1^2 + u_2^2)\,dt$ over a fixed interval (Cauchy–Schwarz), +and for normal extremals the PMP maximisation then +delivers the **normal (sub-Riemannian) Hamiltonian** $$\mathcal{H}_n = \tfrac{1}{2}\bigl(h_1^2 + h_2^2\bigr).$$ @@ -191,15 +195,14 @@ delivered for the SR-on-SE(2) costate phase. The pendulum equation lives on the costate, but the figures below plot the **curvature $\kappa(s)$ of the projected plane curve in its own Euclidean -arc length $s$**. These two are tied together by a quick computation. +arc length $s$** — and $s$ is not the costate time $t$: the two are tied +together by $ds = |h_1|\,dt$. -In unit-speed parametrisation $h_1^{2}+h_2^{2}=1$, the projected speed is -$|h_1| = |\sin(\varphi/2)|$ and the heading derivative is $\dot\theta = h_2 = \cos(\varphi/2)$. -Reparametrising by Euclidean arc length $s$ with $ds = |h_1|\,dt$, and using -$\dot\varphi = \pm 2 h_3$, a short calculation gives an explicit closed-form -$\kappa(s)$ in each regime — the three Jacobi-elliptic curvature profiles -below. They satisfy the **elastica curvature ODE** (the Duffing form -equivalent to the pendulum) +Pushing the pendulum solution through that reparametrisation is the technical +heart of the reduction, and it is carried out in full in Appendix A3. The +outcome is what matters here: in each regime $\kappa(s)$ comes out in closed +form as one of the three Jacobi-elliptic profiles below, and each one solves +the **elastica curvature ODE** — the Duffing form equivalent to the pendulum, $$\kappa''(s) + \tfrac{1}{2}\kappa(s)^{3} - \mu\,\kappa(s) \;=\; 0,$$ @@ -374,8 +377,9 @@ The browser figures on this page use the same AGM algorithm implemented in For the Jacobi functions themselves: ```python -from elliptic import ellipj +from elliptic import ellipj, ellipticK +K = ellipticK(k**2) # quarter-period of sn, cn s = np.linspace(-2*K, 2*K, 800) sn, cn, dn = ellipj(s, k**2) kappa = 2 * k * cn # curvature of inflectional elastica @@ -417,9 +421,9 @@ In the `elliptic` package: ```python from elliptic import elliptic12 -phi = np.arcsin(k * sn) -E_vals, F_vals = elliptic12(phi, k**2) # E(φ|k²) and F(φ|k²) -x = 2 * (E_vals - F_vals / 2) # exact formula from Sachkov (2011) +am = np.arcsin(sn) # Jacobi amplitude am(s | k²) +F_vals, E_vals = elliptic12(am, k**2) # F(am | k²) = s and E(am | k²) +x = 2 * E_vals - F_vals # x(s) = 2 E(am(s)|k²) − s (Sachkov 2011) ``` The full closed-form expressions — due to Sachkov (2011) — express every @@ -435,8 +439,9 @@ A few landmarks worth noting: - **$k = 0.1$** (inflectional): nearly straight, very gentle curvature oscillation. The curve barely bends before straightening again. -- **$k \approx 0.71$** (inflectional): the "figure-eight" lemniscate — the curve - crosses itself once per period and the endpoints of one period coincide. +- **$k \approx 0.909$** (inflectional): the "figure-eight" lemniscate — at the + modulus where $2E(k^2) = K(k^2)$, the curve crosses itself once per period and + the endpoints of one period coincide. This is the **Maxwell stratum** for the symmetric geodesics (Part 3). - **$k \to 1^-$** (inflectional → Euler spiral): the period $4K(k^2)$ diverges and diff --git a/_posts/2026-05-01-geometry-of-seeing-A1-lie-groups.md b/_posts/2026-05-01-geometry-of-seeing-A1-lie-groups.md index def3949434c..d8b9e2e21e5 100644 --- a/_posts/2026-05-01-geometry-of-seeing-A1-lie-groups.md +++ b/_posts/2026-05-01-geometry-of-seeing-A1-lie-groups.md @@ -32,7 +32,7 @@ published: false Part 1 of the series uses the language of Lie groups and Lie algebras as if it were standard furniture: $\mathrm{SE}(2)$, $\mathfrak{se}(2)$, left-invariant vector fields $X_1 = \cos\theta\,\partial_x + -\sin\theta\,\partial_y$, the bracket $[X_1, X_2] = X_3$, the exponential map. +\sin\theta\,\partial_y$, the bracket $[X_1, X_2] = -X_3$, the exponential map. This appendix builds those objects from scratch. Read it once and Part 1 becomes a calmer text. All three figures below are powered by the same $\mathrm{SE}(2)$ matrix exponential routine that the @@ -124,31 +124,40 @@ $G$ tangent to $E_i$ at the identity, multiply it on the left by $g$, and read off its velocity at $g$. -There is a canonical way to push the Lie algebra around the group: at any -$g \in G$ define +There is a canonical way to push any algebra element $E \in \mathfrak g$ +around the group: at any $g \in G$ define -$$X_i(g) \;:=\; (dL_g)_e (E_i),$$ +$$\widetilde E(g) \;:=\; (dL_g)_e (E),$$ -where $L_g(h) = gh$ is left-multiplication. $X_i(g)$ is the pushforward of -the abstract algebra element $E_i$ to the tangent space at $g$ along the -left-translation. These are the **left-invariant vector fields**. They -satisfy $X_i(gh) = (dL_g)_h X_i(h)$, i.e. they look the same in every -left-translated frame. +where $L_g(h) = gh$ is left-multiplication. $\widetilde E(g)$ is the +pushforward of +$E$ to the tangent space at $g$ along the left-translation. The vector +fields obtained this way are the **left-invariant vector fields**; they +satisfy $\widetilde E(gh) = (dL_g)_h \widetilde E(h)$, i.e. they look the +same in every left-translated frame. -Computing them in the chart $(x, y, \theta)$ is mechanical: +Computing the pushforwards in the chart $(x, y, \theta)$ is mechanical: $L_g(x', y', \theta') = (\,x + x'\cos\theta - y'\sin\theta,\;\; y + x'\sin\theta + y'\cos\theta,\;\; - \theta + \theta')$, so + \theta + \theta')$. Following Part 1's index +convention — $X_1$ forward, $X_2$ rotation, $X_3$ sideways — the three +basis fields are $$\boxed{\; - X_1 \;=\; \cos\theta\,\partial_x + \sin\theta\,\partial_y, \qquad - X_2 \;=\; \partial_\theta, \qquad - X_3 \;=\; -\sin\theta\,\partial_x + \cos\theta\,\partial_y .\;}$$ + X_1 = (dL_g)_e E_1 = \cos\theta\,\partial_x + \sin\theta\,\partial_y, \quad + X_2 = (dL_g)_e E_3 = \partial_\theta, \quad + X_3 = (dL_g)_e E_2 = -\sin\theta\,\partial_x + \cos\theta\,\partial_y .\;}$$ + +Note the index shuffle: $X_2$ comes from the *rotation* generator $E_3$ and +$X_3$ from the *translation* generator $E_2$, because Part 1 numbers the +frame (forward, rotation, sideways) while the $E$-basis is numbered +(translate-$x$, translate-$y$, rotate). These are the same three vector fields Part 1 §3 introduced. Now you know -where they come from: they are the basis of $\mathfrak{se}(2)$, parallel- -transported across the group by left-multiplication. They form an -orthonormal frame for the *Cartan-Killing geometry* on $\mathrm{SE}(2)$. +where they come from: they are the basis of $\mathfrak{se}(2)$, +parallel-transported across the group by left-multiplication. Declaring +them orthonormal is exactly how Part 1 puts its left-invariant metric on +$\mathrm{SE}(2)$. @@ -229,14 +238,16 @@ Two warnings worth absorbing. - This basis $\{E_1, E_2, E_3\}$ of $\mathfrak{se}(2)$ is *different* from the left-invariant frame $\{X_1, X_2, X_3\}$ used in Part 1. At the - identity $X_i(e) = E_i$, but at a generic $g \in \mathrm{SE}(2)$, - $X_i(g) \neq E_i$ — the LI vector fields are not constant in + identity $X_1(e) = E_1$, $X_2(e) = E_3$, $X_3(e) = E_2$ (the index + shuffle above); and at a generic $g \in \mathrm{SE}(2)$ the $X_i(g)$ are + no longer the constant matrices $E_j$ — the LI vector fields vary in coordinates. -- For the LI vector fields the relevant bracket is the **vector-field - commutator** below, *not* the matrix commutator of their constant - identity-values. The vector-field bracket of $X_1$ with $X_2$ produces - $\pm X_3$ — the missing sideways direction — even though the matrix - bracket of $E_1$ with $E_2$ vanishes. +- For left-invariant fields the vector-field commutator and the matrix + commutator *agree*: $[X_i, X_j]_{\text{v.f.}}(e) = [\,X_i(e), X_j(e)\,]$ + in $\mathrm{Mat}_3$. What you must *not* do is take the matrix commutator + of the coordinate expressions $X_i(g)$ as if they were constant — they + are not. Done correctly, $[X_1, X_2] = -X_3$ picks out the missing + sideways direction, matching the matrix bracket $[E_1, E_3] = -E_2$. **(ii) Vector-field commutator.** For two vector fields acting on smooth functions $f$, @@ -257,10 +268,10 @@ $$[X_1, X_2] f \;=\; -(\partial_\theta\cos\theta)\partial_x f \;=\; \sin\theta\,\partial_x f - \cos\theta\,\partial_y f \;=\; -X_3 f.$$ -So $[X_1, X_2] = -X_3$ as left-invariant vector fields. The sign is a -convention: Part 1 used the opposite sign convention to land at -$[X_1, X_2] = +X_3$, and we will switch to that in §4 below where it matters. -Either way, the bracket is $\pm X_3$ — non-zero, in the missing sideways +So $[X_1, X_2] = -X_3$ as left-invariant vector fields — the same sign +Part 1 carries. (Writing it the other way round, $[X_2, X_1] = +X_3$, is +the only freedom here; it is a bookkeeping choice, not a real one.) Either +way the bracket is non-zero and points along $X_3$, the missing sideways direction. **(iii) Closing-defect interpretation.** Flow along $X$ for time @@ -377,16 +388,18 @@ through Figure A1.1: $\omega = 0$ gives a straight line, $\omega \neq 0$ gives a circle whose centre is offset from the origin by the screw "axis" $a / \omega$. -### The 1-parameter subgroups *are* the geodesics of the Cartan-Killing metric +### 1-parameter subgroups and bi-invariant geodesics -If you put a left-invariant Riemannian metric on $G$ that is also right- -invariant (a "bi-invariant" metric, which on $\mathrm{SE}(2)$ exists), then -the geodesics through $e$ are exactly the 1-parameter subgroups -$t \mapsto \exp(tX)$. This is **not** the situation in Part 1: Part 1 uses -a *sub-Riemannian* metric that is left-invariant but not right-invariant, -and in that geometry the geodesics are **not** generally -$\exp(tX)$ — they are Euler's elastica. Appendix A5 explains how the SR -exponential map differs from this group exponential. +On a group that carries a *bi-invariant* metric — one invariant under both +left and right translation — the geodesics through $e$ are exactly the +1-parameter subgroups $t \mapsto \exp(tX)$. Compact groups and $\mathbb R^n$ +have such metrics; $\mathrm{SE}(2)$ does **not** — its adjoint action is +non-compact, which is the standard obstruction. So even a Riemannian story +on $\mathrm{SE}(2)$ would not make $\exp(tX)$ geodesic. Part 1 goes further +still: its metric is *sub-Riemannian*, left-invariant but not +right-invariant, and its geodesics are not $\exp(tX)$ at all — they are +Euler's elastica. Appendix A5 explains how the SR exponential map differs +from this group exponential. ## Adjoint and coadjoint actions @@ -421,10 +434,13 @@ $$\mathcal O_c \;:=\; \{(h_1, h_2, h_3) : h_1^2 + h_2^2 = c\},$$ i.e. **vertical cylinders** in $(h_1, h_2, h_3)$-space (plus a degenerate 1-point orbit at $h_1 = h_2 = 0$ for each value of $h_3$). Part 2 §1 -discovered this structure organically: the costate of the SR geodesic -problem evolves on the cylinder $h_1^2 + h_2^2 = 2\mathcal H$, with $h_3 = -\omega_0$ constant — exactly Lie–Poisson dynamics on $\mathfrak{se}(2)^{\ast}$. -Appendix A3 will derive that flow from the PMP. +discovered this structure organically. A caution on indices: Part 2 labels +the costate in the left-invariant frame $\{X_1, X_2, X_3\}$, so its +$(h_1, h_2, h_3)$ are this appendix's $(h_1, h_3, h_2)$. In the $E$-basis +the costate stays on a fixed coadjoint cylinder $h_1^2 + h_2^2 = c$ — the +Casimir is conserved — while the SR Hamiltonian +$\mathcal H = \tfrac12(h_1^2 + h_3^2)$ drives it around that cylinder. +Appendix A3 derives the flow from the PMP. @@ -449,14 +465,13 @@ Appendix A3 will derive that flow from the PMP.
Figure A1.3. The coadjoint orbits of $\mathrm{SE}(2)$ are cylinders $h_1^2 + h_2^2 = c$ in the dual space $\mathfrak{se}(2)^{\ast}$. - The blue curve is the trajectory of the costate $(h_1(t), h_2(t), - h_3(t))$ under the Lie–Poisson flow generated by the SR Hamiltonian - $\mathcal H = \tfrac12 (h_1^2 + h_2^2)$ — derived in Appendix A3, used - without proof in Part 2 §1. The trajectory winds around the cylinder - at constant $h_3$, with angular rate $-h_3$ (slide $h_3$ to vary). This - is a symplectic structure visualised: cylinders for non-trivial orbits, - pinched-off points along the axis $h_1 = h_2 = 0$ for the degenerate - orbits. + The slider $c$ sets the cylinder radius; $h_3$ runs along its axis. + Every Lie–Poisson flow — in particular the SR Hamiltonian flow of + Part 2, derived in Appendix A3 — is confined to one such cylinder, + because the Casimir $h_1^2 + h_2^2$ is conserved. The blue curve shows + a sample trajectory riding the chosen orbit. This is a symplectic + structure visualised: cylinders for the non-trivial orbits, pinched-off + points along the axis $h_1 = h_2 = 0$ for the degenerate ones.
@@ -488,10 +503,10 @@ points to remember: $X + Y$ pieces cancel, leaving $\varepsilon^2 [X, Y]$ at the leading surviving order. See Figure A1.2. - Even when $[X, Y] \neq 0$, the higher commutators - $[X, [X, Y]], [Y, [X, Y]]$ may vanish — and for $\mathrm{SE}(2)$, the - algebra is "step-2 nilpotent at infinity" in a sense, meaning many - identities truncate quickly. This is what makes the Sachkov closed forms - in Part 2 manageable. + $[X, [X, Y]], [Y, [X, Y]]$ are forced back into a small set of + directions. $\mathfrak{se}(2)$ is not nilpotent, but its bracket table is + simple enough that many BCH-type expansions truncate quickly in practice + — this is what keeps the Sachkov closed forms in Part 2 manageable. ## Connection to the elliptic project @@ -554,7 +569,7 @@ on them. Appendix A2 will use this language to give *the right* definition of a contact structure and prove Chow–Rashevskii. Appendix A3 will derive the Lie–Poisson equations $\dot h_1 = h_2 h_3, \dot h_2 = -h_1 h_3, \dot h_3 = -0$ — the equations Part 2 §1 asserts without proof — directly from the +-h_1 h_2$ — the equations Part 2 §1 asserts without proof — directly from the Pontryagin Maximum Principle. diff --git a/_posts/2026-05-02-geometry-of-seeing-A2-distributions-contact.md b/_posts/2026-05-02-geometry-of-seeing-A2-distributions-contact.md index 70f4007dd4c..4d7216b82d5 100644 --- a/_posts/2026-05-02-geometry-of-seeing-A2-distributions-contact.md +++ b/_posts/2026-05-02-geometry-of-seeing-A2-distributions-contact.md @@ -150,7 +150,7 @@ then any two points of $M$ can be joined by a piecewise-horizontal curve. **Sketch.** Bracket-generating means iterated brackets of the frame fields span $T_p M$. Concatenating short flows along a frame field $X_i$ for -times $\pm \varepsilon$ in the pattern of A1 Figure 1.2 produces a net +times $\pm \varepsilon$ in the pattern of Appendix A1's Figure A1.2 produces a net displacement of order $\varepsilon^2$ in the bracket direction, of order $\varepsilon^3$ in iterated-bracket directions, etc. Show that the smooth map $\mathbb R^n \to M$, $(t_1, \ldots, t_n) \mapsto \Phi^{X_{i_1}}_{t_1} @@ -159,7 +159,7 @@ origin under bracket-generation, hence is locally surjective by the inverse function theorem. Compose enough hops and you can reach any point. **For SE(2) the Hörmander condition is satisfied at depth 1.** We computed -$[X_1, X_2] = \pm X_3$, and $X_1, X_2, X_3$ already span $T_g\mathrm{SE}(2)$ +$[X_1, X_2] = -X_3$, and $X_1, X_2, X_3$ already span $T_g\mathrm{SE}(2)$ (three linearly independent fields). No deeper brackets are needed. This is the *minimal* possible depth and is what makes the SR Carnot–Carathéodory distance well-behaved on $\mathrm{SE}(2)$. @@ -229,7 +229,7 @@ $$\alpha \wedge d\alpha \;=\; (\sin\theta\,dx - \cos\theta\,dy) \wedge (\cos\theta\,d\theta \wedge dx + \sin\theta\,d\theta \wedge dy)$$ $$\quad - \;=\; -\sin^2\theta\,dx \wedge d\theta \wedge dy + \;=\; \sin^2\theta\,dx \wedge d\theta \wedge dy - \cos^2\theta\,dy \wedge d\theta \wedge dx \;=\; -dx \wedge dy \wedge d\theta \;\neq\; 0.$$ @@ -325,16 +325,16 @@ admissible $\gamma$ writes as $\dot\gamma = u_1(t) X_1 + u_2(t) X_2$ and $$L_{\mathrm{SR}}(\gamma) \;=\; \int_0^T \sqrt{u_1^2 + u_2^2}\,dt.$$ -After the $u_1 = 1$ unit-speed reduction, this collapses to the elastica +After the $u_1 = 1$ reduction, $s$ is arc length and $\kappa = \dot\theta$, +and the SR extremals project onto the critical curves of the elastica functional $\int \kappa^2(s)\,ds$ — exactly the integrand Euler minimised. Appendix A3 takes this and runs it through the Pontryagin Maximum Principle. -**Mitchell's compactness theorem** (1985) shows that on a connected -bracket-generating SR manifold the infimum is attained: a length-minimising -geodesic exists between any two points. This is the SR analogue of -Hopf–Rinow. For $\mathrm{SE}(2)$ it means every pair of V1 neurons is -connected by an *actual* shortest horizontal curve — not just an -approachable one. +The **sub-Riemannian Hopf–Rinow theorem** shows that on a connected, +complete bracket-generating SR manifold the infimum is attained: a +length-minimising geodesic exists between any two points. For +$\mathrm{SE}(2)$ it means every pair of V1 neurons is connected by an +*actual* shortest horizontal curve — not just an approachable one. ## Connection to the elliptic project @@ -369,28 +369,30 @@ def vf_bracket(X, Y, vars_): out[i] = sp.simplify(s) return out -X3 = vf_bracket(X1, X2, [x, y, th]) -print("[X1, X2] =", X3.T) -# → [-sin(θ), cos(θ), 0] (this is -X3 from the text; sign convention) +br12 = vf_bracket(X1, X2, [x, y, th]) +print("[X1, X2] =", br12.T) +# → [sin(θ), -cos(θ), 0] (this is -X3 from the text) -# In span{X1, X2}? Solve a·X1 + b·X2 = X3 for constants a, b +# In span{X1, X2}? Solve a·X1 + b·X2 = [X1, X2] for constants a, b a, b = sp.symbols('a b') -sol = sp.solve(a * X1 + b * X2 - X3, [a, b], dict=True) -print("solution (None means not in span):", sol) +sol = sp.solve(a * X1 + b * X2 - br12, [a, b], dict=True) +print("solution ([] means not in span):", sol) # → [] — no solution; distribution is NOT integrable ``` ```python -# Contact-form check: α ∧ dα ≠ 0 for SE(2) -alpha = [sp.sin(th), -sp.cos(th), 0] # α = sin θ dx - cos θ dy as coefficients -# d α: coefficients of dx∧dy, dx∧dθ, dy∧dθ -# dα = cos θ dθ ∧ dx + sin θ dθ ∧ dy (verify by direct exterior diff) -# α ∧ dα picks out the volume form coefficient on dx ∧ dy ∧ dθ: -vol_coeff = sp.simplify( - alpha[0] * sp.cos(th) * (-1) # sin θ · cos θ · (dx ∧ dθ ∧ dx) - vanishes; only the dy∧dθ term survives ↘ - + (-alpha[1]) * sp.sin(th) -) -# → 1, i.e. α ∧ dα = -dx ∧ dy ∧ dθ ≠ 0 (SE(2) is contact) +# Contact-form check: α ∧ dα ≠ 0 for SE(2). +# For a 1-form α = αx dx + αy dy + αθ dθ on a 3-manifold, +# α ∧ dα = (α · curl α) · dx ∧ dy ∧ dθ, +# so the manifold is contact iff α · curl α never vanishes. +ax, ay, ath = sp.sin(th), -sp.cos(th), 0 # α = sin θ dx − cos θ dy + +curl = (sp.diff(ath, y) - sp.diff(ay, th), # ∂y αθ − ∂θ αy + sp.diff(ax, th) - sp.diff(ath, x), # ∂θ αx − ∂x αθ + sp.diff(ay, x) - sp.diff(ax, y)) # ∂x αy − ∂y αx +vol_coeff = sp.simplify(ax * curl[0] + ay * curl[1] + ath * curl[2]) +print("α ∧ dα coefficient:", vol_coeff) +# → -1, i.e. α ∧ dα = −dx ∧ dy ∧ dθ ≠ 0 (SE(2) is contact) ``` ## What we covered, and what comes next diff --git a/_posts/2026-05-03-geometry-of-seeing-A3-pmp.md b/_posts/2026-05-03-geometry-of-seeing-A3-pmp.md index 3de49a8d436..938bd462a0f 100644 --- a/_posts/2026-05-03-geometry-of-seeing-A3-pmp.md +++ b/_posts/2026-05-03-geometry-of-seeing-A3-pmp.md @@ -5,9 +5,9 @@ title: "Appendix A3 — Calculus of Variations and the Pontryagin Maximum Princi subtitle: > From Euler–Lagrange to the PMP, then Lie–Poisson reduction on $\mathfrak{se}(2)^{\ast}$. Why the equations $\dot h_1 = h_2 h_3$, - $\dot h_2 = -h_1 h_3$, $\dot h_3 = 0$ that Part 2 §1 used as a starting - point are exactly what you get when you do the optimal-control problem - carefully on a Lie group. + $\dot h_2 = -h_1 h_3$, $\dot h_3 = -h_1 h_2$ that Part 2 §1 used as a + starting point are exactly what you get when you do the optimal-control + problem carefully on a Lie group. date: 2026-05-03 09:00:00 categories: [mathematics] tags: [sub-riemannian, SE2, optimal-control, pmp, hamiltonian, lie-groups] @@ -35,7 +35,7 @@ Part 2 §1 begins: "The Pontryagin Maximum Principle introduces a covector $\lambda$ in the cotangent bundle... (After several lines of unjustified algebra)... the Hamiltonian equations on $\mathfrak{se}(2)^{\ast}$ read -$\dot h_1 = h_2 h_3, \dot h_2 = -h_1 h_3, \dot h_3 = 0$." +$\dot h_1 = h_2 h_3, \dot h_2 = -h_1 h_3, \dot h_3 = -h_1 h_2$." This appendix supplies the missing derivation. It assumes Appendix A1 @@ -282,22 +282,30 @@ Lie–Poisson flow on $\mathfrak{se}(2)^{*}$. The pendulum angle $\varphi = 2\alpha$ is twice the phase of $(h_1, h_2)$ on the unit circle; the pendulum energy $E$ is fixed by the Casimir on the coadjoint orbit. -The reconstruction equation then ties the planar curve's heading $\theta$ -back to the costate. With $u_1^{\ast} = h_1 = \sin(\varphi/2)$ and +The reconstruction equation ties the planar curve's heading $\theta$ back to +the costate. With $u_1^{\ast} = h_1 = \sin(\varphi/2)$ and $u_2^{\ast} = h_2 = \cos(\varphi/2)$, the SE(2) ODE $(\dot x, \dot y, \dot\theta) = (u_1^{\ast}\cos\theta, u_1^{\ast}\sin\theta, u_2^{\ast})$ -gets reparametrised by Euclidean arc length $s$ (with $ds = |u_1^{\ast}|\,dt$); -the resulting projected curve has curvature - -$$\kappa(s) = \frac{u_2^{\ast}}{u_1^{\ast}} = \cot(\varphi/2).$$ - -Tracking $\kappa$ through the three pendulum regimes recovers the three -elastica families of Part 2 §3 — inflectional ($E < 1$, $\kappa = 2k\,\mathrm{cn}$), -separatrix ($E = 1$, $\kappa = 2\,\mathrm{sech}$), non-inflectional -($E > 1$, $\kappa = 2\,\mathrm{dn}$). Appendix A4 §3 tracks that -substitution in detail; the **elastica curvature ODE** +reparametrised by Euclidean arc length $s$ (with $ds = |u_1^{\ast}|\,dt$) +gives the projected curve the curvature + +$$\kappa_{\mathrm{SR}}(s) = \frac{u_2^{\ast}}{u_1^{\ast}} = \cot(\varphi/2).$$ + +This blows up wherever $u_1^{\ast} = 0$: the sub-Riemannian projections are +allowed to have **cusps** — the reversal points of the parallel-parking +trajectories of Appendix A2. + +The smooth elastica that Part 2 plots — $\kappa = 2k\,\mathrm{cn}$, +$2\,\mathrm{sech}$, $2\,\mathrm{dn}$ — are the **Euler elastic problem**: the +sister problem in which the curve carries *unit forward speed* and the +*heading* $\theta$ itself is the pendulum variable, so $\kappa = \dot\theta$ +stays bounded. Both problems are driven by the *same* pendulum equation +$\ddot\varphi + \sin\varphi = 0$ — that shared vertical subsystem is the real +content of the reduction — but their projected curves differ, and it is the +elastic representatives the figures draw. Appendix A4 §3 carries the elliptic +substitution through; the **elastica curvature ODE** $\kappa''(s) + \tfrac12\kappa^3 - \mu\kappa = 0$ is the Duffing form -equivalent to the pendulum and is what the figures in Part 2 actually plot. +equivalent to the pendulum. @@ -355,11 +363,12 @@ $(x, y, \theta)$ chart this is exactly Part 1's Frenet–Serret integration: $$\dot x = u_1^{\ast} \cos\theta, \qquad \dot y = u_1^{\ast} \sin\theta, \qquad \dot \theta = u_2^{\ast}.$$ -Setting $u_1^{\ast} = 1$ (unit-speed normalisation), we get $u_2^{\ast} = \dot\theta = -\kappa$, which is the **curvature** of the projected plane curve. The -relation $\kappa = h_2 / \sqrt c$ ties the costate to the curvature -directly: $\kappa$ inherits the $h_2$-pendulum dynamics and so becomes a -Jacobi sn (libration), sech (separatrix), or dn (rotation) — Part 2 §2. +In the Euler elastic problem one instead fixes $u_1^{\ast} = 1$ (unit forward +speed); then $s = t$, and $u_2^{\ast} = \dot\theta = \kappa$ is the +**curvature** of the projected plane curve. Here the heading $\theta$ is +itself the pendulum, so its curvature is the pendulum velocity — a Jacobi +cn (libration), sech (separatrix), or dn (rotation), i.e. the +$\kappa = 2k\,\mathrm{cn}(s\mid k^2)$ family of Part 2 §2. @@ -382,15 +391,16 @@ Jacobi sn (libration), sech (separatrix), or dn (rotation) — Part 2 §2.
- Figure A3.3. Left: the costate $(h_1, h_2, h_3)$ winds - on a cylinder $h_1^2 + h_2^2 = c$ at the constant rate - $\dot\phi = -h_3$. Right: the plane projection of the resulting - SE(2) geodesic, integrated by the reconstruction equation $\dot g = - g\cdot\xi(t)$. Vary $h_3$ and the curve interpolates between - near-circular (small $h_3$) and the elastica regime; vary $\sqrt c$ - and the curve scales without changing shape — confirming Part 2's - observation that only the dimensionless ratio $h_3 / \sqrt c$ - determines which elastica family you land in. + Figure A3.3. A schematic of the reconstruction step, + drawn in the constant-$h_3$ approximation. Left: a curve on the + cylinder $h_1^2 + h_2^2 = c$, with $(h_1, h_2)$ circling at rate + $-h_3$. Right: the plane curve obtained by feeding that costate into + the reconstruction equation $\dot g = g\cdot\xi(t)$. Vary $h_3$ and the + curve interpolates between near-circular (small $h_3$) and the elastica + regime; vary $\sqrt c$ and it scales without changing shape — the curve + family depends only on the dimensionless ratio $h_3 / \sqrt c$. (In the + full flow $h_3$ varies too, by $\dot h_3 = -h_1 h_2$; freezing it keeps + this picture readable.)
@@ -416,32 +426,28 @@ map). ## Code ```python -# Verify the Lie–Poisson equations on se(2)* numerically -# from a random initial costate (h1, h2, h3) and check that: -# - h3(t) is constant -# - h1²+h2² is constant -# - h1(t) = sqrt(c) cos(-h3 t + φ0) +# Verify the Lie–Poisson equations on se(2)* numerically. +# Conserved quantities: H = ½(h1²+h2²) and Casimir C = h1²+h3². +# Note: h3 is NOT conserved — the third equation is dh3/dt = -h1·h2. import numpy as np from scipy.integrate import solve_ivp def lie_poisson_se2(t, h): h1, h2, h3 = h - return [h2*h3, -h1*h3, 0.0] + return [h2*h3, -h1*h3, -h1*h2] h0 = [0.6, 0.4, 0.7] sol = solve_ivp(lie_poisson_se2, [0, 8], h0, rtol=1e-10, atol=1e-12, t_eval=np.linspace(0, 8, 400)) -# Check Casimirs -c = sol.y[0]**2 + sol.y[1]**2 -print(f"max |c - c0| / |c0| = {np.max(np.abs(c - c[0]))/c[0]:.2e}") # ~ 1e-10 -print(f"max |h3 - h30| = {np.max(np.abs(sol.y[2] - h0[2])):.2e}") # ~ 1e-12 +h1, h2, h3 = sol.y +H = 0.5 * (h1**2 + h2**2) # Hamiltonian +C = h1**2 + h3**2 # Casimir +print(f"max |H - H0| / H0 = {np.max(np.abs(H - H[0]))/H[0]:.2e}") # ~ 1e-10 +print(f"max |C - C0| / C0 = {np.max(np.abs(C - C[0]))/C[0]:.2e}") # ~ 1e-10 -# Match the closed form -phi0 = np.arctan2(h0[1], h0[0]) -phi_t = phi0 - h0[2] * sol.t -h1_pred = np.sqrt(c[0]) * np.cos(phi_t) -print(f"max |h1 - prediction| = {np.max(np.abs(sol.y[0] - h1_pred)):.2e}") # ~ 1e-10 +# h3 itself drifts — confirm it is genuinely not conserved +print(f"h3 range = [{h3.min():.3f}, {h3.max():.3f}] h3(0) = {h0[2]}") ``` ```python @@ -473,15 +479,14 @@ length problem, maximisation over $u_1, u_2$ on the unit circle gives the normal Hamiltonian $\mathcal H_n = \tfrac12(h_1^2 + h_2^2)$. Lie–Poisson reduction on $\mathfrak{se}(2)^{\ast}$ collapses the $T^{\ast}\mathrm{SE}(2)$ flow to the costate equations -$\dot h_1 = h_2 h_3, \dot h_2 = -h_1 h_3, \dot h_3 = 0$ — exactly the -equations Part 2 §1 wrote down. The substitution $h_1 = \sqrt c \cos\phi, -h_2 = \sqrt c \sin\phi$ turns the costate into a uniformly-rotating phase, -and differentiating once more gives the nonlinear pendulum equation for -the curvature. +$\dot h_1 = h_2 h_3, \dot h_2 = -h_1 h_3, \dot h_3 = -h_1 h_2$ — exactly the +equations Part 2 §1 wrote down. Writing $h_1 = \sin\alpha$, $h_2 = \cos\alpha$ +on the unit Hamiltonian level and differentiating $\varphi = 2\alpha$ once +more gives the nonlinear pendulum equation $\ddot\varphi + \sin\varphi = 0$. Appendix A4 will solve the pendulum equation in closed form using Jacobi elliptic functions and the AGM, recovering the period $4K(k^2)$ and the -explicit $\kappa(s) = 2k\,\mathrm{sn}(s\mid k^2)$ formula of Part 2. +explicit $\kappa(s) = 2k\,\mathrm{cn}(s\mid k^2)$ formula of Part 2. diff --git a/_posts/2026-05-04-geometry-of-seeing-A4-jacobi-elliptic.md b/_posts/2026-05-04-geometry-of-seeing-A4-jacobi-elliptic.md index 65ff226de0d..9629e45827a 100644 --- a/_posts/2026-05-04-geometry-of-seeing-A4-jacobi-elliptic.md +++ b/_posts/2026-05-04-geometry-of-seeing-A4-jacobi-elliptic.md @@ -180,11 +180,8 @@ so $\ddot\varphi + \sin\varphi = 0$. ✓ This is **the** identity behind the elastica: the pendulum solution is literally the Jacobi-am function, and the curvature -$\kappa(s) = \dot\varphi(s) = 2k\,\mathrm{cn}(s)$ is one Jacobi function's -worth — exactly the formula Part 2 uses. (An alternative convention takes -$\kappa = 2k\,\mathrm{sn}$, related by a quarter-period shift $s \to s + K$; -either describes the same family of geodesics, just starting at a different -arc-length offset.) +$\kappa(s) = \dot\varphi(s) = 2k\,\mathrm{cn}(s\mid k^2)$ is one Jacobi +function's worth — exactly the formula Part 2 uses. ## The complete elliptic integrals @@ -228,8 +225,7 @@ Gauss (1799, unpublished) proved the miracle: $$K(m) \;=\; \frac{\pi}{2\,\mathrm{AGM}\!\bigl(1,\;\sqrt{1 - m}\bigr)}.$$ -Equivalently $K(m) = \pi / (2\,\mathrm{AGM}(\sqrt{1+\sqrt{1-m}}, -\sqrt{1-\sqrt{1-m}}))/\sqrt 2$ via the descending Landen transformation. +Six AGM steps from $(1, \sqrt{1-m})$ pin $K(m)$ down to full double precision. This is what `ellipticK` in the Python `elliptic` package computes — and @@ -302,7 +298,7 @@ under the rescaling. energy $E$, blue solid; logarithmic divergence at $E \to 1^-$ (red dashed). At small amplitude $E \to -1$, $T \to 2\pi$ (the harmonic limit); at $E = 0$, $T \approx 7.42$, already noticeably longer than - $2\pi$. Asymptotic prediction $T \sim 2\log(16/(1-E))$ near + $2\pi$. Asymptotic prediction $T \sim 2\log(32/(1-E))$ near separatrix overlaid (grey dotted). This is the same plot drawn by the elliptic project's physical-pendulum example — modulo axis labels, the @@ -314,20 +310,20 @@ under the rescaling. ## Identities used in Part 2 §4 -Part 2 uses the closed-form integral +Part 2 uses the closed-form heading integral $$\theta(s) \;=\; \theta_0 + 2\arcsin\bigl(k\,\mathrm{sn}(s\mid k^2)\bigr).$$ -Differentiating with respect to $s$: -$\dot\theta = 2k\,\mathrm{cn}\,\mathrm{dn} / \sqrt{1 - k^2 \mathrm{sn}^2} - = 2k\,\mathrm{cn}\,\mathrm{dn} / \mathrm{dn} - = 2k\,\mathrm{cn}$. -But Part 2 writes $\kappa = 2k\,\mathrm{sn}$. The reconciliation: there -are *two* parametrisations of the inflectional family by Jacobi functions, -related by a quarter-period shift $s \to s + K(k^2)$. Under that shift -$\mathrm{sn}(s + K) = \mathrm{cn}(s) / \mathrm{dn}(s)$ and the two -conventions translate. Both the Sachkov closed form (with sn) and the -"angle-of-pendulum" form (with cn) are used in the literature. +Differentiating with respect to $s$, and using +$\sqrt{1 - k^2\,\mathrm{sn}^2} = \mathrm{dn}$, + +$$\kappa = \dot\theta + = \frac{2k\,\mathrm{cn}\,\mathrm{dn}}{\sqrt{1 - k^2\,\mathrm{sn}^2}} + = \frac{2k\,\mathrm{cn}\,\mathrm{dn}}{\mathrm{dn}} + = 2k\,\mathrm{cn}(s\mid k^2),$$ + +exactly the boxed curvature of Part 2 §3 — the heading integral and the +curvature formula are one statement, differentiated once. The plane curve integration uses the **second-kind incomplete integral** @@ -380,7 +376,7 @@ def pend(t, y): T_numeric = [] for E in E_vals: phi0 = 0.0 - phidot0 = np.sqrt(2 * (E - np.cos(phi0))) # initial velocity from E + phidot0 = np.sqrt(2 * (E + np.cos(phi0))) # from E = ½φ̇² − cos φ sol = solve_ivp(pend, [0, 30], [phi0, phidot0], rtol=1e-12, atol=1e-14, dense_output=True) # Find first return to phi = 0 with phidot > 0 @@ -644,8 +640,8 @@ function drawPeriod() { .attr('font-family', 'Source Sans 3').attr('font-size', 11).attr('fill', '#888') .text('T = 2π (harmonic)'); - // Asymptotic prediction T ~ 2 log(16/(1-E)) at E ~ 1 - const asy = Earr.filter(E => E > 0).map(E => ({ E, T: 2 * Math.log(16 / Math.max(1e-3, 1 - E)) })); + // Asymptotic prediction T ~ 2 log(32/(1-E)) at E ~ 1 + const asy = Earr.filter(E => E > 0).map(E => ({ E, T: 2 * Math.log(32 / Math.max(1e-3, 1 - E)) })); g.append('path') .attr('d', d3.line().x(p => xS(p.E)).y(p => yS(Math.min(p.T, 60)))(asy)) .attr('fill', 'none').attr('stroke', '#888').attr('stroke-width', 1).attr('stroke-dasharray', '2,3'); diff --git a/_posts/2026-05-05-geometry-of-seeing-A5-sr-exponential.md b/_posts/2026-05-05-geometry-of-seeing-A5-sr-exponential.md index 300961ef053..9f79714aaac 100644 --- a/_posts/2026-05-05-geometry-of-seeing-A5-sr-exponential.md +++ b/_posts/2026-05-05-geometry-of-seeing-A5-sr-exponential.md @@ -71,8 +71,8 @@ Two parametrisations of the same costate space are convenient: - **Cylinder coordinates** $(c, \omega_0, \phi_0)$: $h_1 = \sqrt c \cos\phi_0, h_2 = \sqrt c \sin\phi_0, h_3 = \omega_0$. By rescaling, set $c = 1$ on the energy surface. -- **Pendulum-energy coordinates** $(E, \phi_0)$: $E = \tfrac12 \omega_0^2 - - \cos\phi_0$ (with the rescaling of A3 §6). Here $E < 1$ is libration +- **Pendulum-energy coordinates** $(E, \phi_0)$: $E = 2C - 1$ with the + Casimir $C = h_1^2 + h_3^2$ (Appendix A3). Here $E < 1$ is libration (inflectional), $E = 1$ is the separatrix (Euler spiral), $E > 1$ is rotation (non-inflectional). @@ -148,8 +148,8 @@ Three takeaways: Figure 4 — and indeed this figure is the same one, lifted to a more controllable form. As $k \to 1^-$ the inflectional family's period $4K(k^2)$ diverges (Appendix A4) and the curve spirals; for - $k > 1$ (non-inflectional) the curve is a deformed circle that - closes after $T = 2\sqrt{m}\,K(m)$. The endpoint dot is + $k > 1$ (non-inflectional) the curvature is one-signed with spatial + period $T = 2K(m)$ (a closed circle only in the $m \to 0$ limit). The endpoint dot is $\mathrm{Exp}_T(\mu_0)$ for the specified $\mu_0$ and $T$. @@ -216,8 +216,10 @@ For sufficiently symmetric SR problems (and SE(2) is one of them), $$\boxed{\;t_{\mathrm{cut}} \;=\; t_{\mathrm{Maxwell}}^{(1)},\;}$$ -i.e. the cut locus equals the closure of the Maxwell locus. This is the -content of Moiseev–Sachkov (2010, arXiv:0807.4731) for $\mathrm{SE}(2)$. +i.e. the cut locus equals the closure of the Maxwell locus. Characterising +the Maxwell strata for $\mathrm{SE}(2)$ is the work of Moiseev–Sachkov +(2010, arXiv:0807.4731); pinning the cut locus down exactly is Sachkov +(2011, arXiv:0903.0727). For the inflectional family with modulus $k$, the first Maxwell time is exactly @@ -244,7 +246,7 @@ Part 3. 5.00 - γ_A: κ = +2k·sn — γ_B: κ = −2k·sn — coincidence requires y_A(s) = 0 + γ_A: κ = +2k·cn — γ_B: κ = −2k·cn — coincidence requires y_A(s) = 0 @@ -252,7 +254,7 @@ Part 3.
Figure A5.2. The σ-symmetric pair: two geodesics $\gamma_A, \gamma_B$ leaving the origin with curvatures - $\pm 2k\,\mathrm{sn}(s\mid k^2)$ respectively. By the + $\pm 2k\,\mathrm{cn}(s\mid k^2)$ respectively. By the $y \to -y$ reflection symmetry of the pendulum equation, $\gamma_B(s) = (x_A(s),\, -y_A(s),\, -\theta_A(s))$ — they trace mirror-image curves. As SE(2) configurations, they coincide exactly @@ -260,11 +262,15 @@ Part 3. The right panel plots $|\gamma_A(s) - \gamma_B(s)|$ over $s \in [0, T]$; its first zero crossing (red marker, if any) is the first Maxwell time of this pair. - For "lemniscate" values of $k$ (where the elastica closes into a - figure-eight — numerically $k_c \approx 0.84$) this first zero - coincides with the curvature period $4K(k^2)$. For generic $k$ the - Maxwell time is k-dependent and given by Sachkov's transcendental - equation in the modulus. Slide $k$ to find lemniscate values. + Because $y_A(s) = 2k\bigl(1 - \mathrm{cn}(s\mid k^2)\bigr) \ge 0$ returns to + zero only at $s = 4K(k^2)$ (and its multiples), the pair first re-coincides + in full — position and heading — at $s = 4K(k^2)$, independent of + $k$: exactly the first Maxwell time $t_{\mathrm{Maxwell}}^{(1)} = 4K(k^2)/\omega_0$ + boxed above. At the special "figure-eight" modulus + $k_c \approx 0.909$ (root of $2E(k^2) = K(k^2)$) that shared endpoint sits + back at the origin, so the closed curve is itself a single self-crossing + lemniscate; for other $k$ the two curves still meet at $s = 4K(k^2)$, just + away from the origin.
@@ -274,16 +280,17 @@ Part 3. Fix $T$ and vary the initial costate over a 1-parameter slice of the unit- energy surface — concretely, signed initial curvature $k \in [-0.95, 0.95]$ -in $\kappa(s) = 2k\,\mathrm{sn}(s\mid k^2)$. Each $k$ launches a distinct +in $\kappa(s) = 2k\,\mathrm{cn}(s\mid k^2)$. Each $k$ launches a distinct geodesic from the origin. The set of *positions* reached at exact arc length $T$ — one position per geodesic — is the **wavefront** at time $T$: $$\mathcal W_T \;:=\; \bigl\{\,(x(T;k),\; y(T;k)) : k \in [-1, 1]\,\bigr\} \;\subset\; \mathbb R^2.$$ It is a *continuous curve* in the plane (because $k \mapsto $ trajectory is -continuous), and as $T$ grows it sweeps outward. At small $T$ the -wavefront is approximately a circle of radius $T$ around the origin -(everything moves at unit speed in approximately straight lines). As $T$ +continuous), and as $T$ grows it sweeps outward. At small $T$ — since every +geodesic leaves the origin heading the same way and curves only gently — +the wavefront is a short, almost-straight arc near $(T, 0)$, transverse to +the launch direction. As $T$ approaches the first Maxwell time $T_{\mathrm M} = 4K(k^2)/\omega_0$, neighbouring trajectories begin to converge and the wavefront develops **cusps** — these are the projections of conjugate points, where @@ -327,12 +334,13 @@ be locally surjective along a critical curve. $\mathcal W_T$ moves outward and reshapes as time grows. Slide $T$ forward and watch:
    -
  1. at small $T$ the wavefront is nearly a circle;
    2. +
  2. at small $T$ the wavefront is a short arc near $(T, 0)$;
  3. at $T \approx \pi$ it lengthens and starts to flatten;
    4. -
  4. at $T \approx 4K(k_c^2) \approx 6.5$–$7$ cusps appear at - the four corners of the wavefront — these are the first conjugate - points (marked with red rings);
    5. -
  5. past the cusps, the wavefront self-intersects: those crossings +
  6. around $T \approx 2\pi$ — the smallest period $4K(0)$ in the swept + family (the near-straight $k \to 0$ geodesics) — the first + cusps appear at the corners of the wavefront: these are the + first conjugate points (marked with red rings);
    7. +
  7. as $T$ grows the wavefront self-intersects: those crossings are the Maxwell stratum drawn in Figure A5.2.
Press play to animate $T$ continuously. The four-fold @@ -378,10 +386,8 @@ def inflectional_geodesic(k, omega0, phi0, T, N=1500): m = k * k s = np.linspace(0, T, N) sn, cn, dn = ellipj(s, m) # vectorised - # Pendulum: sin(φ/2) = k·sn(s|m), so: - half_phi = np.arcsin(k * sn) - phi = 2 * half_phi - # Curvature κ = dθ/ds = ω0 sin(φ) — but normalise + # Heading half-angle: sin(θ/2) = k·sn(s|m), so θ = 2·arcsin(k·sn). + # Curvature κ = dθ/ds = 2k·cn·dn / sqrt(1 − k²sn²) = 2k·cn (since dn = sqrt). kappa = 2 * k * cn * dn / np.sqrt(1 - k * k * sn * sn) theta = np.zeros_like(s) x = np.zeros_like(s) @@ -508,7 +514,7 @@ fluently, and Parts 3 and 4 will be approachable when they ship. // ── shared helpers ───────────────────────────────────────────────────── function inflectionalGeodesic(k, T, N) { - // Curvature κ(s) = 2k·sn(s|k²) (Part 2 convention) + // Curvature κ(s) = 2k·cn(s|k²) (Part 2 convention) // dθ/ds = κ, dx/ds = cos θ, dy/ds = sin θ const m = k * k; const ds = T / N; @@ -517,7 +523,7 @@ function inflectionalGeodesic(k, T, N) { for (let i = 0; i < N; i++) { const s = i * ds; const j = ellipj(s + ds / 2, m); - const k1 = 2 * k * j.sn; + const k1 = 2 * k * j.cn; const tmid = theta + 0.5 * k1 * ds; x += Math.cos(tmid) * ds; y += Math.sin(tmid) * ds; @@ -633,7 +639,7 @@ function drawExpMap() { } else { const m = Math.max(0.01, Math.min(0.99, 2 - kVal)); const Km = ellipticK(m); - info = `non-inflectional m = ${m.toFixed(2)} period 2√m·K = ${(2*Math.sqrt(m)*Km).toFixed(2)} T = ${T.toFixed(2)}`; + info = `non-inflectional m = ${m.toFixed(2)} period 2K(m) = ${(2*Km).toFixed(2)} T = ${T.toFixed(2)}`; } g.append('text').attr('x', margin.l).attr('y', margin.t + 12) .attr('font-family', 'JetBrains Mono').attr('font-size', 11).attr('fill', '#555') @@ -674,7 +680,7 @@ function drawMaxwell() { for (let i = 0; i < N; i++) { const s = i * ds; const j = ellipj(s + ds / 2, m); - const k1 = 2 * k * j.sn; + const k1 = 2 * k * j.cn; const tMid = theta + 0.5 * k1 * ds; x += Math.cos(tMid) * ds; y += Math.sin(tMid) * ds; @@ -691,17 +697,11 @@ function drawMaxwell() { // Per-s separation: |γA(s) − γB(s)| = 2|y_A(s)| const sep = ptsA.map(p => ({ s: p.s, d: 2 * Math.abs(p.y) })); - // Find first zero crossing of y_A after s = ds (skip the initial s = 0) - let firstZero = null; - for (let i = 8; i < ptsA.length; i++) { - if ((ptsA[i - 1].y > 0) !== (ptsA[i].y > 0)) { - // Linear interpolation - const y0 = ptsA[i - 1].y, y1 = ptsA[i].y; - const t0 = ptsA[i - 1].s, t1 = ptsA[i].s; - firstZero = t0 + (t1 - t0) * (-y0) / (y1 - y0); - break; - } - } + // The pair coincides in full (position AND heading) at s = 4K(k²): there + // y_A = 2k(1 − cn) returns to 0 and θ_A returns to 0 (mod 2π). Since + // y_A ≥ 0 only *touches* zero (never crosses), a sign-change search would + // miss it — use the analytic first Maxwell time directly. + const firstZero = (4 * Km <= T_full) ? 4 * Km : null; // Layout: left panel = trajectories up to T; right panel = sep(s) plot const split = W * 0.55; @@ -850,7 +850,7 @@ function buildConjTrajectories() { for (let i = 0; i < conjState.N; i++) { const s = i * ds; const j = ellipj(s + ds / 2, m); - const kappaMid = sign * 2 * ka * j.sn; + const kappaMid = sign * 2 * ka * j.cn; const tMid = theta + 0.5 * kappaMid * ds; x += Math.cos(tMid) * ds; y += Math.sin(tMid) * ds; @@ -994,11 +994,11 @@ function drawConjugate() { g.append('circle').attr('cx', xS(0)).attr('cy', yS(0)).attr('r', 4).attr('fill', '#222'); // Status text + scale legend - // Use k ≈ 0.85 as a reference for "first Maxwell time" - const Tmax_ref = 4 * ellipticK(0.85 * 0.85); + // Reference: the figure-eight modulus k_c ≈ 0.909 (root of 2E(k²) = K(k²)). + const Tmax_ref = 4 * ellipticK(0.9089 * 0.9089); g.append('text').attr('x', margin.l).attr('y', margin.t + 12) .attr('font-family', 'JetBrains Mono').attr('font-size', 11).attr('fill', '#555') - .text(`T = ${conjState.T.toFixed(2)} (1st Maxwell @ k = 0.85: T₁ = ${Tmax_ref.toFixed(2)})`); + .text(`T = ${conjState.T.toFixed(2)} (1st Maxwell @ k = 0.909: T₁ = ${Tmax_ref.toFixed(2)})`); // Colour-bar mini-legend const lx = margin.l + 4, ly = H - margin.b - 14; diff --git a/projects/geometry-of-seeing.md b/projects/geometry-of-seeing.md index 194a7697a38..3a6c773c156 100644 --- a/projects/geometry-of-seeing.md +++ b/projects/geometry-of-seeing.md @@ -165,11 +165,11 @@ elliptic cosine, with spatial period $T = 4K(k^2)$. - **Petitot's contact model** (Part 1): why V1 lifts the image to SE(2) and why modal completion is a geodesic problem. -- **Complete parametrisation** (Part 2): all three inflectional families, the - Euler separatrix, and the non-inflectional family, written in closed form using - $\mathrm{sn}, \mathrm{cn}, \mathrm{dn}$. +- **Complete parametrisation** (Part 2): all three families — the inflectional + family, the Euler separatrix, and the non-inflectional family — written in + closed form using $\mathrm{sn}, \mathrm{cn}, \mathrm{dn}$. - **Maxwell strata** (Part 3): discrete symmetry group $$\mathbb{Z}_2 \times \mathbb{Z}_2$$, - the first Maxwell time $$t_{\max}^1 = 2\pi/\sqrt{H}$$, loss of optimality. + the first Maxwell time $$t_{\max}^1 = 4K(k^2)/\omega_0$$, loss of optimality. - **The open problem** (Part 4): the cut time is bounded above by $$t_{\max}^1$$ but the exact value for $k \in (0,1)$ remains unproved.