The Length-Preserving Curvature Flow: the Andrews–Bryan Theorem, Numerically

MATH 8230 — Discrete Differential Geometry. Companion to the lecture notes Round limits of the length-preserving curvature flow: the Andrews–Bryan argument (lecture_notes/andrews-bryan-grayson.pdf), which follow Andrews & Bryan, J. reine angew. Math. 653 (2011) 179–187.

What this notebook does. The lecture notes prove that an embedded closed plane curve of length  evolving by the length-preserving curvature flow stays embedded and converges exponentially, in , to a unit circle. Nobody should read that proof before seeing the flow, so the notebook opens with the show: after the machinery is loaded, Section 3 runs five curves — dumbbell, yeti, octopus, Cthulhu, and a tight double spiral — as scrubbable movies on the flow's natural logarithmic clock. Then the mathematics: the theorem (Section 5), its comparison function verified symbolically (Section 6), the chord–arc scale made computable (Section 7), and the distance comparison rendered as a moving chord envelope squeezed between the circle profile and the sliding barrier (Section 8), with the curvature bound and the monotone scale following (Section 9). The closing sections compare the smooth world with the discrete one: the "obvious" polygonal translations of the flow, their exact structure, their convergence, and interactive demonstrations.

Evaluate Sections 2 and 3 first (a couple of minutes: five spectral solves plus their stage clocks); everything afterwards is fast, except where marked. The first cell below restarts the kernel, so a full evaluation always starts clean — stale definitions and stale memoized values from an earlier session are the most confusing bugs a notebook can have:

–Jason Cantarella and Claude (Fable 5)

Exit[];   (* fresh kernel: a full evaluation of this notebook always starts clean *)

1. The length-preserving flow

Let  be a family of closed curves,  arclength,  the unit tangent,  the inward unit normal for a positively oriented embedded curve,  the signed curvature, so  and a convex curve has . The object of study is the initial value problem

for a smooth embedded initial curve of length . The term  is curvature-driven smoothing; the nonlocal dilation  is tuned to freeze total length: the arclength element evolves by , so  and  throughout.

By Cauchy–Schwarz and the theorem of turning tangents ( for an embedded curve),

with equality iff , i.e. iff the curve is a unit circle. The unit circle centred at the origin is a stationary solution.

The flow is not translation-invariant. The dilation term uses the position vector, so the flow depends on the choice of origin. If  solves the flow, so does  for any  with  — and since , translations are amplified at least like . The geometric quantities the theorem is about (, , , ) never see the difference, but the images of two solutions with different origins fly apart exponentially. Section 4 shows this happening; it is not a numerical artifact.

Relation to curve shortening flow. If  solves CSF  and we set , then  with the reparametrized clock  solves the length-preserving flow: our flow is CSF viewed at constant magnification, in its intrinsic clock. CSF's finite extinction time  becomes . This dictionary is exact and reversible; Section 4 uses it as an independent numerical check, and it is how the Andrews–Bryan theorem proves Grayson's theorem.

The equation we actually solve.  is parabolic only modulo reparametrization. As with CSF, the standard gauge fix replaces  by , which differs from it by a tangential vector field, so it traces the same family of curves while being uniformly parabolic — and its tangential term conveniently pulls the parametrization toward constant speed.

2. Setup: machinery, the solver, and the menagerie

(* ------------------------------------------------------------------- *)
(*  A discrete curve is an n x 2 matrix of vertices, read cyclically:   *)
(*  gamma_0, ..., gamma_{n-1}, with gamma_n = gamma_0.                  *)
(* ------------------------------------------------------------------- *)

edges[P_] := RotateLeft[P] - P;                      (* e_i = gamma_{i+1} - gamma_i *)
edgeLengths[P_] := Sqrt[Total[edges[P]^2, {2}]];
tangents[P_] := edges[P]/edgeLengths[P];             (* unit tangent T_i of edge i *)
dualLengths[P_] := With[{h = edgeLengths[P]}, (RotateRight[h] + h)/2];   (* ell_i *)
polyLength[P_] := Total[edgeLengths[P]];
polyArea[P_] := Total[P[[All, 1]] RotateLeft[P][[All, 2]]
                    - P[[All, 2]] RotateLeft[P][[All, 1]]]/2;   (* shoelace; CCW > 0 *)

(* exterior (turning) angle theta_i at vertex i, in (-Pi, Pi] *)
turningAngles[P_] := Module[{T = tangents[P], Tp},
   Tp = RotateRight[T];
   ArcTan[Tp[[All, 1]] T[[All, 1]] + Tp[[All, 2]] T[[All, 2]],
          Tp[[All, 1]] T[[All, 2]] - Tp[[All, 2]] T[[All, 1]]]];

(* inward unit normal at a vertex: the angle bisector of T_{i-1}, T_i rotated by +90 *)
vertexNormals[P_] := Module[{T = tangents[P], s},
   s = RotateRight[T] + T;
   Transpose[{-s[[All, 2]], s[[All, 1]]}]/Sqrt[Total[s^2, {2}]]];

(* resample a closed polygon at n equally spaced arclength points *)
resample[P_, n_] := Module[{Q = Append[P, First[P]], s, L, f},
   s = Prepend[Accumulate@Sqrt[Total[(Rest[Q] - Most[Q])^2, {2}]], 0.];
   L = Last[s];
   f = Interpolation[Transpose[{s, Q}], InterpolationOrder -> 1];
   Table[f[k L/n], {k, 0, n - 1}]];

(* least-squares circle fit (Kasa): returns {center, radius} *)
circleFit[P_] := Module[{A, b, cx, cy, d},
   A = Transpose[{2 P[[All, 1]], 2 P[[All, 2]], ConstantArray[1., Length[P]]}];
   b = Total[P^2, {2}];
   {cx, cy, d} = LeastSquares[A, b];
   {{cx, cy}, Sqrt[d + cx^2 + cy^2]}];

(* scale-free measure of "how far from concyclic": 0 iff the vertices lie on a circle *)
circleResidual[P_] := Module[{c, r}, {c, r} = circleFit[P];
   StandardDeviation[Sqrt[Total[(# - c)^2]] & /@ P]/r];

(* exact distance from a point to a closed polyline, vectorized over segments *)
polylineDistance[Q_][p_] := Module[{A = Q, B = RotateLeft[Q], d, w, s, proj, pp},
   pp = ConstantArray[p, Length[A]];
   d = B - A; w = pp - A;
   s = Clip[Total[w d, {2}]/Total[d^2, {2}], {0., 1.}];
   proj = A + s d;
   Sqrt[Min[Total[(pp - proj)^2, {2}]]]];

drawPolygon[P_, opts : OptionsPattern[Graphics]] :=
  Graphics[{Line[Append[P, First[P]]], Red, PointSize[Medium], Point[P]},
   opts, AspectRatio -> Automatic];

(* sample a parametrized curve as a polygon with equally spaced arclength *)
sampleCurve[f_, n_, oversample_ : 20] :=
  resample[N@Table[f[u], {u, 0., 2 Pi (1 - 1/(oversample n)), 2 Pi/(oversample n)}], n];

(* rescale a parametrized curve to length exactly 2 Pi, measuring length on a fine
   polygon rather than by NIntegrate: accurate to ~10^-6 and free of convergence
   messages on curves whose speed has sharp minima (the dumbbell) *)
normalizeCurve[f_] := Module[{pts, L},
   pts = N@Table[f[u], {u, 0., 2 Pi (1 - 1/4096), 2 Pi/4096}];
   L = polyLength[pts];
   With[{r = 2 Pi/L}, Function[u, r f[u]]]];

The flow preserves length , so every initial curve is rescaled to that length once and for all. These are the same test curves as in the curve shortening flow notebook, now normalized.

circleCurve[u_] := {Cos[u], Sin[u]};
ellipseCurve = normalizeCurve[{1.6 Cos[#], 0.5 Sin[#]} &];
blobCurve = normalizeCurve[(1 + 0.30 Cos[3 #] + 0.18 Sin[5 #] - 0.12 Cos[2 #]) {Cos[#], Sin[#]} &];
lopsidedCurve = normalizeCurve[(1 + 0.35 Cos[#] + 0.25 Cos[3 # + 0.7] + 0.12 Sin[2 #]) {Cos[#], 1.15 Sin[#]} &];
dumbbellCurve = normalizeCurve[{1.5 Cos[#], 0.55 Sin[#] (1 + 0.9 Cos[#]^2)} &];

GraphicsRow[drawPolygon[sampleCurve[#, 60], ImageSize -> 165] & /@
  {circleCurve, ellipseCurve, blobCurve, lopsidedCurve, dumbbellCurve}]

And the yeti: the first 32 harmonics of Wolfram|Alpha's (Entity["yeti curve", Entity["PopularCurve","YetiCurve"], "Entity[\"PopularCurve\", \"YetiCurve\"]", "popular curve"]), an embedded curve (turning number 1) that is spectacularly non-convex. Rescaled to length  it is a stern test: the flow must unwind the fingers and feet without ever changing the total length.

yetiHarmonics = {
   {{0.03982, -1.02498}, {0.78499, 0.14545}}, {{0.00376, 0.16583}, {-0.02991, 0.00484}},
   {{0.13306, -0.13196}, {0.29221, -0.0513}}, {{-0.01663, 0.41687}, {-0.00118, -0.19014}},
   {{0.00611, 0.04793}, {-0.09758, -0.04151}}, {{-0.12329, -0.06126}, {-0.01566, -0.03933}},
   {{0.04561, 0.04979}, {-0.00148, -0.05911}}, {{0.01826, 0.04786}, {-0.01173, -0.03728}},
   {{0.02115, -0.04354}, {0.05949, 0.00381}}, {{0.02467, -0.02023}, {0.03302, 0.01009}},
   {{0.01833, -0.02521}, {0.03565, -0.01619}}, {{-0.00346, 0.00315}, {-0.00336, 0.00104}},
   {{-0.03492, -0.00742}, {-0.02242, 0.01282}}, {{-0.01925, 0.00019}, {0.01502, 0.0029}},
   {{-0.01367, -0.01352}, {0.01962, -0.01873}}, {{-0.0141, 0.00718}, {-0.01513, -0.01049}},
   {{-0.00526, -0.03014}, {-0.00175, 0.01201}}, {{0.00304, 0.00071}, {-0.01499, 0.00389}},
   {{0.00155, -0.00163}, {-0.00922, -0.00639}}, {{-0.00317, -0.00209}, {-0.00222, 0.00708}},
   {{-0.01201, -0.00106}, {0.00105, 0.00731}}, {{0.00187, -0.00131}, {-0.00541, -0.00511}},
   {{0.00165, -0.00685}, {-0.0015, 0.00439}}, {{-0.00318, 0.00538}, {-0.00582, 0.00333}},
   {{0.00629, 0.01104}, {-0.00369, 0.00258}}, {{0.00776, -0.00202}, {0.01038, 0.00683}},
   {{-0.0001, -0.00124}, {0.00392, 0.00252}}, {{0.00253, -0.00086}, {0.00446, -0.00769}},
   {{0.01163, -0.0063}, {0.01111, 0.00403}}, {{-0.00298, -0.00165}, {0.00042, 0.00375}},
   {{0.00127, -0.00079}, {0.00082, -0.00468}}, {{0.00451, 0.00446}, {-0.00032, 0.0036}}
};
yetiCurve = normalizeCurve[Function[u, Total[Table[
     yetiHarmonics[[k, 1]] Cos[k u] + yetiHarmonics[[k, 2]] Sin[k u],
     {k, Length[yetiHarmonics]}]]]];
drawPolygon[sampleCurve[yetiCurve, 200], ImageSize -> 240]

The fast numeric Fourier basis (the symbolic fourierBasis below is fine for building the solver's fixed collocation matrices, but interactive scrubbing wants this vectorized form):

basisEval[M_, ts_] := Transpose[Prepend[
    Flatten[Table[{Cos[k ts], Sin[k ts]}, {k, M}], 1],
    ConstantArray[1., Length[ts]]]];

We integrate the gauge-fixed flow  as an ODE system for Fourier coefficients , evaluating the right-hand side by collocation on  points. The nonlocal coefficient is computed on the same grid, with one refinement: we use  with a small gain . On the exact flow  and the extra term vanishes identically, so this is the same equation; but for the truncated dynamics it acts as a restoring force, , that pins the length to  instead of letting spectral truncation error accumulate. (Without it, the violent initial transient of a spiky curve like the yeti costs a permanent 0.3% length error — and a curve of the wrong length converges to a circle of the wrong radius, which quietly wrecks every diagnostic that compares  with 1.)

(* d-th derivative rows via the closed form  d^n cos(ku)/du^n = k^n cos(ku + n Pi/2):
   purely numeric, so no global symbol needs to be valueless (a symbolic-D version
   here once broke the moment a Manipulate slider left a value in a global) *)
fourierBasis[M_, uu_List, d_ : 0] := Transpose@Join[
    {If[d == 0, ConstantArray[1., Length[uu]], ConstantArray[0., Length[uu]]]},
    Flatten[Table[{k^d Cos[k uu + d Pi/2.], k^d Sin[k uu + d Pi/2.]}, {k, 1, M}], 1]];

$lengthGain = 5.;   (* the constraint-restoring gain lambda *)

abFourier[curve_, M_ : 32, tEnd_ : 6.] :=
  Module[{n = 4 M + 2, uu, du, E0, E1, E2, Proj, c0, rhs, sol, cc, t},
   uu = N@Table[2 Pi j/n, {j, 0, n - 1}];
   du = 2 Pi/n;
   {E0, E1, E2} = Table[fourierBasis[M, uu, d], {d, 0, 2}];
   Proj = PseudoInverse[E0];
   c0 = Proj . N[curve /@ uu];
   rhs[cf_List] := Module[{c = Partition[cf, 2], g1, g2, sp2, sp, k, L, kbar2},
     g1 = E1 . c; g2 = E2 . c;
     sp2 = Total[g1^2, {2}]; sp = Sqrt[sp2];
     k = (g1[[All, 1]] g2[[All, 2]] - g1[[All, 2]] g2[[All, 1]])/sp^3;
     L = du Total[sp];
     kbar2 = (du Total[k^2 sp] - $lengthGain (L - 2 Pi))/L;
     Flatten[Proj . (g2/sp2)] + kbar2 cf];   (* the dilation is linear in the coefficients *)
   sol = NDSolveValue[{cc'[t] == rhs[cc[t]], cc[0] == Flatten[c0]}, cc, {t, 0, tEnd},
     AccuracyGoal -> 10, PrecisionGoal -> 10];
   <|"Coefficients" -> sol, "Modes" -> M, "EndTime" -> tEnd,
     "Curve" -> Function[{t, m}, fourierBasis[M, N@Table[2 Pi j/m, {j, 0, m - 1}]] .
        Partition[sol[t], 2]],
     "Kappa" -> Function[{t, m}, Module[{g1, g2, sp},
        g1 = fourierBasis[M, N@Table[2 Pi j/m, {j, 0, m - 1}], 1] . Partition[sol[t], 2];
        g2 = fourierBasis[M, N@Table[2 Pi j/m, {j, 0, m - 1}], 2] . Partition[sol[t], 2];
        sp = Sqrt[Total[g1^2, {2}]];
        (g1[[All, 1]] g2[[All, 2]] - g1[[All, 2]] g2[[All, 1]])/sp^3]]|>];

The chord–arc machinery (the mathematics behind it is developed in Sections 6 and 7; the movies of Section 3 already need it, because the flow's natural clock is built from it):

(* all vertex pairs: shorter-arc lengths and chord lengths *)
chordArcPairs[P_] := Module[{n = Length[P], h = edgeLengths[P], s, L, ii, jj, dd, ll, ds},
   s = Prepend[Most@Accumulate[h], 0.];
   L = Total[h];
   {ii, jj} = Transpose@Flatten[Table[{i, j}, {i, n - 1}, {j, i + 1, n}], 1];
   dd = Sqrt[Total[(P[[ii]] - P[[jj]])^2, {2}]];
   ds = Abs[s[[ii]] - s[[jj]]];
   ll = Min /@ Transpose[{ds, L - ds}];
   {ll, dd}];

(* abar from the pairs: smallest b with d >= f(l, -Log[b]) for every pair,
   by bisection on Log[b].  Pairs with d >= 2 Sin[l/2] never bind, since
   f < 2 Sin[x/2] always, so they are harmless to include. *)
abarPairs[P_] := Module[{ll, dd, lo = -30., hi, viol},
   {ll, dd} = chordArcPairs[P];
   viol[lgb_] := Min[dd - 2 E^-lgb ArcTan[E^lgb Sin[ll/2]]];
   If[viol[lo] >= 0, Return[0.]];
   hi = 1.;
   While[viol[hi] < 0 && hi < 60, hi = hi + 1.];
   Do[With[{mid = (lo + hi)/2}, If[viol[mid] >= 0, hi = mid, lo = mid]], {50}];
   E^hi];

(* the diagonal contribution a(p,p) = Sqrt[(k^2-1)/2], using the osculating-circle
   curvature of the polygon *)
osculatingK[P_] := 2 Sin[turningAngles[P]]/
    Sqrt[Total[(RotateLeft[P] - RotateRight[P])^2, {2}]];
abarDiag[P_] := Sqrt[Max[0., (Max[osculatingK[P]^2] - 1.)/2.]];

abar[P_] := Max[abarPairs[P], abarDiag[P]];

Finally, three heavier beasts join the menagerie — the octopus, Cthulhu, and the tight double spiral, imported from the phantom-curve experiment (PhantomFlow.nb; 64-harmonic truncations of Wolfram PopularCurve entities, plus the constructed spiral). Their Fourier coefficients are embedded below as iconized data, so this notebook is a single-file download:

monsterHarm=;
spiralHarm=;

harmCurve[harm_] := normalizeCurve[Function[u, Total[Table[
      harm[[k, 1]] Cos[k u] + harm[[k, 2]] Sin[k u], {k, Length[harm]}]]]];
octopusCurve = harmCurve[monsterHarm["Octopus"]];
cthulhuCurve = harmCurve[monsterHarm["Cthulhu"]];
spiralCurve = harmCurve[spiralHarm];
GraphicsRow[Table[With[{P = sampleCurve[c[[2]], 800]},
   Graphics[{GrayLevel[0.1], AbsoluteThickness[1.1], Line[Append[P, First[P]]]},
    AspectRatio -> Automatic, ImageSize -> 210, PlotLabel -> c[[1]]]],
  {c, {{"octopus", octopusCurve}, {"cthulhu", cthulhuCurve}, {"spiral", spiralCurve}}}], ImageSize -> 660]

3. The flow, watched

If you have never seen a curvature flow move, nothing below will land until you have. So before any theorems: here are five curves — the dumbbell, the yeti, the octopus, Cthulhu, and a tight double spiral — evolving under the length-preserving flow, computed by the spectral solver above. Two viewing notes. First, every frame is drawn comoving (recentred on the curve's centroid); Section 4 shows why the raw frames would sail off the screen exponentially. Second, the time slider is **not uniform in **. A feature of curvature  is erased in time , so the interesting motion of a wiggly curve is crammed into the first instants; uniform  shows a twitch and then hours of nothing. The natural clock is the stage , built from the running chord–arc scale  of Section 7: it starts at  (the constant, determined by the initial curve, that the Andrews–Bryan barrier runs on), ticks like  through the transient — one tick per dying feature scale — and settles into lockstep with  as the curve rounds. The slider is uniform in ; the label reports , , and .

(* ClearAll before every memoized definition: re-evaluating an edited definition
   without it leaves memoized values computed by the OLD definition in place *)
ClearAll[flowRun, stageMap, extentOf];
flowCurves = <|"dumbbell" -> {dumbbellCurve, 32}, "yeti" -> {yetiCurve, 64},
   "octopus" -> {octopusCurve, 80}, "cthulhu" -> {cthulhuCurve, 80},
   "spiral" -> {spiralCurve, 80}|>;
flowRun[nm_] := flowRun[nm] = abFourier[flowCurves[nm][[1]], flowCurves[nm][[2]], 6.];
curveAt[nm_, t_, m_] := Module[{c = Partition[flowRun[nm]["Coefficients"][t], 2]},
   basisEval[flowCurves[nm][[2]], N@Table[2 Pi j/m, {j, 0, m - 1}]] . c];

(* the stage theta(t) = -Log[abar(t)], sampled log-densely, forced monotone,
   truncated when abar falls below measurability *)
stageMap[nm_] := stageMap[nm] = Module[{ts, acc = -Infinity, out = {}},
   ts = Join[{0.}, 10.^Range[-4., Log10[6.], (Log10[6.] + 4.)/27]];
   Do[Module[{ab = abar[resample[curveAt[nm, tt, 900], 300]]},
     If[ab < 10.^-3, Break[]];
     acc = Max[acc, -Log[ab]];
     AppendTo[out, {tt, acc}]], {tt, ts}];
   out];
timeOfStage[nm_, th_] := Module[{sm = stageMap[nm]},
   If[th <= sm[[1, 2]], sm[[1, 1]],
    If[th >= sm[[-1, 2]], sm[[-1, 1]],
     Interpolation[Map[Reverse, sm], InterpolationOrder -> 1][th]]]];
tbarOf[nm_] := -stageMap[nm][[1, 2]];
(* fixed frame for the whole movie: the largest comoving extent the curve ever
   attains, sampled on the stage grid (mid-flow elongated shapes are bigger than
   both the initial curve and the final unit circle) *)
extentOf[nm_] := extentOf[nm] = 1.05 Max[Table[
     Module[{P = curveAt[nm, tt, 500]}, Max[Abs[Map[# - Mean[P] &, P]]]],
     {tt, stageMap[nm][[All, 1]]}]];

(* warm everything now so the movies scrub instantly *)
Grid[Prepend[Table[{nm, flowCurves[nm][[2]],
    NumberForm[tbarOf[nm], {4, 2}],
    NumberForm[stageMap[nm][[-1, 2]], {4, 2}]}, {nm, Keys[flowCurves]}],
  {"curve", "modes M", "tbar", "final stage"}], Frame -> All]

The flow, live. Pick a curve and drag the stage.

Manipulate[
 Module[{sm = stageMap[nm], th, tt, P, ctr, half},
  th = sm[[1, 2]] + stg (sm[[-1, 2]] - sm[[1, 2]]);
  tt = timeOfStage[nm, th];
  P = curveAt[nm, tt, 520]; ctr = Mean[P]; half = extentOf[nm];
  Graphics[{GrayLevel[0.1], AbsoluteThickness[1.9],
    Line[Append[#, First[#]] &@Map[# - ctr &, P]]},
   PlotRange -> {{-half, half}, {-half, half}}, AspectRatio -> 1,
   ImageSize -> 460, Background -> White,
   PlotLabel -> Style[Row[{"θ = ", NumberForm[th, {5, 2}],
       "    t = ", ScientificForm[tt, 3],
       "    t - tbar = ", NumberForm[tt - tbarOf[nm], {5, 2}]}], 12,
     GrayLevel[0.2]]]],
 {{nm, "yeti", "curve"}, Keys[flowCurves]},
 {{stg, 0., "stage"}, 0., 1.},
 ControlPlacement -> Top, SaveDefinitions -> True]

And the same five films as a contact sheet, six stages each — equal steps of , so each column represents the same amount of geometric progress:

GraphicsGrid[Table[Module[{sm = stageMap[nm], th, tt, P, ctr, half},
    th = sm[[1, 2]] + q (sm[[-1, 2]] - sm[[1, 2]]);
    tt = timeOfStage[nm, th];
    P = curveAt[nm, tt, 400]; ctr = Mean[P]; half = extentOf[nm];
    Graphics[{GrayLevel[0.1], AbsoluteThickness[1.2],
      Line[Append[#, First[#]] &@Map[# - ctr &, P]]},
     PlotRange -> {{-half, half}, {-half, half}}, AspectRatio -> 1,
     ImageSize -> 130,
     PlotLabel -> Style[Row[{"θ ", NumberForm[th, {4, 1}]}], 8]]],
   {nm, Keys[flowCurves]}, {q, {0., 0.2, 0.4, 0.6, 0.8, 1.}}],
 ImageSize -> 800]

Watch for the shared choreography: whatever the initial anatomy — lobes, fingers, tentacles, or three turns of spiral — the flow erases features from the smallest scale upward (that is the -clock's whole point), passes through a brief "generic blob" phase, and then rounds out in shape while, invisibly here, the whole picture drifts (Section 4). The theorem this notebook is about says the roundness is not luck: it comes with quantitative, exponential-in- certificates.

4. Trusting the solver: conservation, the runaway centre, and the dictionary

First validation: the unit circle at the origin is stationary.

solCirc = abFourier[circleCurve, 8, 4.];
Max@Table[Abs[Sqrt[Total[#^2]] & /@ solCirc["Curve"][t, 64] - 1.], {t, 0., 4., 0.5}]

Now the lopsided curve. Length is conserved to  after the initial spectral-truncation offset of a few times , and the circle-fit residual decays — the shape becomes round:

solLop = abFourier[lopsidedCurve, 32, 6.];
Grid[Prepend[Table[{t, NumberForm[polyLength[solLop["Curve"][t, 800]], {10, 8}],
     ScientificForm[circleResidual[solLop["Curve"][t, 400]], 3]}, {t, 0., 6.}],
   {"t", "length", "circle-fit residual"}], Frame -> All]

The residual falls by a factor of about  per unit time, the exponential rate the theorem predicts. But look at where the curve went:

circleFit[solLop["Curve"][6., 400]]

The shape is a unit circle to nine digits — centred more than a hundred units from the origin. This is the translation instability of Section 1 in action, and it deserves its own subsection.

The runaway centre

In the dictionary of Section 1,  where  is the shrinking CSF solution and . CSF shrinks the curve to its extinction point , so  and since , the centre runs off like  — unless the origin is the extinction point. The flow converges to the unit circle only modulo translation, and the translation is not a small correction: it is exponentially divergent for every choice of origin but one. (This corrects a parenthetical in the proof of Theorem 6.2 in the lecture notes, which asserts that the pointwise velocity  decays and hence the centre of mass converges: the numerics here show the centre of the lopsided solution at  by  and growing like . The statement "modulo translation" in the theorem is what is actually true.)

ListLogPlot[Table[{t, Norm[First@circleFit[solLop["Curve"][t, 400]]]}, {t, 0., 6., 0.25}],
 Joined -> True, Mesh -> All, AxesLabel -> {"t", "|centre|"}, ImageSize -> 380,
 PlotLabel -> "the fitted-circle centre runs away like E^t",
 GridLines -> {None, None},
 Epilog -> {Dashed, Line[{{0, Log[0.25]}, {6, Log[0.25 E^6]}}]}]

We can even recover the extinction point from the runaway:  where  is the fitted centre and , computable from the solution. Two rounds of this shooting correction (each round divides the residual drift by ) pin the origin to the extinction point, and the same flow then stays within half a unit of the origin through  instead of sailing 124 units away — the residual drift is just  times the remaining error in , so every extra correction round buys another factor of :

kbar2At[sol_, t_] := Module[{k = sol["Kappa"][t, 800], P = sol["Curve"][t, 800]},
   Total[k^2 dualLengths[P]]/polyLength[P]];
logRho[sol_, T_] := Module[{ts = Range[0., T, T/120.], kb},
   kb = kbar2At[sol, #] & /@ ts;
   Total[(T/120.) (Most[kb] + Rest[kb])/2]];

x0 = {0., 0.};
Do[Module[{sol = abFourier[Function[u, lopsidedCurve[u] - x0], 32, 6.]},
   x0 = x0 + First[circleFit[sol["Curve"][6., 400]]]/E^logRho[sol, 6.]], {2}];
x0

solCtr = abFourier[Function[u, lopsidedCurve[u] - x0], 32, 6.];
Grid[Prepend[Table[{t, NumberForm[Norm[First@circleFit[solLop["Curve"][t, 400]]], {6, 3}],
     ScientificForm[Norm[First@circleFit[solCtr["Curve"][t, 400]]], 2]}, {t, 0., 6.}],
   {"t", "|centre|, original origin", "|centre|, origin at x0"}], Frame -> All]

The visual counterpart of these tables is Section 3: every movie there is drawn in this comoving frame, because the raw frame exits stage right at speed .

And the yeti, which needs more modes () to resolve its fingers through the violent early transient — watch the length column: even while the curvature is enormous,  never strays more than a few times  from :

solYeti = abFourier[yetiCurve, 64, 4.];
Grid[Prepend[Table[{t, NumberForm[polyLength[solYeti["Curve"][t, 1600]], {10, 7}],
     ScientificForm[circleResidual[solYeti["Curve"][t, 400]], 3]},
    {t, {0., 0.1, 0.25, 0.5, 1., 2., 4.}}],
   {"t", "length", "circle-fit residual"}], Frame -> All]

The dictionary as an independent check

The dictionary of Section 1 is exact, so we can cross-examine the solver: run classical CSF with a completely different discretization ((TextRefLink[NDSolve, paclet:ref/NDSolve, https://reference.wolfram.com/language/ref/NDSolve.html])'s method of lines, pseudospectral in space), rescale each time slice by , reparametrize time by , and compare pointwise with the direct solution. Both solvers use the same tangential gauge, and the gauge commutes with the dictionary, so the two parametrized curves should agree point by point, not merely as point sets.

csfMOL[curve_, tEnd_, np_ : 129] := Module[{x, y, t, u}, NDSolveValue[{
    D[x[t, u], t] == D[x[t, u], {u, 2}]/(D[x[t, u], u]^2 + D[y[t, u], u]^2),
    D[y[t, u], t] == D[y[t, u], {u, 2}]/(D[x[t, u], u]^2 + D[y[t, u], u]^2),
    x[0, u] == curve[u][[1]], y[0, u] == curve[u][[2]],
    x[t, 0] == x[t, 2 Pi], y[t, 0] == y[t, 2 Pi]},
   {x, y}, {t, 0, tEnd}, {u, 0, 2 Pi},
   Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
      "DifferenceOrder" -> "Pseudospectral", "MinPoints" -> np, "MaxPoints" -> np}}]];

A0 = polyArea[sampleCurve[lopsidedCurve, 800]];
tauMax = 0.93 A0/(2 Pi);       (* CSF dies at A0/2Pi; stop at 93% of the way there *)
{xf, yf} = csfMOL[lopsidedCurve, tauMax, 129];
molCurve[tau_, m_] := With[{uu = N@Table[2 Pi j/m, {j, 0, m - 1}]},
   Transpose[{xf[tau, uu], yf[tau, uu]}]];

(* the dictionary: rho(tau) and the intrinsic clock t(tau) *)
LtildeTab = Table[{tau, polyLength[molCurve[tau, 800]]}, {tau, 0., tauMax, tauMax/400}];
rho2 = Interpolation[{#[[1]], (2 Pi/#[[2]])^2} & /@ LtildeTab];
tOfTau = Interpolation[Module[{ts = LtildeTab[[All, 1]]},
    Transpose[{ts, Prepend[Accumulate@Table[
         NIntegrate[rho2[s], {s, ts[[i]], ts[[i + 1]]}], {i, Length[ts] - 1}], 0.]}]]];

Grid[Prepend[Table[Module[{t = tOfTau[tau], rho = Sqrt[rho2[tau]]},
     {NumberForm[tau, {4, 3}], NumberForm[t, {4, 3}],
      ScientificForm[Max@Abs[solLop["Curve"][t, 400] - rho molCurve[tau, 400]], 2]}],
    {tau, 0.1 tauMax, tauMax, 0.166 tauMax}],
   {"CSF time tau", "intrinsic time t", "max pointwise difference"}], Frame -> All]

Agreement to a few times  or better over the whole comparable range — and note the two clocks: the entire CSF lifetime  unrolls to , and by  the intrinsic clock has only reached . The direct solver cruises to  without noticing that, in CSF terms, it is resolving the last  of the curve's life. That is the practical payoff of the intrinsic clock, and the numerical face of the equivalence between the Andrews–Bryan theorem and Grayson's theorem. As a bonus, the CSF solution exhibits the extinction point directly — the curve near its death huddles around , and it is the same point the shooting correction found above:

{Mean[molCurve[tauMax, 400]], x0}

5. The Andrews–Bryan theorem

For  let  be the chord (extrinsic) distance and  the length of the shorter arc between them. The engine of the proof is the comparison function

For each fixed  it is strictly concave in  with ; it is strictly increasing in ; as  it tends to 0 and as  it tends to , the chord–arc profile of the unit circle. It is a supersolution machine: it satisfies the key identity

identically, which is what makes the maximum-principle argument close (Section 6 checks this symbolically).

The chord–arc scale. For a fixed embedded curve of length  define, for ,

where  is the unique solution of . Equivalently, . The function  extends continuously to the diagonal with  — the cubic Taylor coefficient of  encodes curvature, which is the mechanism converting a chord–arc estimate into a pointwise curvature bound. Set

Theorem (Andrews–Bryan 2011; distance comparison). *If  is a smooth embedded solution of the flow with length  and  is computed from the initial curve, then for all  and all ,*

Since  off the diagonal, the flow preserves embeddedness with a quantitative modulus, and since , the chord–arc profile of the evolving curve converges to that of the unit circle.

Theorem (exponential curvature bound). *Letting  in the distance comparison,*

*Moreover the normalized scale  is non-increasing along the flow.*

Theorem (convergence). *The flow exists for all time, and  exponentially in : e.g. , and the curves converge exponentially in every  norm, modulo translation, to a unit circle.*

Corollary (Grayson's theorem). *Under CSF, any smooth embedded closed curve stays embedded, shrinks to a point  at time , and the rescaled curves  converge in  to a round unit circle.* (Run the dictionary of Section 1 backwards; the classical statement is exactly the length-preserving theorem in the intrinsic clock.)

What is remarkable about the argument — and what makes it such good numerical material — is what it does not use: no monotonicity formula, no blow-up analysis, no classification of self-similar solutions. Everything is controlled by one two-point function sliding along one explicit barrier, and every object in it can be evaluated on a polygon.

6. The comparison function, symbolically

The two "verification lemmas" that the Andrews–Bryan paper leaves to the reader are a perfect job for the computer. First the function itself and its -derivative,  with :

fAB[x_, t_] := 2 E^t ArcTan[E^-t Sin[x/2]];

Block[{x, t},
 Simplify[D[fAB[x, t], x] - Cos[x/2]/(1 + E^(-2 t) Sin[x/2]^2)]]

The key identity (Lemma 4.5 of the notes): , *identically in  and * — this is the entire reason  works, and the notes call its discovery "serendipitous":

Block[{x, t},
 FullSimplify[
  4 D[fAB[x, t], {x, 2}] + fAB[x, t] -
    (4 D[fAB[x, t], x]/Sin[x/2]) (D[fAB[x, t], x] - Cos[x/2]) -
    D[fAB[x, t], t],
  Assumptions -> 0 < x < 2 Pi && Element[t, Reals]]]

Monotonicity in :  with , and , , so  increases in  — this is what makes  uniquely solvable and the barrier "slide":

Block[{x, t, w},
 With[{z = E^-t Sin[x/2]},
  {Simplify[D[fAB[x, t], t] - 2 E^t (ArcTan[z] - z/(1 + z^2))],
   Simplify[D[ArcTan[w] - w/(1 + w^2), w] - 2 w^2/(1 + w^2)^2]}]]

The Taylor expansion that converts chords into curvature (Lemma 3.1 of the notes): , so matching against the expansion  of a genuine curve near the diagonal reads off , i.e.  where positive:

(* Clear, not Block: this cell's OUTPUT contains x and a, and a Block result
   re-evaluates after the saved values are restored *)
Clear[x, a];
Assuming[a > 0, Series[fAB[x, -Log[a]], {x, 0, 4}] // Simplify]

And the convexity of  (Lemma 4.4 of the notes, used to linearize the curvature integral): substituting , the second derivative reduces to , which is nonnegative because :

(* Clear, not Block: symbolic output containing p *)
Clear[z, p];
Simplify[(D[ArcCos[z]^2, {z, 2}] /. z -> Cos[p]) Sin[p]^3/2,
 Assumptions -> 0 < p < Pi/2]

Finally the picture worth having in mind:  interpolates monotonically between 0 (as ) and the unit-circle profile  (as ). Huisken's 1998 comparison theorem is a static lower bound of the form ; the Andrews–Bryan refinement replaces it with this sliding barrier, which tightens to the circle profile as the flow runs:

Plot[Evaluate[Append[Table[fAB[x, t], {t, -2, 3}], 2 Sin[x/2]]], {x, 0, 2 Pi},
 PlotStyle -> Append[Table[ColorData["DeepSeaColors"][(t + 2)/5.5], {t, -2, 3}],
    Directive[Black, Dashed]],
 AxesLabel -> {"x", "f(x,t)"}, ImageSize -> 420,
 PlotLabel -> "the sliding barrier, t = -2, ..., 3; dashed: 2 Sin[x/2]"]

7. The chord–arc scale, computationally

Everything in the theorem is driven by one number computed from the initial curve: , the sup of . On a polygon this is genuinely computable. The code lives in the Setup section (the movies already used it for the stage clock); the pairs sweep uses the equivalence :  is the smallest  making every pair satisfy the inequality, which is one bisection over  with a fully vectorized evaluation over all  pairs inside.

Sanity checks. The unit circle has  for every pair and , so . For the ellipse, the diagonal formula  and the pairs sweep should agree — on an ellipse the binding pairs are the near-diagonal ones at the tips, so the two computations are estimating the same limit:

Grid[Prepend[
  Table[Module[{P = sampleCurve[pair[[2]], 400]},
    {pair[[1]], NumberForm[abarPairs[P], {6, 3}], NumberForm[abarDiag[P], {6, 3}],
     NumberForm[Max@Abs[osculatingK[P]], {6, 3}]}],
   {pair, {{"circle", circleCurve}, {"ellipse", ellipseCurve},
     {"lopsided", lopsidedCurve}, {"yeti", yetiCurve}}}],
  {"curve", "abar (pairs)", "abar (diagonal)", "max |k|"}], Frame -> All]

Note the yeti: its  comes from the pairs, not the diagonal — a pair of points on opposite banks of a finger, close in space but far along the curve, beats the maximum curvature by a factor of four. The chord–arc scale is a genuine embeddedness scale, not just a curvature bound; comparison quantities that only see curvature cannot play its role in the proof. (For the gentler curves the two contributions nearly coincide, because the binding pairs are the near-diagonal ones at the point of highest curvature.)

8. The distance comparison: the chord envelope and the sliding barrier

Theorem 4.1 (Section 5) is a statement about every pair of points at once, which makes it hard to see: a scatter of all  chord–arc pairs is a cloud, and clouds move badly. Here is a rendering that moves well. Fix an arc length  and slide the basepoint around the curve: the chord  varies continuously, sweeping out exactly the interval  between two continuous envelope functions. The whole chord–arc data of the curve at one instant is the band between  and  — and the distance-comparison theorem says precisely that this band rides above the sliding barrier , for every  and every , while both band and barrier converge to the circle profile . The squeeze is the theorem.

(* chord envelope of an equal-arclength polygon: for each arc offset, the min and
   max chord over all basepoints *)
chordEnvelope[P_] := Module[{n = Length[P], ells, ms, Ms},
   ells = Table[2. Pi k/n, {k, 1, Floor[n/2]}];
   {ms, Ms} = Transpose@Table[
      Module[{d = Sqrt[Total[(RotateLeft[P, k] - P)^2, {2}]]}, {Min[d], Max[d]}],
      {k, 1, Floor[n/2]}];
   {ells, ms, Ms}];

envelopePanel[nm_, tt_] := Module[{tb = tbarOf[nm], e},
   e = chordEnvelope[resample[curveAt[nm, tt, 900], 360]];
   Graphics[{
     {RGBColor[0.30, 0.50, 0.80, 0.30], EdgeForm[],
      Polygon[Join[Transpose[{e[[1]], e[[3]]}], Reverse@Transpose[{e[[1]], e[[2]]}]]]},
     {RGBColor[0.16, 0.35, 0.75], AbsoluteThickness[1.4],
      Line[Transpose[{e[[1]], e[[3]]}]], Line[Transpose[{e[[1]], e[[2]]}]]},
     First@Plot[fAB[x, tt - tb], {x, 0.001, Pi},
       PlotStyle -> Directive[RGBColor[0.95, 0.55, 0.10], AbsoluteThickness[2.4]]],
     First@Plot[2 Sin[x/2], {x, 0, Pi},
       PlotStyle -> Directive[GrayLevel[0.1], Dashed, AbsoluteThickness[1.2]]]},
    Axes -> False, Frame -> True,
    FrameLabel -> {"arc ℓ", "chord"},
    PlotRange -> {{0, Pi}, {0, 2.12}}, AspectRatio -> 0.9, ImageSize -> 330]];

The flow and its envelope, side by side. Left: the curve, comoving. Right: the chord band (blue, between  and ), the sliding barrier  (orange), and the circle profile (dashed). Drag the stage and watch the band get squeezed onto the dashed curve with the barrier rising underneath — never touching the band, exactly as Theorem 4.1 requires:

Manipulate[
 Module[{sm = stageMap[nm], th, tt, P, ctr, half},
  th = sm[[1, 2]] + stg (sm[[-1, 2]] - sm[[1, 2]]);
  tt = timeOfStage[nm, th];
  P = curveAt[nm, tt, 520]; ctr = Mean[P]; half = extentOf[nm];
  GraphicsRow[{
    Graphics[{GrayLevel[0.1], AbsoluteThickness[1.8],
      Line[Append[#, First[#]] &@Map[# - ctr &, P]]},
     PlotRange -> {{-half, half}, {-half, half}}, AspectRatio -> 1,
     ImageSize -> 330],
    envelopePanel[nm, tt]},
   ImageSize -> 700,
   PlotLabel -> Style[Row[{"θ = ", NumberForm[th, {5, 2}],
       "    t = ", ScientificForm[tt, 3],
       "    t - tbar = ", NumberForm[tt - tbarOf[nm], {5, 2}]}], 12,
     GrayLevel[0.2]]]],
 {{nm, "octopus", "curve"}, Keys[flowCurves]},
 {{stg, 0., "stage"}, 0., 1.},
 ControlPlacement -> Top, SaveDefinitions -> True]

A contact sheet of the squeeze for the spiral — the curve above, its envelope below, at four stages:

Module[{nm = "spiral", sm, stages},
 sm = stageMap[nm];
 stages = sm[[1, 2]] + {0.05, 0.35, 0.6, 0.95} (sm[[-1, 2]] - sm[[1, 2]]);
 GraphicsGrid[Transpose@Table[Module[{tt = timeOfStage[nm, th], P, ctr, half},
     P = curveAt[nm, tt, 420]; ctr = Mean[P]; half = extentOf[nm];
     {Graphics[{GrayLevel[0.1], AbsoluteThickness[1.2],
        Line[Append[#, First[#]] &@Map[# - ctr &, P]]},
       PlotRange -> {{-half, half}, {-half, half}}, AspectRatio -> 1,
       ImageSize -> 180,
       PlotLabel -> Style[Row[{"θ = ", NumberForm[th, {4, 1}]}], 9]],
      Show[envelopePanel[nm, tt], ImageSize -> 190]}],
    {th, stages}], ImageSize -> 800]]

And the theorem as one number per curve: the worst slack  over the whole flow. It should never be negative; the values below sit at zero to within the  chord error of the 360-gon (a few times ), which is the resolution of this check:

Grid[Prepend[Table[Module[{worst = Infinity, tb = tbarOf[nm]},
    Do[Module[{e = chordEnvelope[resample[curveAt[nm, tt, 900], 360]]},
      worst = Min[worst, Min[e[[2]] - fAB[e[[1]], tt - tb]]]],
     {tt, Join[{0.}, 10.^Range[-4., Log10[6.], 0.25]]}];
    {nm, ScientificForm[worst, 3]}], {nm, Keys[flowCurves]}],
  {"curve", "worst band slack over the flow"}], Frame -> All]

9. The curvature bound and the monotone scale

Everything in the Andrews–Bryan theorem is now on the table: compute  from the initial curve, then watch the inequalities hold along the flow. For the smooth flow we upgrade the diagonal part of  to use the spectral curvature (the polygon estimator's  bias would otherwise swamp  once it decays below ):

abarSmooth[sol_, t_, m_ : 600] := Module[{P = sol["Curve"][t, m], kmax},
   kmax = Max@Abs[sol["Kappa"][t, 1600]];
   Max[abarPairs[P], Sqrt[Max[0., (kmax^2 - 1.)/2.]]]];

tbarLop = Log[abarSmooth[solLop, 0.]]

The exponential curvature bound, . The lopsided curve's  is attained on the diagonal, so at  the bound is an equality by construction. The actual curvature then plunges far below the envelope — the flow destroys the initial curvature spike much faster than the worst case — before settling to decay at the theorem's asymptotic rate , parallel to the bound. The two panels pair each quantity with its own bound and shade the slack between them; the left is the sup bound, the right the  bound of Proposition 6.1:

Module[{ts = Range[0., 6., 0.25], f1, b1, f2, b2},
 f1 = Table[{t, Max[Abs[solLop["Kappa"][t, 1600]]]^2 - 1}, {t, ts}];
 b1 = Table[{t, 2 E^(-2 (t - tbarLop))}, {t, ts}];
 f2 = Table[{t, Module[{P = solLop["Curve"][t, 1600], k = solLop["Kappa"][t, 1600]},
     Total[(k - 1)^2 dualLengths[P]]]}, {t, ts}];
 b2 = Table[{t, 4 Pi E^(-2 (t - tbarLop))}, {t, ts}];
 GraphicsRow[{
   ListLogPlot[{b1, f1}, Joined -> True, Filling -> {1 -> {2}},
    FillingStyle -> RGBColor[0.95, 0.68, 0.10, 0.18],
    PlotStyle -> {Directive[RGBColor[0.95, 0.55, 0.10], AbsoluteThickness[2.]],
      Directive[RGBColor[0.16, 0.35, 0.75], AbsoluteThickness[1.6]]},
    AxesLabel -> {"t", None}, ImageSize -> 360, PlotRange -> All,
    PlotLabel -> Style[
      "max k² - 1 (blue) under its bound 2E^(-2(t-tbar)) (orange)", 10]],
   ListLogPlot[{b2, f2}, Joined -> True, Filling -> {1 -> {2}},
    FillingStyle -> RGBColor[0.95, 0.68, 0.10, 0.18],
    PlotStyle -> {Directive[RGBColor[0.95, 0.55, 0.10], AbsoluteThickness[2.]],
      Directive[RGBColor[0.16, 0.35, 0.75], AbsoluteThickness[1.6]]},
    AxesLabel -> {"t", None}, ImageSize -> 360, PlotRange -> All,
    PlotLabel -> Style[
      "∫(k-1)² ds (blue) under its bound 4πE^(-2(t-tbar)) (orange)", 10]]},
  ImageSize -> 760]]

The distance comparison itself has its own moving picture — the chord envelope of Section 8 — so here we only record the two pointwise consequences.

The monotone quantity. The theorem can be restarted at any time, which says the normalized chord–arc scale  is non-increasing (Remark 5.2 of the notes). It collapses steeply while the flow erases the initial concavities, then settles onto a perfect  decay — the flow carries a Lyapunov-like quantity of pointwise-extremal type, found nowhere near the usual integral monotonicity formulas:

ListLinePlot[Table[{t, E^t abarSmooth[solLop, t]}, {t, 0., 6., 0.25}],
 Mesh -> All, ScalingFunctions -> "Log", AxesLabel -> {"t", "E^t abar(t)"},
 ImageSize -> 420, PlotLabel -> "Remark 5.2: E^t abar(t) is non-increasing"]

10. Polygonal translations of the flow

So much for the smooth world; the last part of the notebook is about the discrete one. Now the discrete side. A polygonal translation of  needs two ingredients, and each has several "obvious" choices.

**Ingredient 1: the curvature vector .* We take the four discrete curvatures from Crane's Lecture 1 (see the curve shortening flow notebook for their derivations and their inequivalence), in the pointwise* scaling that is consistent with the smooth : turning angle , length variation , Steiner  — each along the bisector normal — and the osculating circle, for which both the curvature  and the normal come from the circumcircle of : .

Ingredient 2: the dilation coefficient. Three translations of , in increasing order of cunning:

Quadrature:  — just discretize the integral. (We use  rather than  throughout, for the same self-correcting reason as in Section 2's solver: it makes length preservation exact-in-principle whatever the current length is.)

ExactLength: choose  so that the discrete length is exactly conserved, . The gradient of polygon length is , which is exactly the (integrated) length-variation curvature normal  — Crane's definition B repays a debt here. So

and  holds identically, for any polygon and any curvature model, before any time discretization.

RescaledCSF: translate the dictionary instead of the equation — take a plain discrete CSF step , then rescale the whole polygon about the origin to length exactly . This is Grayson's theorem's view of the flow, discretized.

discreteK[P_, "TurningAngle"] := turningAngles[P]/dualLengths[P];
discreteK[P_, "LengthVariation"] := 2 Sin[turningAngles[P]/2]/dualLengths[P];
discreteK[P_, "Steiner"] := 2 Tan[turningAngles[P]/2]/dualLengths[P];
discreteK[P_, "OsculatingCircle"] := osculatingK[P];

circumcenters[P_] := Module[{a = RotateRight[P], b = P, c = RotateLeft[P], na, nb, nc, det},
   {na, nb, nc} = Total[#^2, {2}] & /@ {a, b, c};
   det = 2 (a[[All, 1]] (b[[All, 2]] - c[[All, 2]])
           + b[[All, 1]] (c[[All, 2]] - a[[All, 2]])
           + c[[All, 1]] (a[[All, 2]] - b[[All, 2]]));
   Transpose[{(na (b[[All, 2]] - c[[All, 2]]) + nb (c[[All, 2]] - a[[All, 2]])
               + nc (a[[All, 2]] - b[[All, 2]]))/det,
              (na (c[[All, 1]] - b[[All, 1]]) + nb (a[[All, 1]] - c[[All, 1]])
               + nc (b[[All, 1]] - a[[All, 1]]))/det}]];

(* kappa_i N_i for each model; OsculatingCircle uses the circumcenter normal *)
curvatureVector[P_, model_] := If[model === "OsculatingCircle",
   With[{d = circumcenters[P] - P}, d/Total[d^2, {2}]],
   discreteK[P, model] vertexNormals[P]];

dilation[P_, kv_, "Quadrature"] := With[{l = dualLengths[P]},
   Total[Total[kv^2, {2}] l]/Total[l]];
dilation[P_, kv_, "ExactLength"] := With[{gradL = tangents[P] - RotateRight[tangents[P]]},
   -Total[gradL kv, 2]/Total[gradL P, 2]];

abVelocity[P_, model_, dil_] := With[{kv = curvatureVector[P, model]},
   If[dil === "RescaledCSF", kv, kv + dilation[P, kv, dil] P]];

Three small propositions, each provable in a few lines and each checkable to machine precision.

Proposition 1 (ExactLength does what it says). By construction  identically. The other two dilations only preserve length to  per step (Quadrature) or restore it by fiat (RescaledCSF).

Proposition 2 (regular polygons are stationary). For a regular -gon centred at the origin, every model's velocity is radial and identical at all vertices by symmetry; the ExactLength constraint then forces it to vanish — a nonzero radial velocity changes the length. So regular polygons of length , at any , are exact fixed points, for all four curvature models. (They are the discrete unit circles of this flow.)

Proposition 3 (the osculating-circle model inherits ROUND — and the runaway). Let the vertices lie, arbitrarily spaced, on a circle of radius  centred at . Every circumcircle is that circle, so , and Quadrature gives the dilation ; the velocity of every vertex is then : concyclic polygons move by exact rigid translation with velocity , staying concyclic forever, and are exactly stationary iff centred at the origin. The discrete flow reproduces the smooth flow's translation instability  literally, uneven spacing and all.

(* Prop 2: all four models, ExactLength, regular 24-gon of length 2 Pi *)
reg = With[{P = Table[{Cos[u], Sin[u]}, {u, 0., 2 Pi (1 - 1/24), 2 Pi/24}]},
   (2 Pi/polyLength[P]) P];
Grid[Prepend[Table[{m, ScientificForm[Max@Abs[abVelocity[reg, m, "ExactLength"]], 2]},
    {m, {"TurningAngle", "LengthVariation", "Steiner", "OsculatingCircle"}}],
  {"model", "max |velocity| on a regular 24-gon"}], Frame -> All]

(* Prop 3: unevenly spaced concyclic vertices; origin-centred vs shifted *)
SeedRandom[3];
angs = Sort@Table[2 Pi k/24 + 0.5 (2 Pi/24) RandomReal[{-1, 1}], {k, 0, 23}];
Circ0 = Transpose[{Cos[angs], Sin[angs]}];
{Max@Abs[abVelocity[Circ0, "OsculatingCircle", "Quadrature"]],
 abVelocity[Map[# + {0.7, -0.3} &, Circ0], "OsculatingCircle", "Quadrature"][[1 ;; 3]]}

The first number is machine zero; the shifted circle's velocity is exactly  at every vertex.

The flow loop. Explicit Euler with a step set from the shortest edge (, as for any parabolic vertex flow), the same (TextRefLink[TrueQ, paclet:ref/TrueQ, https://reference.wolfram.com/language/ref/TrueQ.html]) guards against folded vertices as in the curve shortening flow notebook, optional arclength resampling (vertices still cluster in concave regions — the dilation term does not redistribute them), and optional exact length renormalization after each step, which plays the role of the -feedback of Section 2's solver.

abFlow[P0_, model_, dil_, opts___] := Module[{
    P = N[P0], t = 0., steps, safety, tEnd, minEdge, red, stride, renorm,
    n = Length[P0], v, h, sp, dt, i = 0, o},
   o = Association[{"Steps" -> 4000000, "Safety" -> 0.25, "EndTime" -> 1.,
      "MinEdge" -> 10.^-7, "Redistribute" -> False, "FrameStride" -> 20,
      "RenormalizeLength" -> True}, opts];
   steps = o["Steps"]; safety = o["Safety"]; tEnd = o["EndTime"];
   minEdge = o["MinEdge"]; red = o["Redistribute"]; stride = o["FrameStride"];
   renorm = o["RenormalizeLength"];
   Reap[Sow[{0., P}];
     Do[
      h = Min[edgeLengths[P]];
      If[! TrueQ[h > minEdge] || ! TrueQ[t < tEnd], Break[]];
      v = abVelocity[P, model, dil];
      sp = Sqrt[Max[Total[v^2, {2}]]];
      If[! TrueQ[0 < sp < Infinity], Break[]];
      dt = Min[safety Min[h^2, h/sp], tEnd - t];
      P = P + dt v; t += dt; i++;
      If[dil === "RescaledCSF" || renorm, P = (2 Pi/polyLength[P]) P];
      If[red, P = resample[P, n]];
      If[Mod[i, stride] == 0, Sow[{t, P}]], steps];
     If[Mod[i, stride] != 0, Sow[{t, P}]]][[2, 1]]];

How much does the choice of dilation matter? Switch the per-step renormalization off and measure the length drift over  — and the honest answer is: for a resolved polygon, less than the time discretization error. ExactLength conserves length exactly for the vertex ODE, but explicit Euler steps along its tangent, so both it and Quadrature drift at ; RescaledCSF is exact by fiat. The real distinctions are structural (the propositions above), not quantitative:

P80 = With[{P = sampleCurve[lopsidedCurve, 80]}, (2 Pi/polyLength[P]) P];
Grid[Prepend[Table[Module[{run = abFlow[P80, "TurningAngle", dil, "EndTime" -> 1.,
       "RenormalizeLength" -> False]},
     {dil, ScientificForm[polyLength[run[[-1, 2]]] - 2 Pi, 2]}],
    {dil, {"Quadrature", "ExactLength", "RescaledCSF"}}],
   {"dilation", "length drift at t = 1, no renormalization"}], Frame -> All]

Finally, the origin-dependence of Section 1, on the polygon: run the same initial polygon from two origins two units apart. The curvatures agree to solver precision — the two runs are the same curve evolution — while the images separate by far more than the original two units:

runO = abFlow[P80, "TurningAngle", "ExactLength", "EndTime" -> 1.];
runS = abFlow[Map[# + {2., 0.} &, P80], "TurningAngle", "ExactLength", "EndTime" -> 1.];
{Max@Abs[Sort[osculatingK[runO[[-1, 2]]]] - Sort[osculatingK[runS[[-1, 2]]]]],
 Norm[Mean[runS[[-1, 2]]] - Mean[runO[[-1, 2]]]]}

11. Discrete vs smooth

Because both sides preserve length, the discrete and smooth flows share an honest clock — we can compare states at equal , with no area-matching contortions. Sample the initial curve at  vertices, flow to , and measure the largest distance from any vertex to the spectral solution. All combinations are second-order accurate. (The table takes a few minutes, nearly all of it in the  rows: with no resampling, vertices crowd into the concavities, the shortest edge shrinks, and the explicit step  shrinks with its square. Resampling would fix the step count but degrade the observed order to one, exactly as in the curve shortening notebook — so we pay.)

ref1 = solLop["Curve"][0.5, 3000];
dist1 = polylineDistance[ref1];
convErr[n_, model_, dil_] := Module[{P = sampleCurve[lopsidedCurve, n]},
   Max[dist1 /@ abFlow[(2 Pi/polyLength[P]) P, model, dil,
       "EndTime" -> 0.5][[-1, 2]]]];

convTab = Table[{n, convErr[n, "TurningAngle", "ExactLength"],
    convErr[n, "OsculatingCircle", "Quadrature"],
    convErr[n, "LengthVariation", "RescaledCSF"]}, {n, {20, 40, 80, 160}}];
Grid[Prepend[Map[If[NumberQ[#] && # < 1, ScientificForm[#, 3], #] &, convTab, {2}],
  {"n", "TurningAngle+ExactLength", "OscCircle+Quadrature", "LengthVar+RescaledCSF"}],
 Frame -> All]

(* observed order of accuracy between successive n *)
Grid[Prepend[Table[Prepend[Table[NumberForm[
      Log2[convTab[[k, j]]/convTab[[k + 1, j]]], {4, 2}], {j, 2, 4}],
    Row[{convTab[[k, 1]], " -> ", convTab[[k + 1, 1]]}]],
   {k, Length[convTab] - 1}],
  {"n", "TA+EL", "OC+Q", "LV+RCSF"}], Frame -> All]

And the theorem holds on the polygon's own terms: compute  from the initial polygon's chord–arc data, run the discrete flow for a long time, and check the discrete curvature against the envelope. (The polygon's  is a faithful stand-in for the smooth one — the pairs sweep sees exactly the same chords the theorem does.)

tbarPoly = Log[abar[P80]];
runLong = abFlow[P80, "TurningAngle", "ExactLength", "EndTime" -> 6.,
   "Redistribute" -> True, "FrameStride" -> 50];
Grid[Prepend[Table[With[{fr = runLong[[j]]},
     {NumberForm[fr[[1]], {4, 2}],
      ScientificForm[Max[osculatingK[fr[[2]]]]^2 - 1, 3],
      ScientificForm[2 E^(-2 (fr[[1]] - tbarPoly)), 3],
      ScientificForm[circleResidual[fr[[2]]], 2]}],
    {j, 1, Length[runLong], Ceiling[Length[runLong]/10]}],
   {"t", "max k² - 1 (discrete)", "bound 2 E^(-2(t-tbar))", "circle residual"}],
 Frame -> All]

One more lesson from the runaway centre: compare the long discrete run's final circle fit with the smooth one. The shapes agree to many digits, but the positions differ by order one — the exponential translation instability amplifies discretization error in the centre of mass like , even while the shape error stays tiny. Position under this flow is numerically meaningless at large ; shape is rock solid:

{circleFit[runLong[[-1, 2]]], circleFit[solLop["Curve"][runLong[[-1, 1]], 400]]}

12. Demonstrations

Shared scaffolding, rebuilt around three lessons from the earlier sections. First, the interesting motion happens early, so the demonstration runs store their frames log-spaced in time and the sliders scrub uniformly through those frames — the left end of every slider is the first instants of the flow, magnified. Second, the menagerie now includes the octopus and Cthulhu, plus two genuinely low-poly carpenter's-rule configurations after Cantarella–Demaine–Iben–O'Brien: teeth (a 17-gon comb whose spikes nearly touch the bottom rail) and a pentagonal spiral (a 24-gon thickened spiral of 72-degree turns). Third, low-poly polygons expose a real phenomenon: with resampling off, the flow shrinks the edges at sharp corners to zero in finite time, so the runner performs edge collapse — when an edge drops below a threshold, its vertex is merged away and the flow continues with one vertex fewer. Watch the vertex count in the demonstrations: the discrete flow can eat most of a low-poly curve on its way to convexity. (As everywhere, every memoized definition carries its (TextRefLink[ClearAll, paclet:ref/ClearAll, https://reference.wolfram.com/language/ref/ClearAll.html]), memoized functions have fixed arity, and tunables are $-globals.)

ClearAll[smoothRun, discRun, demoTbar];
demoCurves = <|"circle" -> circleCurve, "ellipse" -> ellipseCurve, "blob" -> blobCurve,
   "lopsided" -> lopsidedCurve, "dumbbell" -> dumbbellCurve, "yeti" -> yetiCurve,
   "octopus" -> octopusCurve, "cthulhu" -> cthulhuCurve|>;

(* carpenter's-rule polygons, normalized to length 2 Pi *)
teethPoly[k_ : 7, W_ : 4., H_ : 1., delta_ : 0.08] := Module[{xs, top},
   xs = Range[0., W, W/(2 k)];
   top = Table[{W - xs[[j]], If[OddQ[j], H, delta]}, {j, Length[xs]}];
   Join[{{0., 0.}, {W, 0.}}, top]];
pentSpiralPoly[turns_ : 2, w_ : 0.28, growth_ : 1.13, L0_ : 0.8] := Module[
   {nE = 5 turns + 1, pts = {{0., 0.}}, dir = {1., 0.}, len = L0, rot, core, nv,
    left, right},
   rot = RotationMatrix[72. Degree];
   Do[AppendTo[pts, Last[pts] + len dir]; dir = rot . dir; len *= growth, {nE}];
   core = pts;
   nv[i_] := Module[{tv},
     tv = If[i == 1, core[[2]] - core[[1]],
       If[i == Length[core], core[[-1]] - core[[-2]],
        Normalize[core[[i + 1]] - core[[i]]] + Normalize[core[[i]] - core[[i - 1]]]]];
     tv = Normalize[tv]; {-tv[[2]], tv[[1]]}];
   left = Table[core[[i]] + (w/2) nv[i], {i, Length[core]}];
   right = Table[core[[i]] - (w/2) nv[i], {i, Length[core]}];
   Reverse[Join[left, Reverse[right]]]];
demoPolys = <|"teeth" -> (2 Pi/polyLength[#]) # &@teethPoly[7],
   "pentagon spiral" -> (2 Pi/polyLength[#]) # &@pentSpiralPoly[2]|>;
demoNames = Join[Keys[demoCurves], Keys[demoPolys]];
GraphicsRow[Table[With[{P = demoPolys[nm]},
   Graphics[{GrayLevel[0.1], AbsoluteThickness[1.4], Line[Append[P, First[P]]],
     Red, PointSize[0.014], Point[P]}, AspectRatio -> Automatic, ImageSize -> 300,
    PlotLabel -> Row[{nm, "  (n = ", Length[P], ")"}]]],
  {nm, Keys[demoPolys]}], ImageSize -> 640]

The demonstration runner: the polygonal flow of Section 10 plus edge collapse and log-spaced frame storage. Frames are recorded when  crosses a geometric grid ( per frame from ), so a run of hundreds of thousands of Euler steps stores a few hundred frames covering five decades of time.

demoFlow[P0_, model_, dil_, opts___] := Module[{
    P = N[P0], t = 0., steps, safety, tEnd, red, renorm, collapse,
    v, h, hi, sp, dt, i = 0, nextT, rho = 1.02, o},
   o = Association[{"Steps" -> 400000, "Safety" -> 0.25, "EndTime" -> 4.,
      "CollapseBelow" -> 0., "Redistribute" -> False,
      "RenormalizeLength" -> True, "FirstFrame" -> 10.^-4}, opts];
   steps = o["Steps"]; safety = o["Safety"]; tEnd = o["EndTime"];
   collapse = o["CollapseBelow"]; red = o["Redistribute"];
   renorm = o["RenormalizeLength"]; nextT = o["FirstFrame"];
   Reap[Sow[{0., P}];
     Do[
      Module[{els = edgeLengths[P]}, h = Min[els]; hi = First@Ordering[els, 1]];
      While[collapse > 0 && h < collapse && Length[P] > 5,
       P = Delete[P, Mod[hi, Length[P]] + 1];
       Module[{els = edgeLengths[P]}, h = Min[els]; hi = First@Ordering[els, 1]]];
      If[! TrueQ[h > 10.^-7] || ! TrueQ[t < tEnd], Break[]];
      v = abVelocity[P, model, dil];
      sp = Sqrt[Max[Total[v^2, {2}]]];
      If[! TrueQ[0 < sp < Infinity], Break[]];
      dt = Min[safety Min[h^2, h/sp], tEnd - t];
      P = P + dt v; t += dt; i++;
      If[renorm, P = (2 Pi/polyLength[P]) P];
      If[red, P = resample[P, Length[P]]];
      If[t >= nextT, Sow[{t, P}]; While[nextT <= t, nextT *= rho]], steps];
     Sow[{t, P}]][[2, 1]]];

demoInit[name_, nv_] := If[KeyExistsQ[demoPolys, name], demoPolys[name],
   With[{P = sampleCurve[demoCurves[name], nv]}, (2 Pi/polyLength[P]) P]];
discRun[name_, nv_, model_, dil_, red_] := discRun[name, nv, model, dil, red] =
   demoFlow[demoInit[name, nv], model, dil, "Redistribute" -> red,
    "CollapseBelow" -> If[red, 0., 0.02]];
smoothRun[name_] := smoothRun[name] = abFourier[demoCurves[name],
    Which[name === "yeti", 64, MemberQ[{"octopus", "cthulhu"}, name], 80, True, 32], 4.];
demoTbar[name_, nv_] := demoTbar[name, nv] = Log[abar[demoInit[name, nv]]];
frameOf[run_, u_] := run[[1 + Round[u (Length[run] - 1)]]];
(* one fixed frame per run: the largest comoving extent over the whole trajectory *)
runExtent[run_] := 1.05 Max[Table[Module[{c = Mean[fr[[2]]]},
     Max[Abs[Map[# - c &, fr[[2]]]]]], {fr, run[[1 ;; -1 ;; Max[1, Floor[Length[run]/60]]]]}]];
modelColors = <|"TurningAngle" -> RGBColor[0.85, 0.37, 0.01],
   "LengthVariation" -> RGBColor[0.00, 0.62, 0.45],
   "Steiner" -> RGBColor[0.80, 0.47, 0.65],
   "OsculatingCircle" -> RGBColor[0.60, 0.06, 0.35]|>;
smoothColor = RGBColor[0.10, 0.30, 0.70];
models = {"TurningAngle", "LengthVariation", "Steiner", "OsculatingCircle"};
dilations = {"Quadrature", "ExactLength", "RescaledCSF"};

The flow explorer. Pick a curve (or one of the carpenter's-rule polygons — the vertex slider is ignored for those), a curvature model, a dilation, and scrub the log-time slider. With resampling off, watch two things at once: whether the polygon ever stops being simple (it never has, in our experiments — that is the conjecture at work) and how the vertex count falls as corners collapse. First use of a parameter combination computes the whole trajectory: a few seconds usually, up to a minute for the teeth with resampling off.

Manipulate[
 Module[{run, fr, t, P, c, half, simple},
  run = discRun[nm, n, model, dil, red];
  fr = frameOf[run, uu]; t = fr[[1]]; P = fr[[2]];
  c = Mean[P]; half = runExtent[run];
  simple = SimplePolygonQ[Polygon[P]];
  Column[{
    Graphics[{modelColors[model], AbsoluteThickness[1.8],
      Line[Append[#, First[#]] &@Map[# - c &, P]],
      GrayLevel[0.15], PointSize[0.012], Point[Map[# - c &, P]]},
     PlotRange -> {{-half, half}, {-half, half}}, AspectRatio -> 1,
     ImageSize -> 430, Frame -> True, FrameTicks -> None],
    Grid[{{"t", "vertices", "L", "circle residual", "embedded?"},
      {ScientificForm[t, 3], Length[P],
       NumberForm[polyLength[P], {7, 5}],
       ScientificForm[circleResidual[P], 2],
       Style[If[simple, "yes", "NO"], If[simple, Darker[Green], Red], Bold]}},
     Frame -> All]}, Center]],
 {{nm, "teeth", "curve"}, demoNames},
 {{model, "TurningAngle", "curvature"}, models},
 {{dil, "ExactLength", "dilation"}, dilations},
 {{n, 30, "vertices (sampled curves)"}, 10, 200, 2, Appearance -> "Labeled",
  ContinuousAction -> False},
 {{red, False, "resample each step"}, {True, False}},
 {{uu, 0., "log time"}, 0., 1.},
 ControlPlacement -> Left, SaveDefinitions -> True]

The conjecture, live. The chord–arc cloud of the evolving polygon against the sliding barrier  — and the retitling is deliberate: Theorem 4.1 is a theorem about smooth curves, and nobody has proved a discrete Andrews–Bryan comparison for these polygonal flows. Every scrub of this demonstration is an experiment: the cloud staying above the barrier for a 10-gon, or for the teeth, is evidence for a conjecture, not an illustration of a theorem. (Section 13's open threads start exactly here.)

Manipulate[
 Module[{run, fr, t, ll, dd, tb},
  run = discRun[nm, n, "TurningAngle", "ExactLength", red];
  tb = demoTbar[nm, n];
  fr = frameOf[run, uu]; t = fr[[1]];
  {ll, dd} = chordArcPairs[fr[[2]]];
  Show[
   ListPlot[Transpose[{ll, dd}], PlotStyle -> {PointSize[0.008], GrayLevel[0.4]}],
   Plot[{fAB[x, t - tb], 2 Sin[x/2]}, {x, 0.001, Pi},
    PlotStyle -> {Directive[Thick, ColorData[97, 4]], Directive[Black, Dashed]}],
   PlotRange -> {{0, Pi}, {0, 2.1}}, Axes -> False, Frame -> True,
   FrameLabel -> {"arc ℓ", "chord d"}, ImageSize -> 460,
   PlotLabel -> Row[{"t = ", ScientificForm[t, 3], "    vertices ", Length[fr[[2]]],
      "    min(d - f) = ", ScientificForm[Min[dd - fAB[ll, t - tb]], 2]}]]],
 {{nm, "teeth", "curve"}, demoNames},
 {{n, 40, "vertices (sampled curves)"}, 10, 160, 2, Appearance -> "Labeled",
  ContinuousAction -> False},
 {{red, False, "resample each step"}, {True, False}},
 {{uu, 0., "log time"}, 0., 1.},
 ControlPlacement -> Left, SaveDefinitions -> True]

Smooth against discrete. The spectral solution (blue) overlaid with any subset of the discrete models at the same time, each recentred on its own centroid. One striking feature is not a bug: the shape distances of the different models agree to four digits or more (check the table), because the dominant error is the polygon-inscription sag  against the same smooth curve — identical for every model — while the models' genuine differences are an order smaller and show up only in the centre separation, where the flow's translation instability amplifies them exponentially.

Manipulate[
 Module[{run0, t, sm, cs, ds, rows},
  run0 = discRun[nm, n, First[ms /. {} -> {"TurningAngle"}], dil, True];
  t = frameOf[run0, uu][[1]];
  sm = smoothRun[nm]["Curve"][t, 300]; cs = Mean[sm];
  ds = Association@Table[m -> frameOf[discRun[nm, n, m, dil, True], uu], {m, ms}];
  rows = Table[Module[{P = ds[m][[2]], cp},
     cp = Mean[P];
     {Style["■", modelColors[m]], m,
      ScientificForm[Max[polylineDistance[Map[# - cs &, sm]] /@ Map[# - cp &, P]], 3],
      ScientificForm[Norm[cp - cs], 3]}], {m, ms}];
  Column[{
    Graphics[{smoothColor, AbsoluteThickness[2.2],
      Line[Append[#, First[#]] &@Map[# - cs &, sm]],
      Sequence @@ Table[Module[{P = ds[m][[2]], cp}, cp = Mean[P];
         {modelColors[m], AbsoluteThickness[1.3],
          Line[Append[#, First[#]] &@Map[# - cp &, P]],
          PointSize[0.010], Point[Map[# - cp &, P]]}], {m, ms}]},
     PlotRange -> {{-runExtent[run0], runExtent[run0]},
       {-runExtent[run0], runExtent[run0]}}, AspectRatio -> 1,
     ImageSize -> 430, Frame -> True, FrameTicks -> None,
     PlotLabel -> Style[Row[{"t = ", ScientificForm[t, 3]}], 11]],
    If[ms === {}, "select a model",
     Grid[Prepend[rows, Style[#, Italic] & /@
        {"", "model", "shape distance to smooth", "|centre - centre_smooth|"}],
      Frame -> All, Alignment -> Left]]}, Center]],
 {{nm, "lopsided", "curve"}, Keys[demoCurves]},
 {{n, 60, "vertices"}, 10, 160, 2, Appearance -> "Labeled", ContinuousAction -> False},
 {{ms, {"TurningAngle"}, "models"}, models, ControlType -> CheckboxBar},
 {{dil, "ExactLength", "dilation"}, dilations},
 {{uu, 0., "log time"}, 0., 1.},
 ControlPlacement -> Left, SaveDefinitions -> True]

13. Notes and open threads

1. The discrete monotone quantity. Is  non-increasing along any of the discrete flows, exactly? The smooth proof restarts Theorem 4.1, which has no discrete counterpart yet; numerically it looks true after the initial transient. This would be a genuinely structure-preserving discretization result, and it smells tractable, at least for the osculating-circle model whose ROUND property is exact.

2. A discrete distance comparison? Proposition 3 shows the osculating-circle translation gets the equality case (concyclic polygons) exactly right. Whether a full discrete analogue of  holds — with vertex pairs in place of point pairs — is open here; the first-variation structure (the equal-angle condition of Section 4.2 of the notes) makes sense verbatim for polygons.

3. The centre is not the story. All position information under this flow is exponentially unstable (Sections 4, 11), so any numerical claim about where the limit circle sits is meaningless without recentring; every claim about shape, length, curvature, and chords is robust. When restating the convergence theorem, "modulo translation" is doing real work — the parenthetical justification for translational convergence in Theorem 6.2 of the current lecture notes draft needs the extinction-point caveat.

4. Semi-implicit stepping would remove the  restriction: for the length-variation model  with current-configuration cotangent-free weights is the natural candidate, exactly as in the curve shortening notebook's notes. Not implemented.

5. Crane's no-free-lunch table (TOTAL/DRIFT/ROUND and the inequivalence of the four discrete curvatures) lives in the companion curve shortening flow notebook; this notebook is its sequel, and Proposition 3 is what its ROUND analysis becomes when the dilation term is added.

14. References

1. B. Andrews and P. Bryan, Curvature bound for curve shortening flow via distance comparison and a direct proof of Grayson's theorem, J. Reine Angew. Math. 653 (2011) 179–187; arXiv:0908.2682.

2. J. Cantarella, Round limits of the length-preserving curvature flow: the Andrews–Bryan argument, MATH 8230 lecture notes, 2026 (lecture_notes/andrews-bryan-grayson.pdf). All theorem numbers cited above refer to these notes.

3. M. Gage and R. S. Hamilton, The heat equation shrinking convex plane curves, J. Differential Geom. 23 (1986) 69–96.

4. M. A. Grayson, The heat equation shrinks embedded plane curves to round points, J. Differential Geom. 26 (1987) 285–314.

5. G. Huisken, A distance comparison principle for evolving curves, Asian J. Math. 2 (1998) 127–133.

6. B. Andrews, The affine curve-lengthening flow, J. Reine Angew. Math. 506 (1999) 43–83. (Source of the interpolation inequalities behind the  convergence.)

7. P. Bryan, M. Langford, and J. J. Zhu, Sharp distance comparison for curve shortening flow on the round sphere, arXiv:2310.02649. (Where the comparison-function method goes next.)

8. K. Crane, Discrete Differential Geometry: An Applied Introduction, CMU 15-458/858 course notes and slides, Lecture 1. (Source of the four discrete curvatures.)