Jim Propp once asked: "Is there a way to relax an approximation to a space-filling curve in continuous time so that it works out its kinks and regresses to simpler approximations? No interim self-intersections please!" Julian Ziegler Hunts developed a Fourier matrix product that can produce the analogous animation.

Julian wrote a nifty Fourier expander for recursive Koch polygons. E.g.,

ptsgnlst2Fouriermat[{0, 1, I^(2/3), 1 + I^(2/3), 2}, {1, -1, 1, -1}]

spacefills the triangle joining 0 to 2 via 1+i?3. Actually, it makes an infinite 3x3 matrix product for the coefficients a(k). Then Sum a(k+1/m) exp(2 i ? (k+1/m)) repeats the fractal on the sides of an m-gon. The GIF just accumulates consecutive harmonics. These are the definitions of functions:

ClearAll[ptsgnlst2Fouriermat]; ptsgnlst2Fouriermat[points_List, signs_List] := Block[{Amat, Bmat, M, klists = {(3 - signs)/2, Reverse[(signs + 3)/2]}, pt0 = points[[1]], pt1 = points[[-1]], n = Length[signs], m = 2, c,x}, Amat = {#, Reverse[#]} &[ signs*(Rest[points] - Most[points])/(pt1 - pt0)]; Bmat = {#, Reverse[#]} &[-signs*pt0/(pt1 - pt0) + points[[Range[n] + (1 - signs)/2]]]; Function[x, c] /. c -> Append[Table[Append[ Table[Exp[-I*x*(i - 1)/n]*Amat[[j, i]], {i, 1, n}].IdentityMatrix[m][[klists[[j]]]]/n, Sum[Exp[-I*x*(i - 1)/n]*Bmat[[j, i]], {i, 1, n}]*I*(Exp[-I*x/n] - 1)/x], {j, 1, m}], Table[Boole[i == m], {i, 0, m}]]]; ClearAll[mat2series]; mat2series[var_, mat_, n_, m_: 1, terms_, prodterms_, prec_: 60] := Block[{A, x}, A[0] = Limit[mat[x], x -> 0]; A[x_] = mat[x]; Exp[2*I*\[Pi]*Range[-terms + 1/m, terms + 1/m]* var].Table[(Dot @@ Table[A[N[2*\[Pi]*(k + 1/m)/n^r, prec]], {r, 0, prodterms}])[[1 ;; -2, -1]], {k, -terms, terms}]]

Here is an example of 3-fold rotational symmetry:

Block[{start = 0, end = 1, m = -3., f = mat2series[t, %, 7, -3., 33, MachinePrecision][[1]]}, ParametricPlot[ReIm[f + Sum[ E^(2 \[Pi] (k + 1/m) I t) I*(start - Exp[-2*I*\[Pi]/m]*end)/ 2/\[Pi]/(k + 1/m), {k, -33, 33}]], {t, 0, m}, Axes -> False]]

And here it is with higher details:

Block[{start = 0, end = 1, m = -3., f = mat2series[t, %%, 7, -3., 69, MachinePrecision][[1]]}, ParametricPlot[ReIm[f + Sum[ E^(2 \[Pi] (k + 1/m) I t) I*(start - Exp[-2*I*\[Pi]/m]*end)/ 2/\[Pi]/(k + 1/m), {k, -69, 69}]], {t, 0, m}, Axes -> False]]

And given enough time the code below will produce the animation at the top of this post:

Table[Block[{start = 0, end = 2, m = 6., f = mat2series[t, ptsgnlst2Fouriermat[{0, 1, I^(2/3), 1 + I^(2/3), 2}, {1, -1, 1, -1}], 4, m, n, MachinePrecision][[1]]}, ParametricPlot[ReIm[f + Sum[E^(2 \[Pi] (k + 1/m) I t) I*(start - Exp[-2*I*\[Pi]/m]*end)/ 2/\[Pi]/(k + 1/m), {k, -n, n}]], {t, 0, 6}, Axes -> False]], {n, 99}]; Export["/Users/billgosper/Movies/hellodoily.gif", Join[ConstantArray[%[[1]], 6], %, ConstantArray[%[[-1]], 9]]]