跳转至

随机多边形转化为椭圆的过程研究

——基于Mathematica软件进行仿真模拟

一、问题背景

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]]

random-polygon.png

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}]

multi-ellipse.png

从上面的图像可以初步得出下面几点结论:

  • 该算法确实能够较好地将随机多边形转化为椭圆;
  • 在转化过程的初期,由于点的间隔较大,点之间较为分散,可以看出转换速度非常快,而后期的转换速度较为缓慢;
  • 虽然直观地从图像来看,各个图像之间似乎是相同大小的,但实际从坐标来看,这个图形确是在逐步地变小。这一点很好理解,因为对于封闭图形来说,不断地连接边的中点形成新的图形显然会使得图形越来越小。

由于取出部分图形来观察的方法还不够直观,在验证完初步的猜想后,我们进一步地用动态画面来观察这一过程,使该过程更完整地展现出来。

plotlist = Table[Paint[iter[[i]]], {i, 1, loop}];
ListAnimate[plotlist]

polygon-to-ellipse.gif

  • 这个图看上去好厉害的样子,是不是?:)
  • 有强迫症的同学看到后面请一定保持耐心。。。

3.5 探究结论

根据转化过程的过程图像,并结合上述的动画,我们完整地观察到了随机生成多边形转化为椭圆的动态过程,通过Mathematica这一数学软件,我们成功验证了前文所提出的问题,并弥补了理论证明在直观性上的不足,使得这个结论更加直观可感,也达到了之前的研究目的。

四、参考资料

  1. From Random Polygon to Ellipse: An Eigenanalysis.Adam N. Elmachtoub & Charles F. Van Loan.Cornell University.
  2. 随机多边形化为椭圆的数学证明.刘孟白.(链接)

五、后记

这篇文章可以说技术含量并不高,因为只是一篇简单的作业文,说得狠一点就是篇“水文”。 而由于自己数学能力不足且对Mathematica的掌握还不够,所以可以肯定的是上面的东西,绝对做得还不够,希望能多通过这样的小“研究”来帮助提升自己的能力。