离散平均曲率流的一种数值模拟

给定一个$n$-多边形, 假设其顶点满足方程
\[
\dot v_i(t)=\frac{\nu_i(t)}{\| \nu_i(t) \|^2},
\]
其中$\nu_i(t)$是顶点$v_{i-1},v_i,v_{i+1}$构成的三角形之外接圆心。它可以视为连续情形下的曲线平均曲率流的一种离散推广。

我们知道连续情形下，平均曲率流有所谓的Gage-Hamilton-Grayson定理，它表明平均曲率流保持简单曲线为简单曲线。

但下面的数值模拟表面，这个平均曲率流不一定保持曲线的简单性。

\documentclass[border=5mm]{standalone}
\usepackage{luamplib}
\begin{document}
\mplibtextextlabel{enable}
\mplibnumbersystem{double}
\begin{mplibcode}

u:=10pt;
vardef DiscreteFlow(suffix $)(text Pairs) = 
  save i_; i_ := 0; %save makes i_ local
    pair $[]; %define a pair with the variable name $
    %Typical loop over passed elements
    for i = Pairs: 
      $[i_] := i*u;
      i_ := i_ + 1;
      show $[i_]; 
    endfor;
enddef;
%DiscreteFlow(B)((1,2),(3,4),(5,5),(-2,-3));
%draw B[0] for i=1 upto 3: --B[i] endfor --cycle;
 
vardef CircumCenter(expr A,B,C)= 
  save p; pair p;
  p=.5[A,B]+whatever*((A-B) rotated 90)=.5[A,C]+whatever*((A-C) rotated 90);
  p
enddef;
vardef NormSquare(expr A) = 
  xpart(A)**2+ypart(A)**2
enddef;
 
 
 
beginfig(1);
  %pair B[];
  Maxiteration=50; delta=1; N=20; %number of edges
  for j=1 upto N:
    randomseed := uniformdeviate infinity;
    z[j]=8*dir(360j/N)+(normaldeviate,normaldeviate);
    show "z[" & decimal j & "]= (" & decimal xpart(z[j]) & "," & decimal ypart(z[j]) & ")";
  endfor;
 
  z0=z[N];
 
  draw (z[1] for i=2 upto N: --z[i] endfor --cycle) scaled 10u  withcolor blue withpen pencircle scaled 3pt;
  dotlabel("c2", CircumCenter(z[1],z[2],z[3])*10u);
  for i=1 upto N:
    dotlabel(decimal i, z[i] scaled 10u);
  endfor;
  pair p[], c;
  %initial
  for i=0 upto N: p[i]=z[i]; endfor;
  %iteration
  for n=1 upto Maxiteration:
    for i=1 upto N:
      c:=CircumCenter(p[(N*(n-1)+i-1) mod N],p[(N*(n-1)+i) mod N],p[(N*(n-1)+i+1) mod N]);
      p[N*n+i]=p[N*(n-1)+i]+delta*c/NormSquare(c);
      %label(decimal (N*n+i), p[N*n+i]*10u);
	drawarrow (p[N*(n-1)+i]--p[N*n+i]) scaled 10u withpen pencircle scaled .3pt;
    endfor;
      if n=Maxiteration:
	draw (p[N*n+1] for i=2 upto N: -- p[N*n+i] endfor --cycle) scaled 10.1u withpen pencircle scaled 1pt withcolor red;
      elseif (n mod 2=0):
	draw (p[N*n+1] for i=2 upto N: -- p[N*n+i] endfor --cycle) scaled 10u;
      else:
      draw (p[N*n+1] for i=2 upto N: -- p[N*n+i] endfor --cycle) scaled 10u  dashed evenly withcolor .5*red;
      fi;
  endfor;
 
 
 
endfig;

\end{mplibcode}
\end{document}

一个输出为，给定顶点如下

(luamplib)               {randomseed:=164.13048854794684}
(luamplib)               >> "z[1]= (7.1663813354355792,2.115410943910577)"
(luamplib)               {randomseed:=1282.7618269274324}
(luamplib)               >> "z[2]= (6.1198547032965198,3.1294415485934794)"
(luamplib)               {randomseed:=403.39625352345581}
(luamplib)               >> "z[3]= (5.3874308800111868,4.8324648551767311)"
(luamplib)               {randomseed:=3339.1733235391926}
(luamplib)               >> "z[4]= (3.5465950228343885,6.1858621289233282)"
(luamplib)               {randomseed:=1839.1381135784272}
(luamplib)               >> "z[5]= (1.2672314185221965,7.8164049682362586)"
(luamplib)               {randomseed:=3420.0974717036256}
(luamplib)               >> "z[6]= (-3.7256314873696939,6.8671219483403938)"
(luamplib)               {randomseed:=510.18410242830026}
(luamplib)               >> "z[7]= (-5.6650842541666693,6.4155586555366675)"
(luamplib)               {randomseed:=12.746562895524596}
(luamplib)               >> "z[8]= (-7.3959720323999152,5.2194560696781931)"
(luamplib)               {randomseed:=1303.1675809171504}
(luamplib)               >> "z[9]= (-6.7246625396621438,1.5797575825552226)"
(luamplib)               {randomseed:=1827.3952890813218}
(luamplib)               >> "z[10]= (-9.1673328058487815,-0.45325023466673714)"
(luamplib)               {randomseed:=7.8023681259649997}
(luamplib)               >> "z[11]= (-9.510793507129911,-3.9190636005348356)"
(luamplib)               {randomseed:=1418.1949241054936}
(luamplib)               >> "z[12]= (-7.425895076183326,-4.2579680419470947)"
(luamplib)               {randomseed:=1480.3974041179326}
(luamplib)               >> "z[13]= (-3.3223499681157214,-6.4744847060536879)"
(luamplib)               {randomseed:=876.68583251129962}
(luamplib)               >> "z[14]= (-1.8601946468263932,-6.6736668112858331)"
(luamplib)               {randomseed:=2646.6707481852013}
(luamplib)               >> "z[15]= (-0.45124071673288363,-9.4417326012844018)"
(luamplib)               {randomseed:=132.57237941859191}
(luamplib)               >> "z[16]= (2.9170180427407999,-8.5078926786590099)"
(luamplib)               {randomseed:=1047.8092142440955}
(luamplib)               >> "z[17]= (4.2494809218946745,-6.7105559485829076)"
(luamplib)               {randomseed:=2021.5000778674312}
(luamplib)               >> "z[18]= (8.0277937535304655,-3.4644864211184609)"
(luamplib)               {randomseed:=2055.7357077356655}
(luamplib)               >> "z[19]= (10.187937026289996,-3.2670946608207374)"
(luamplib)               {randomseed:=3756.315583007492}
(luamplib)               >> "z[20]= (8.0942754644586294,-0.78244053305555517)"

蓝色多边形经过演化后，变成红色多边形

Leave a Reply Cancel reply