随机多边形转化为椭圆的过程研究
一、问题背景¶
1.1 背景来源¶
在网上看到一篇由 Adam N. Elmachtoub 和 Charles F. Van Loan 合作撰写的论文From Random Polygon to Ellipse: An Eigenanalysis 后，我对这篇文章的结论非常感兴趣。由于原文主要是对问题进行了理论上的证明，而没有对这个过程进行很好的展示，恰巧小学期的Mathematica课程要做个报告作业，于是打算利用Mathematica对这个问题进行具体的仿真模拟研究，使得对这个问题的理解更为透彻。
1.2 问题描述¶
在二维实空间中随机生成\(N\)个点，将它们按生成顺序首尾相连，构成一个封闭的\(N\)边多边形。将这个随机生成的多边形的\(N\)条边的中点按边的生成顺序相连，得到新的封闭多边形。重新对生成的封闭多边形进行相同的操作，并不断重复上述的迭代过程，在经过足够次数的变换后(更具体地讲，记\(n\)为操作次数，则要求\(n\to\infty\))，最终这个封闭图形一定会是一个椭圆。
1.3 问题的初步理解¶
先考虑最简单的情形——三角形，如果是对二维空间中不共线的三个点进行上述操作，可以想象，最终会逐步变成一个非常小的三角形，它甚至会无限趋近于同一个点，显然，这个点可以被理解为一个非常小的圆形，也即是一种特殊的椭圆。
当然，一个很小的点我们可以认为是很小的椭圆，但甚至也可能认为是很小的矩形或者其他图形。所以，上面的话姑且只能当做一个用处不大的“猜想验证”。
而当情况变为一个任意的多边形时，至少可以很容易地理解，这些逐步形成的多边形的边会越来越近似于一条光滑的曲线。 从上述描述可以看出，我们通过简单的想象，至少能够初步对1.2中所提出结论有一个初步的直观感受，而至于结论是否成立，则需要通过理性的证明以及合理的验证来得到结论。
二、问题的数学证明思路¶
【说明】原文中的数学证明较为冗长，且因本报告的重点内容在于结论的过程仿真，即结论的验证过程，故只对问题的数学证明过程作一个较为简略的介绍。
首先，将取中点的操作抽象为矩阵的变换，记该变换的对应矩阵为\(T\)，被变换的矩阵是一个\(N\times 2\)阶矩阵，矩阵是由\(N\)个二维点坐标按行排列得到的，记为\(A\)，则\(A \times T\)可以抽象地表示原图形初次操作后的结果，进一步地，不难理解，\(A\times T^k\)可以抽象地表示原图形经过\(k\)次变换后的图形。
所需要进一步理解的就是这个变换矩阵的\(k\)次幂，事实上，这样的矩阵是一个循环矩阵(把一个行向量循环右移一位，作为第二行，再循环右移一位，作为第三行，以此类推，最终所得到的矩阵即为循环矩阵了)。
记上述所得到的循环矩阵为\(M\)，根据线性代数中的相关定理可知，\(M\)可以表示成形如\(P \Lambda P^{-1}\)的形式，其中\(\Lambda\)是由\(M\)的特征值构成的对角矩阵，则\(M^k=P \Lambda^k P^{-1}\)。
在原文中，作者引入一个循环矩阵的性质，即循环矩阵的特征值就是矩阵首行的傅里叶变换。接下来将傅里叶变换矩阵代入，并且发现当\(k\to\infty\)的时候，只有一项不为\(0\)，最终对该项进行验算发现满足椭圆性质，进而得到了原命题，即通过足够多次的变换后，原多边形可以变换为一个椭圆。
三、随机生成的多边形转化为椭圆过程的Mathematica仿真研究¶
3.1 随机多边形的生成¶
首先随机生成30组二元数组，用来表示二维平面上的30个点。
n = 30;
PointList = RandomReal[10, {n, 2}]
以二元坐标形式记录生成的点，如下
{{3.80897, 8.86568}, {5.16511, 1.04105}, {2.20855, 5.62259},
{9.999, 7.61291}, {9.93663, 1.7641}, {3.70494, 5.25041},
{6.99509, 8.61465}, {7.74337, 7.56245}, {3.75603, 6.21289},
{4.83932, 8.79755}, {9.42899, 2.48778}, {8.10076, 2.56023},
{8.19965, 3.89118}, {1.88996, 0.200388}, {7.17207, 2.40234},
{8.44527, 5.75125}, {2.24597, 6.67419}, {9.44965, 0.763671},
{0.748501, 0.54816}, {4.56093, 3.60478}, {5.5144, 6.29168},
{7.9264, 7.02054}, {1.30604, 5.32975}, {6.03109, 8.22584},
{4.31784, 4.02132}, {0.776391, 4.3736}, {1.17457, 7.57537},
{5.60569, 0.112303}, {2.76632, 0.70948}, {4.2963, 3.65391}}
然后根据这30个点绘出其对应的随机多边形，用以初步观察其性质。
Graphics[Polygon[PointList]]
3.2 绘图函数的构建¶
为了能够在迭代环境下进行图形绘制，我们必须要构建出合适的绘图函数。
下面的绘图函数中，pt参量指点坐标集，而next参量指迭代后的后续点坐标，通过不断地向pt中加入next点，我们可以使用该函数针对本问题在后续中画出相应图形。
在该函数中，使用了具体的点格式参数，便于后续观察。
Paint[pt_] :=
Module[{next = pt},AppendTo[next, pt[[1]]];
ListLinePlot[next, Epilog -> { PointSize[Large], Orange, Point[next]}]];
3.3 功能函数的构建¶
首先，我们要构建出取中点的功能函数。
在该函数中，以centerpt为中点参量，通过表达式centerpt = (pt + RotateLeft[pt, 1])/2确保centerpt能够不断地在迭代过程中取得各边的中点。
GetCenterList[pt_] :=
Module[{centerpt},centerpt = ((pt + RotateLeft[pt, 1])/2)];
接下来即是程序设计中最为核心的一个部分，即为点的迭代过程，通过看似简单的语句，它能够将前面几个函数巧妙地联系在一起，通过迭代过程达到我们的预期目的。
其中，需要指出的是，loop参量是指迭代次数也即循环次数，根据前面的描述可知，理论上这个数值设定得越大越好，但经过本人的实际测试，经过几百次的迭代后所得到的图形已经足够近似于椭圆，如果将数值设定得更大则意义不大，反倒会减慢程序的运行速度。
loop = 300;
iter = NestList[GetCenterList, PointList, loop];
3.4 程序调试¶
上述的工作都是大部分的准备工作，接下来我们需要利用上述成果，重新切入本次报告的核心主题——研究随机多边形转化为椭圆的动态过程。因此，我们需要对上述的一些函数进行进一步的调试使用，进而达到观察这个过程的目的。
首先，我们需要观察一下整个过程的图像，但由于数量太多，我们等间隔地取出一部分来作观察。
Table[Paint[iter[[i]]], {i, 1, loop, 20}]
从上面的图像可以初步得出下面几点结论：
- 该算法确实能够较好地将随机多边形转化为椭圆；
- 在转化过程的初期，由于点的间隔较大，点之间较为分散，可以看出转换速度非常快，而后期的转换速度较为缓慢；
- 虽然直观地从图像来看，各个图像之间似乎是相同大小的，但实际从坐标来看，这个图形确是在逐步地变小。这一点很好理解，因为对于封闭图形来说，不断地连接边的中点形成新的图形显然会使得图形越来越小。
由于取出部分图形来观察的方法还不够直观，在验证完初步的猜想后，我们进一步地用动态画面来观察这一过程，使该过程更完整地展现出来。
plotlist = Table[Paint[iter[[i]]], {i, 1, loop}];
ListAnimate[plotlist]
- 这个图看上去好厉害的样子，是不是？:)
- 有强迫症的同学看到后面请一定保持耐心。。。
3.5 探究结论¶
根据转化过程的过程图像，并结合上述的动画，我们完整地观察到了随机生成多边形转化为椭圆的动态过程，通过Mathematica这一数学软件，我们成功验证了前文所提出的问题，并弥补了理论证明在直观性上的不足，使得这个结论更加直观可感，也达到了之前的研究目的。
四、参考资料¶
- From Random Polygon to Ellipse: An Eigenanalysis.Adam N. Elmachtoub & Charles F. Van Loan.Cornell University.
- 随机多边形化为椭圆的数学证明.刘孟白.(链接)
五、后记¶
这篇文章可以说技术含量并不高，因为只是一篇简单的作业文，说得狠一点就是篇“水文”。 而由于自己数学能力不足且对Mathematica的掌握还不够，所以可以肯定的是上面的东西，绝对做得还不够，希望能多通过这样的小“研究”来帮助提升自己的能力。