均分田地问题

Sat, 19th October 2019Edit on Githubbubbles

绪论

均分田地问题中,gxqcn提出了如下问题

有一块田地需要分给n户家庭,要求各户分得面积都相等。田地内部不同家庭分得的区域将建田埂以分隔(原待分地已有田埂圈定)。现为实现耕地面积最大化,要求新建田埂总长度最小。请问如何规划?

关于这个问题的三维问题,存在Plateau's laws

在这块田地是单位正方形,对于n=3,4,5, KeyTo9_Fans分别给出了很不错的结果:

n=3时, 田埂总长度可以仅为23+π6+3=1.623278144157\frac23+\frac{\pi}6+\sqrt{3}=1.623278144157 (我们可以搜索到法语论坛相关讨论得出了相同的结果) s3c

n=4时,田埂总长度可以仅为2+13+π3=1.975592884782\sqrt{2}+\sqrt{1-\sqrt{3}+\frac{\pi}3}=1.975592884782 (我们可以搜索到意大利语相关讨论得出了相同的结果) s4c

n=5时, 田埂总长度可以为2+((1+3)π6)215((2+3)π3(1+3))=2.5021129304272+((1+\sqrt{3})\pi-6) \sqrt{\frac{2}{15((2+\sqrt{3})\pi-3(1+\sqrt{3}))}}=2.502112930427 (可以网络找到Martin Gardner包含相同的结果) s5c

对于单位圆形田地,在n=3时,从圆心发出三条两两夹角互为120°的半径就可以等分面积,总长度为3。

在n=4时,我们给出了田埂总长为3.945702967267的方案 c4

在n=5的圆形田地,我们找出了田埂总长为4.833846643527的方案 c5

在n=6的圆形田地,我们找出了田埂总长为5.406796929952的对称方案 c6

在n=7的圆形田地zgg__和数学星空指出在田地正中心使用正六边形,每个顶点向圆周引垂线可得田埂总长为6的方案

在帖子中,很早就有人发现每段田埂应该是直线段或者圆弧;在区域内部,三条田埂相交于一点时,交点处切向量两两夹角为120°;而且田埂和田地边界相遇时在相遇点的切线要和田地边界在相遇点的切线垂直。后文中,我们会直接用俩田埂的夹角为120°来表示俩田埂相遇于一个交点,而且在交点处的切向量夹角为120°;同样我们会直接用田埂和田地边界垂直代表田埂和田地边界相遇并且在相遇点两者切线相互垂直。

KeyTo9_Fans最早指出田埂两两夹角为120°。他的思路是如下图,在仅考虑三条田埂相交于点P的情况,分别在三条田埂上各自选择充分接近P点的A、B、C三点。为了达到田埂总长最小值,P点应该移动到三角形ABC的费马点才会让AP+BP+CP取到最小值,由此得出AP、BP、CP两两夹角应该为120°。但是这里稍微有点不严密,因为图中这种移动P到P’的方式会改变周围三个区域的面积,而我们的问题中要求所有区域的面积要求是相等的。 fermat

zgg__最早指出,本问题等价于二维肥皂泡稳态问题,由此得出每段田埂都必须是圆弧或直线段,而且两两夹角为120°;hujunhua进一步指出利用肥皂泡模型可以说明和田地边界相遇的田埂必然垂直于田地边界;gxqcn紧接着指出在区域数目n充分大时,内部区域应该出现“蜂窝”状。
我们对上面情况进行总结,并且尝试用数学方法进行证明:

i)田埂每一段必然是圆弧或直线段(可以看成圆弧的退化情况)。
ii)在内部,最多三个不同的田埂共点,这时田埂之间必然两两夹角相等
iii)田埂和边界接触的地方必然同边界垂直

其中性质i)可以直接用初等方法证明。我们可以利用一个大家所熟知的结论:面积一定的简单平面图形中,圆的周长最短。利用这个结论可以得出一侧为固定线段,另一侧为任意简单曲线围成的面积一定的图形中,在另一侧任意简单曲线为圆弧时,这条简单曲线长度达到最短。

arc 如上图所示,对于围住的面积固定的曲线BAC和固定直线段BC,我们先以线段BC为弦作出和这个图形等面积的弓形BA’C并且补充弓形BA”C形成一个完整的圆BA’CA”。于是图形BACA”和图形BA’CA”等面积,于是我们根据等面积图形中圆的周长最小可以得出(BA"C) ̂+(BA'C) ̂≤(BA"C) ̂+(BAC),所以我们可以得出(BA'C) ̂≤(BAC),得出结论i)。 于是我们得出极值条件所有田埂都是圆弧或直线段(直线段是圆弧的退化情况)。 后面我们均只需要考虑所有田埂都是直线段或圆弧的情况。

对于其中任意一个三条田埂汇聚于一点的P(xp,yp)P(x_p,y_p )的情况,假设这三段田埂另一个端点分别为A(xa,ya),B(xb,yb),C(xc,yc)A(x_a,y_a ),B(x_b,y_b ),C(x_c,y_c), 而三段圆弧和对应端点连线的有向夹角分别为θa,θb,θc\theta_a,\theta_b,\theta_c (对应三段圆弧的弧度分别为2θa,2θb,2θc2\theta_a,2\theta_b,2\theta_c)。 其中圆弧关于P点向顺时针方向凸出(也就是对应圆心落在对应直线段的逆时针位置时)对应的θs\theta_s看成是正角度;而圆弧逆顺时针凸出时,对应的θs\theta_s看成是负角度,比如图中所示θa,θb\theta_a,\theta_b为负角度,θc\theta_c为正角度。另外需要注意的是,一段田埂在一个端点处的角度为θ\theta,那么在另外一个端点处的角度必然是θ-\theta

basic

使用拉格朗日乘数法,选择A、B、C三点固定不动,仅让P点保持运动并且让P点周围三个区域的面积保持不变,要求三条田埂长度之和最小,可以得出一些约束方程, 最后可以把极值约束条件化简{F1=sin(θa)La+sin(θb)Lb+sin(θc)Lc=0F2=Da,x+Db,x+Dc,x=0F3=Da,y+Db,y+Dc,y=0\begin{cases} F_1=\frac{\sin⁡(\theta_a)}{L_a} +\frac{\sin⁡(\theta_b)}{L_b} +\frac{\sin⁡(\theta_c)}{L_c} =0\\ F_2=D_{a,x}+D_{b,x}+D_{c,x}=0\\ F_3=D_{a,y}+D_{b,y}+D_{c,y}=0 \end{cases} (I)

其中 {Ds,x=xsxpLscos(θs)+ysypLssin(θs)Ds,y=ysypLscos(θs)xsxpLssin(θs)Ls=(xpxs)2+(ypys)2sa,b,c\begin{cases} D_{s,x}=\frac{x_s-x_p}{L_s}\cos⁡(\theta_s)+\frac{y_s-y_p}{L_s}\sin⁡(\theta_s )\\ D_{s,y}=\frac{y_s-y_p}{L_s}\cos⁡(\theta_s)-\frac{x_s-x_p}{L_s}\sin⁡(\theta_s )\\ L_s=\sqrt{(x_p-x_s )^2+(y_p-y_s )^2 }\\ s\in \text{{a,b,c}} \end{cases}

上面约束条件(I)中第一条方程代表三段圆弧(或直线段)有向曲率之和为0。而后面两条方程代表三段圆弧在P点处单位切线向量之和为零向量,所以它们两两夹角都相等,均为120°(P点同时为三个单位切线向量末端构成三角形的外心和重心,所以这是一个正三角形)。

对于边界为直线段构成的凸田地,可以想象如果某个田埂的一个端点落在田地边界上,那么将整个田地和所有田埂关于这个边界线做对称图形,然后把原图形和对称图形合并在一起并且将仅被这条分界直线段分隔的区域相互合并,于是变成一个关于更大的田地更多区域的固定面积最优化问题,结果同样需要满足类似的约束条件。于是我们可以看出在原问题中,这个田埂必须垂直于田地边界才符合要求。

另外我们推导了田地边界是圆弧构成的情况,同样得出了和边界相遇的田埂必须垂直于圆弧田地边界的条件。

boundary 如图,不妨设田地边界这一段圆弧为单位圆的圆弧,设田埂在边界上点为U(cos(u),sin(u))U(\cos⁡(u),\sin⁡(u) ),田埂另外一个端点为点P(x0,y0)P(x_0,y_0 ),圆弧和俩端点连线夹角为θ\theta。并且设P点极坐标为(r0,θ0)(r_0,θ_0 ), 得出的约束条件为 tan(θ)r0sin(uθ0)1r0cos(uθ0)=0\tan⁡(\theta)-\frac{r_0\sin⁡(u-\theta_0 )}{1-r_0\cos⁡(u-\theta_0 )}=0
E1=sin(θ)(1cos(u)x0sin(u)y0)cos(θ)(sin(u)x0cos(u)y0)=0E_1=\sin(\theta)(1-\cos⁡(u) x_0-\sin⁡(u) y_0 )-\cos⁡(θ)(\sin⁡(u) x_0-\cos⁡(u) y_0)=0 (II)

这个条件正好就是要求田埂(UP) ̂垂直单位圆田地边界于U点。

利用上面分析结论我们使用计算机将正方形田地均分到最多14个区域并且将圆形田地均分到最多19个区域. 但是很快我们发现英国数学家Simon Cox已经计算到42个区域以内的非常不错的结果。

后来我们利用pari/gp计算出了圆形情况不超过43个区域下的50位精度的高精度结果。 点击链接可以下载关键点的坐标和圆弧角度的高精度结果。 另外我们还继续计算得出将圆形划分44~52个区域时的不错结果,如下面是将圆形田地划分为44个区域的划分方案图: c44

第一部分 理论分析

现在我们可以把这个问题重新描绘为:

定理1 给定边界为最多包含有限个不光滑点的简单曲线围成的区域,在其内部通过简单曲线分割成n个分别给定面积的小区域,那么在这些分割线长度之和达到极小值时,必然满足一下条件:

i)每段分割线必然是圆弧或直线段

ii)如果两段曲率不同的分割线相遇于区域内部一点P,那么必然还有第三段分割线同样和这两段分割线相遇于P点

iii)如果两条分割线相遇于一点P,那么它们的夹角不小于120°。特别的,如果P点在区域内部,那么必然是三条分割线在端点相遇于P点,而且两两夹角正好是120°;而且不会有四条或以上的分割线在区域内部相遇于一点。

iv)如果一条分割线相遇区域边界于点P,那么分割线和边界在点P两个方向的切向量夹角都不小于90°。特别的,如果区域边界在P点光滑,那么分割线这时必然垂直区域边界于P点。

v)如果有三条分割线相遇于区域内部一点P,那么这三条分割线的有向曲率和为0。

证明如下:

其中性质i)在绪论中通过反证法已经给出了比较严密的初等证明。

对于性质ii),我们同样采用类似的反证法,假设这时没有第三条分割线经过P点,由于P点是区域内点,这两条分割线两侧属于两个不同的小区域。在两条分割线上分别选择和P点充分接近的A、B两点,类似绪论中的方法,把非圆弧APB替换成一段圆弧并且保持两边小区域面积不变可以使得分割线总长度减少。由于P时区域内点,在A、B、P充分接近时,这种替换不会和其它分割线相交,所以可以得到一个合法的分割线总长度更小的方案,说明这种方法必然无法达到极小值。

对于性质iii), 同样可以采用反证法(下面的结果暂时还需要依赖数值计算结果)。

fermat2

如图假设这时有两条分割线AB,AC相遇于一点A而且AB,AC是同一个小区域的两条相邻的边。这个小区域分别还隔着AB或BC和另外两个小区域相邻。选择分割线上的点B、C使得AB=AC而且点B、C充分接近A点,在局部这时可以近似把AB、AC看成直线段,这时我们看看能否找到两条夹角为120°的圆弧(BD) ̂和(DC) ̂使得添加点D并且用(BD) ̂、(DC) ̂、(AD) ̅替换AB、AC后保持三个周边小区域的面积不变。如果这时有(BD) ̂+(AD) ̅+(DC) ̂<(AB) ̅+(AC) ̅,那么我们就可以通过这种替换保持各个小区域面积不变的同时减少了分割线总长度,从而得出原先的方案不是极小值方案。在∠BAC<120°时,我们让D在BC中垂线上从Fermat点移动到A点,保持圆弧BD和AD夹角为120°。当D移动到Fermat点时,显然BDC围成的小区域面积会偏小(另外两个区域的面积都增加了),而当D移动到A点时,(BD) ̂和(DC) ̂分别移动到直线段AB和AC的外侧,所以BDC围成的小区域面积会偏大(另外两个区域面积都减少了),根据中值定理,这中间必然正好存在一个D点,使得三个区域面积都和原先相等。

anglecurve

如图设BAD=θ,AEB=2ψ\angle BAD=\theta, \angle AEB=2\psi,其中E为圆弧(BD) ̂的圆心,经计算可以得到关系式ctan(θ)=2ψsin(2ψ)+2sin(2π32ψ)sin2(ψ)4sin2(ψ)sin3(π3ψ)=2ψsin(2ψ)4sin2(ψ)sin3(π3ψ)+cos(π3ψ)sin2(π3ψ)=u(ψ)+v(ψ)ctan(\theta)=\frac{2\psi-\sin⁡(2\psi)+2\sin⁡(\frac{2\pi}3-2\psi)\sin^2⁡(\psi)}{4\sin^2⁡(\psi)\sin^3⁡(\frac{\pi}3-\psi)}=\frac{2\psi-\sin⁡(2\psi)}{4\sin^2⁡(\psi)\sin^3⁡(\frac{\pi}3-\psi)}+\frac{\cos⁡(\frac{\pi}3-\psi)}{\sin^2⁡(\frac{\pi}3-\psi)} =u(\psi)+v(\psi),图上红色曲线给出了0<θ<π3,0<ψ<π30\lt \theta\lt\frac{\pi}3,0\lt\psi\lt\frac{\pi}3θ,ψ\theta,\psi的关系图,其中横坐标为θ\theta,纵坐标为ψ\psi,另外绿色曲线给出((AB) ̅+(AC) ̅)-((BD) ̂+(AD) ̅+(DC) ̂)和(AB) ̅的比例22ψsin(θ)sin(ψ)sin(2π3+ψ)sin(π3θψ)sin(2π3+ψ)2-\frac{2\psi\sin⁡(\theta)}{\sin⁡(\psi)\sin⁡(\frac{2\pi}3+\psi)}-\frac{\sin⁡(\frac{\pi}3-\theta-\psi)}{\sin⁡(\frac{2\pi}3+ψ)}的图像,可以看出这种情况进行替换总是可以缩短分割线总长度的。

为了证明上图中的性质,我们计算dv(ψ)dψ=sec(ψ+π6)(tan2(ψ+π6)+sec2(ψ+π6))0\frac{dv(\psi)}{d\psi}=\sec⁡(\psi+\frac{\pi}6)(\tan^2⁡(\psi+\frac{\pi}6)+\sec^2⁡(\psi+\frac{\pi}6))\ge0, du(ψ)dψ=14csc2(ψ)sec3(ψ+π6)(4(1ψtan(ψ))+3(2ψsin(2ψ))tan(ψ+π6))0\frac{du(\psi)}{d\psi}=\frac14\csc^2(\psi)\sec^3⁡(\psi+\frac{\pi}6)(4(1-\frac{\psi}{\tan⁡(\psi)})+3(2\psi-\sin⁡(2\psi))\tan⁡(\psi+\frac{\pi}6))\ge 0,

所以我们确认了θ,ψ\theta,\psi是递减关系。

另外我们还需要证明22ψsin(θ)sin(ψ)sin(2π/3+ψ)sin(π3θψsin(2π3+ψ)02-\frac{2\psi\sin⁡(\theta)}{\sin⁡(\psi)\sin⁡(2\pi/3+\psi)}-\frac{\sin⁡(\frac{\pi}3-\theta-\psi}{\sin⁡(\frac{2\pi}3+\psi)} \ge 02ψcos(π3ψ)sin(ψ)sin(ψ)sin(π3ψ)2cos(θ)sin(θ)\frac{2\psi-\cos⁡(\frac{\pi}3-\psi)\sin⁡(\psi)}{\sin⁡(\psi)\sin⁡(\frac{\pi}3-\psi)}\le \frac{2-\cos⁡(\theta)}{\sin⁡(\theta)}。待解决!

boundarymin

对于性质iv),同样假设现在有一条分割线AP和区域边界PC相遇于P点而且θ=∠CPA<90°,其中P点局部放大以后我们可以近似认为PC和PA都接近直线段。我们现在试着用垂直PC的一段圆弧(AD) ̂替换AP使得替换后保持两边小区域的面积不变。设∠ACP=ψ,我们得到条件sin(ψ)ψ=sin(θ)sin(θ+ψ)\frac{\sin⁡(\psi)}{\psi}=\frac{\sin⁡(\theta)}{\sin⁡(\theta+\psi)} ,而(AD) ̂/(AP) ̅ =(ψ sin⁡(θ))/sin⁡(ψ) =sin⁡(θ+ψ)≤1。其中去等号时要求sin⁡(ψ)/ψ=sin⁡(θ)=cos⁡(ψ)即tan⁡(ψ)=ψ只有角度为0时才可能。所以我们同样得出这时不是极小值情况。

对于性质v),对于三条分割线交于区域内部一点P(xp,yp)P(x_p,y_p ), 假设我们选择不调整点P和三条分割线另外一个端点A(xa,ya),B(xb,yb),C(xc,yc)A(x_a,y_a ),B(x_b,y_b ),C(x_c,y_c )的位置,但是允许调整三条分割线的曲率但是要求保持周边三个小区域的面积不变。直接采用链接拉格朗日乘数法中的方法,记LaL_a=(AP) ̅,LbL_b=(BP) ̅,LcL_c=(CP) ̅, basic 最后可以得出取极小值时的必要条件是sin(θa)La+sin(θb)Lb+sin(θc)Lc=0\frac{\sin⁡(\theta_a)}{L_a} +\frac{\sin⁡(\theta_b)}{L_b} +\frac{\sin⁡(\theta_c)}{L_c} =0,即三条分割线的有向曲率之和为0。

对于边界是圆的情况,由于整个图形旋转对称,一个最优结果绕图形中心任意旋转还是最优结果。为了这种情况能够让结果唯一,我们可以选择固定一个和边界相交的分割线在便加上交点的位置,不允许移动,于是在方程中,关于这个分割线在这个点必须垂直边界的约束条件也可以删除了。

更进一步,我们也可以选择在限制某些内部交点位置的情况下的极小值(但是不限制交于这个点的三条分割线的曲率),于是要求三条分割线在这点两两夹角相等的约束条件可以去除,但是需要保留这三条分割线的曲率和为零这个约束条件,即sin(θa)La+sin(θb)Lb+sin(θc)Lc=0\frac{\sin⁡(\theta_a )}{L_a} +\frac{\sin⁡(\theta_b )}{L_b} +\frac{\sin⁡(\theta_c )}{L_c} =0。 而三条分割线两两夹角相等的公式形式即方程组

{sa,b,cxsxpLscos(θs)+ysypLssin(θs)=0sa,b,cysypLscos(θs)xsxpLssin(θs)=0\begin{cases}\sum_{s\in {a,b,c}}\frac{x_s-x_p}{L_s}\cos⁡(\theta_s)+\frac{y_s-y_p}{L_s}\sin⁡(\theta_s )=0\\ \sum_{s\in{a,b,c}}\frac{y_s-y_p}{L_s}\cos⁡(\theta_s)-\frac{x_s-x_p}{L_s}\sin⁡(\theta_s)=0 \end{cases}

第二部分 边界约束条件方程分析

对于边界为单位圆的情况,绪论里面给出了边界上点的约束方程为 E1=sin(θ)(1cos(u)x0sin(u)y0)cos(θ)(sin(u)x0cos(u)y0)=0E_1=\sin(θ)(1-\cos⁡(u) x_0-\sin⁡(u) y_0 )-\cos⁡(θ)(\sin⁡(u) x_0-\cos⁡(u) y_0)=0

同样对于边界为一般的圆的情况,假设边界上点为P(a+rcos(u),b+rsin(u))P(a+r\cos(u),b+r\sin(u)), 这条分界线另外一个端点为A(x0,y0)A(x_0,y_0 ), 于是PA方向向量为(x0arcos(u),y0brsin(u))(x_0-a-r \cos⁡(u),y_0-b-r \sin⁡(u) ),这个向量顺时针旋转θ度后变为((x0arcos(u))cos(θ)+(y0brsin(u))sin(θ),(y0brsin(u))cos(θ)(x0arcos(u))sin(θ))((x_0-a-r \cos⁡(u) )\cos⁡(θ)+( y_0-b-r\sin⁡(u) )\sin⁡(θ),(y_0-b-r\sin⁡(u) )\cos⁡(θ)-(x_0-a-r\cos⁡(u) )\sin⁡(θ))和切线方向(sin(u),cos(u))(-\sin⁡(u),\cos(u))垂直,所以有(x0arcos(u))sin(u)cos(θ)+(y0brsin(u))sin(u)sin(θ)(y0brsin(u))cos(u)cos(θ)+(x0arcos(u))cos(u)sin(θ)=0(x_0-a-r\cos⁡(u))\sin⁡(u)\cos⁡(θ)+(y_0-b-r\sin⁡(u))\sin⁡(u)\sin⁡(θ)-(y_0-b-r \sin⁡(u))\cos⁡(u)\cos⁡(θ)+(x_0-a-r\cos⁡(u))\cos⁡(u)\sin⁡(θ)=0, 即Ec=cos(θ)((x0a)sin(u)(y0b)cos(u))sin(θ)((y0b)sin(u)+(x0a)cos(u)r)=0E_c=\cos⁡(θ) ((x_0-a)\sin⁡(u)-(y_0-b)\cos⁡(u))-\sin⁡(θ) ((y_0-b)\sin⁡(u)+(x_0-a) \cos⁡(u)-r)=0

straightmin

而对于边界是直线的情况,结果会更简单一些,比如对于下图边界为x=1的情况

tan(θ)+ypy01x0=0\tan⁡(θ)+\frac{y_p-y_0}{1-x_0}=0 , 即 E2=sin(θ)(1x0)+cos(θ)(ypy0)=0E_2=\sin⁡(θ)(1-x_0 )+\cos⁡(θ) (y_p-y_0 )=0 类似,对于边界x=0,y=1,y=0,我们分别有约束公式{E3=x0sin(θ)(ypy0)cos(θ)=0E4=(1y0)sin(θ)(xpx0)cos(θ)=0E5=sin(θ)y0(x0xp)cos(θ)=0\begin{cases}E_3=x_0\sin⁡(θ)-(y_p-y_0 )\cos⁡(θ)=0\\E_4=(1-y_0 )\sin⁡(θ)-(x_p-x_0 )\cos⁡(θ)=0\\ E_5=\sin⁡(θ) y_0-(x_0-x_p )\cos⁡(θ)=0\end{cases}

为了统一起见,对于边界点的x坐标或y坐标,我们后面有时也会统一用参数u表示,而约束方程统一用形式函数E(θ,u,x0,y0)=0E(θ,u,x_0,y_0 )=0表示。

对于一般情况,设直线边界上点为P(au+b,cu+d)P(au+b,cu+d),于是PA方向向量为(x0aub,y0cud)(x_0-au-b,y_0-cu-d),PA顺时针旋转θ度后变为((x0aub)cos(θ)+(y0cud)sin(θ),(y0cud)cos(θ)(x0aub)sin(θ))((x_0-au-b)\cos⁡(θ)+(y_0-cu-d)\sin⁡(θ), (y_0-cu-d)\cos⁡(θ)-(x_0-au-b)\sin⁡(θ))直线切线方向为(a,c),所以(a(x0aub)cos(θ)+a(y0cud)sin(θ)+c(y0cud)cos(θ)c(x0aub)sin(θ))=0(a(x_0-au-b)\cos⁡(θ)+a(y_0-cu-d)\sin⁡(θ)+c(y_0-cu-d)\cos⁡(θ)-c(x_0-au-b) \sin⁡(θ) )=0,即El=cos(θ)(ax0+cy0abcd(a2+c2)u)+sin(θ)(ay0cx0+bcad)=0E_l=\cos⁡(θ) (ax_0+cy_0-ab-cd-(a^2+c^2)u)+\sin⁡(θ) (ay_0-cx_0+bc-ad)=0.

第三部分 面积约束条件

另外我们还需要提供每个区域的面积计算公式,如下图部分面积公式可以写成

areamin

Sa,b,c,d=12a,b(xaybxbyaSDa,b(2θasin(2θa))4sin2(θa)=SaveS_{a,b,c,d}=\frac12\sum_{a,b}(x_ay_b-x_by_a-\frac{SD_{a,b} (2\theta_a-\sin⁡(2\theta_a ))}{4\sin^2⁡(\theta_a)}=S_{ave} (IV), 其中SDa,b=(xaxb)2+(yayb)2SD_{a,b}=(x_a-x_b )^2+(y_a-y_b )^2代表两点间距离平方。

但是某个区域的某条边还有可能不是田埂而是田地的边界,这时其两个端点都在边界上。如果边界这时是单位圆,设两个端点分别为(cos(u),sin(u)),(cos(v),sin(v))(\cos(u), \sin(u)), (\cos(v), \sin(v)), 那么我们相当于可以上面公式中的θaθ_a使用vu2\frac{v-u}2替换即可。

由此,对于一个n个区域,其中有m个内部交叉点(三条田埂相交)和k个边界交叉点(田埂和边界交叉)的配置,使用欧拉定理可以计算出m+k=2(n-1)。而田埂的数目e满足2e=3m+k且e=m+n-1。

对于m个内部交叉点,每个点可以有横坐标,纵坐标两个变量,共2m个变量,而k各边界交叉点,根据落在边界上的位置,每个边界点需要一个变量,共k个变量。另外每条田埂有一个有向弧度参数θ,所以e条田埂共提供e个变量。由此我们总共定义了2m+k+e个变量。

另外一方面,每个内部交叉点提供了3条(I)中的方程,每个边界交叉点提供了一条类似(II),(III)的方程,共3m+k条方程,另外n个区域的面积约束提供了n-1条独立的方程(IV) (n个区域面积之和为已知常数,所以要去掉一个面积约束,总共只有n-1个面积约束)。所以我们总共有3m+k+n-1=2m+k+e条方程。方程和变量数目相同,理论上已经可以用这些方程把这些变量全部计算出来。

但是由于方程为超越方程,很难有系统的方法找出所有解,我们可以尝试使用牛顿迭代法从某些给定的初始近似解求出符合条件的数值解。 为此我们需要进一步各个表达式的各种一阶偏导数,其中在点或边的形式不同时,由于包含的参数数目不同,偏导数形式可能还会略有不同。

比如:

Lsxs=xsxpLs,Lsys=ysypLs\frac{\partial L_s}{\partial x_s}=\frac{x_s-x_p}{L_s} ,\frac{\partial L_s}{\partial y_s}=\frac{y_s-y_p}{L_s} , Lsxp=xsxpLs=Lsxs,Lsyp=ysypLs=Lsys\frac{∂L_s}{∂x_p}=-\frac{x_s-x_p}{L_s} =-\frac{∂L_s}{∂x_s},\frac{∂L_s}{∂y_p}=-\frac{y_s-y_p}{L_s} =-\frac{∂L_s}{∂y_s} Ds,xxs=cos(θsLsDs,x(xsxp)Ls2,Ds,xys=sin(θs)LsDs,x(ysyp)Ls2\frac{∂D_{s,x}}{∂x_s}=\frac{\cos⁡(θ_s}{L_s} -\frac{D_{s,x}(x_s-x_p)}{L_s^2},\frac{∂D_{s,x}}{∂y_s}=\frac{\sin⁡(θ_s)}{L_s} -\frac{D_{s,x} (y_s-y_p )}{L_s^2}, Ds,xxp=Ds,xxs,Ds,xyp=Ds,xys\frac{∂D_{s,x}}{∂x_p}=-\frac{∂D_{s,x}}{∂x_s},\frac{∂D_{s,x}}{∂y_p}=-\frac{∂D_{s,x}}{∂y_s} Ds,yxs=sin(θsLsDs,y(xsxp)Ls2,Ds,yys=cos(θs)LsDs,y(ysyp)Ls2\frac{∂D_{s,y}}{∂x_s}=-\frac{\sin⁡(θ_s}{L_s} -\frac{D_{s,y} (x_s-x_p )}{L_s^2},、\frac{∂D_{s,y}}{∂y_s}=\frac{\cos⁡(θ_s)}{L_s} -\frac{D_{s,y} (y_s-y_p)}{L_s^2}, Ds,yxp=Ds,yxs,Ds,yyp=Ds,yys\frac{∂D_{s,y}}{∂x_p}=-\frac{∂D_{s,y}}{∂x_s},\frac{∂D_{s,y}}{∂y_p}=-\frac{∂D_{s,y}}{∂y_s } Ds,xθs=Ds,y,Ds,yθs=Ds,x\frac{∂D_{s,x}}{∂θ_s}=D_{s,y},\frac{∂D_{s,y}}{∂θ_s}=-D_{s,x}

{F1θs=cos(θs)LsF1xs=sin(θs)(xsxp)Ls3,F1ys=sin(θs)(ysyp)Ls3F1xp=sa,b,cF1xs,F1yp=sa,b,cF1ys\begin{cases}\frac{∂F_1}{∂θ_s}=\frac{cos⁡(θ_s)}{L_s}\\ \frac{∂F_1}{∂x_s}=-\frac{sin⁡(θ_s ) (x_s-x_p )}{L_s^3},\frac{∂F_1}{∂y_s}=-\frac{\sin⁡(θ_s ) (y_s-y_p )}{L_s^3}\\ \frac{∂F_1}{∂x_p}=-\sum_{s\in{a,b,c}}\frac{∂F_1}{∂x_s},\frac{∂F_1}{∂y_p}=-\sum_{s\in{a,b,c}}\frac{∂F_1}{∂y_s}\end{cases}

F2θs=Ds,xθs=Ds,y,F3θs=Ds,yθs=Ds,x\frac{∂F_2}{∂θ_s}=\frac{∂D_{s,x}}{∂θ_s}=D_{s,y},\frac{∂F_3}{∂θ_s}=\frac{∂D_{s,y}}{∂θ_s}=-D_{s,x} F2xs=Ds,xxs=cos(θs)LsDs,x(xsxp)Ls2\frac{∂F_2}{∂x_s}=\frac{∂D_{s,x}}{∂x_s}=\frac{\cos⁡(θ_s )}{L_s} -\frac{D_{s,x} (x_s-x_p )}{L_s^2},F2ys=Ds,xys=sin(θs)LsDs,x(ysyp)Ls2\frac{∂F_2}{∂y_s}=\frac{∂D_{s,x}}{∂y_s}=\frac{\sin⁡(θ_s )}{L_s} -\frac{D_{s,x} (y_s-y_p )}{L_s^2} F3xs=Ds,yxs=sin(θs)LsDs,y(xsxp)Ls2\frac{∂F_3}{∂x_s}=\frac{∂D_{s,y}}{∂x_s}=-\frac{\sin⁡(θ_s )}{L_s} -\frac{D_{s,y} (x_s-x_p )}{L_s^2},F3ys=Ds,yys=cos(θs)LsDs,y(ysyp)Ls2\frac{∂F_3}{∂y_s}=\frac{∂D_{s,y}}{∂y_s}=\frac{\cos⁡(θ_s )}{L_s} -\frac{D_{s,y} (y_s-y_p )}{L_s^2} F2xp=sa,b,cDs,xxp=sa,b,cF2xs\frac{∂F_2}{∂x_p}=\sum_{s\in{a,b,c}}\frac{∂D_{s,x}}{∂x_p}=-\sum_{s∈{a,b,c}}\frac{∂F_2}{∂x_s}, F2yp=sa,b,cDs,xyp=sa,b,cF2ys\frac{∂F_2}{∂y_p}=\sum_{s\in{a,b,c}}\frac{∂D_{s,x}}{∂y_p}=-\sum_{s∈{a,b,c}}\frac{∂F_2}{∂y_s}, F3xp=sa,b,cDs,yxp=sa,b,cF3xs\frac{∂F_3}{∂x_p}=\sum_{s\in{a,b,c}}\frac{∂D_{s,y}}{∂x_p}=-\sum_{s∈{a,b,c}}\frac{∂F_3}{∂x_s}, F3yp=sa,b,cDs,yyp=sa,b,cF3ys\frac{∂F_3}{∂y_p}=\sum_{s\in{a,b,c}}\frac{∂D_{s,y}}{∂y_p}=-\sum_{s∈{a,b,c}}\frac{∂F_3}{∂y_s},

另外我们还要处理各种边界条件和面积约束的一阶偏导数。

{E1θ=cos(θ)(1cos(u)x0sin(u)y0)+sin(θ)(sin(u)x0cos(u)y0)E1u=sin(θ)(sin(u)x0cos(u)y0)cos(θ)(cos(u)x0+sin(u)y0)E1x0=sin(θ)cos(u)cos(θ)sin(u)=sin(θ+u)E1y0=sin(θ)sin(u)+cos(θ)cos(u)=cos(θ+u)E2θ=cos(θ)(1x0)sin(θ)(ypy0)E2yp=cos(θ)E2x0=sin(θ)E2y0=cos(θ)E3θ=cos(θ)x0+sin(θ)(ypy0)E3yp=cos(θ)E3x0=sin(θ)E3y0=cos(θ)E4θ=cos(θ)(1y0)+sin(θ)(xpx0)E4xp=cos(θ)E4x0=cos(θ)E4y0=sin(θ)E5θ=cos(θ)y0sin(θ)(xpx0)E5xp=cos(θ)E5x0=cos(θ)E5y0=sin(θ)\begin{cases} \frac{∂E_1}{∂θ}= \cos(θ)(1-\cos⁡(u) x_0-\sin⁡(u) y_0 )+\sin(θ)(\sin⁡(u) x_0-\cos⁡(u) y_0)\\ \frac{∂E_1}{∂u}= \sin(θ)(\sin⁡(u) x_0-\cos⁡(u) y_0 )-\cos⁡(θ)(\cos⁡(u) x_0+\sin⁡(u) y_0)\\ \frac{∂E_1}{∂x_0}= -\sin(θ)\cos⁡(u)-\cos⁡(θ)\sin⁡(u)=-\sin⁡(θ+u)\\ \frac{∂E_1}{∂y_0}=-\sin(θ)\sin⁡(u)+\cos⁡(θ)\cos⁡(u)=\cos⁡(θ+u)\\ \frac{∂E_2}{∂θ}=\cos⁡(θ)(1-x_0 )-\sin⁡(θ) (y_p-y_0 )\\ \frac{∂E_2}{∂y_p}=\cos⁡(θ)\\ \frac{∂E_2}{∂x_0}=-\sin⁡(θ)\\ \frac{∂E_2}{∂y_0}=-\cos⁡(θ)\\ \frac{∂E_3}{∂θ}=\cos⁡(θ)x_0+\sin⁡(θ) (y_p-y_0 )\\ \frac{∂E_3}{∂y_p}=-\cos⁡(θ)\\ \frac{∂E_3}{∂x_0}=\sin⁡(θ)\\ \frac{∂E_3}{∂y_0}=\cos⁡(θ)\\ \frac{∂E_4}{∂θ}=\cos⁡(θ)(1-y_0)+\sin⁡(θ)(x_p-x_0 )\\ \frac{∂E_4}{∂x_p}=-\cos⁡(θ)\\ \frac{∂E_4}{∂x_0}=\cos⁡(θ)\\ \frac{∂E_4}{∂y_0}=-\sin(θ)\\ \frac{∂E_5}{∂θ}=\cos⁡(θ)y_0-\sin⁡(θ) (x_p-x_0 )\\ \frac{∂E_5}{∂x_p}=\cos⁡(θ)\\ \frac{∂E_5}{∂x_0}=-\cos⁡(θ)\\ \frac{∂E_5}{∂y_0}=\sin(θ)\end{cases}

关于区域面积的偏导数,有

{SDa,bxa=2(xaxb)SDa,bxb=2(xbxa)SDa,bya=2(yayb)SD(a,b)yb=2(ybya)\begin{cases}\frac{∂SD_{a,b}}{∂x_a} =2(x_a-x_b )\\ \frac{∂SD_{a,b}}{∂x_b} =2(x_b-x_a )\\ \frac{∂SD_{a,b}}{∂y_a} =2(y_a-y_b )\\ \frac{∂SD_(a,b)}{∂y_b} =2(y_b-y_a )\end{cases}

得出面积偏导数(区域顶点按顺序排列,d代表a前面一个点,b代表a后面一个点)

{Sθa=SDa,b(1θactg(θa))2sin2(θa)Sxa=12(ybyd(2θasin(2θa))xaxb2sin2(θa)(2θdsin(2θd))xaxd2sin2(θd))Sya=12(xdxb(2θasin(2θa))yayb2sin2(θa)(2θdsin(2θd))yayd2sin2(θd))\begin{cases}\frac{∂S}{∂θ_a} =-\frac{SD_{a,b} (1-θ_a ctg(θ_a ))}{2 \sin^2⁡(θ_a )}\\ \frac{∂S}{∂x_a} =\frac12 (y_b-y_d-(2θ_a-\sin⁡(2θ_a ) )\frac{x_a-x_b }{2\sin^2⁡(θ_a)}-(2θ_d-\sin⁡(2θ_d ) )\frac{x_a-x_d }{2\sin^2⁡(θ_d )})\\ \frac{∂S}{∂y_a} =\frac12 (x_d-x_b-(2θ_a-\sin⁡(2θ_a ))\frac{y_a-y_b}{2\sin^2⁡(θ_a )}-(2θ_d-\sin⁡(2θ_d ))\frac{y_a-y_d}{2\sin^2⁡(θ_d)})\end{cases}

下面我们对m个内点进行编号为N1,N2,,NmN_1,N_2,\dots,N_m;对k个边界点标号为N(m+1),N(m+2),,N(m+k)N_(m+1),N_(m+2),\dots,N_(m+k)

e条田埂编号为E1,E2,,EmE_1,E_2,\dots,E_m, 而对于这m条田埂并且任意确定一个方向,然后使用E1,E2,,EmE_{-1},E_{-2},\dots,E_{-m}分别代表交换这些田埂的起点和终点后得到的有向田埂。

另外还需要n个区域编号为C1,C2,,CnC_1,C_2,\dots,C_n

对于每个内点NiN_i函数aia_i=get_a(Ni)(N_i), bib_i=get_b(Ni)(N_i), cic_i=get_c(Ni)(N_i), 返回三条以内点NiN_i为起点的有向田埂,对于每个边界点NiN_i,函数get_a(Ni)(N_i)返回以NiN_i为起点的田埂,函数get_b和get_c报错。对于每条田埂EiE_i,有函数get_start(Ei)(E_i)和get_end(Ei)(E_i)返回田埂的两个端点。对于每个区域CiC_i,函数get_degree(Ci)(C_i)返回区域的顶点数目,v[i,h]=get_vertex(CiC_i,h)返回区域CiC_i的第h个顶点,其中所有顶点按逆时针顺序依次排列。另外如果从NiN_iNjN_j存在田埂,那么get_edge(Ni,Nj)(N_i,N_j)返回这条田埂,不然报错。

另外函数get_x(Ni)(N_i), get_y(Ni)(N_i) 依次返回NiN_i的横坐标和纵坐标, get_theta(Ei)(E_i)返回EiE_i的角度θiθ_i。而显然有get_theta(Ei)(E_{-i})=- get_theta(Ei)(E_i)=- θiθ_i。另外我们可以使用符号θi,aθ_{i,a}代表get_theta(get_edge(Ni(N_i, get_a(Ni)(N_i)))。

于是在上面约定下,我们有2m+k+e个变量x1,y1,x2,y2,,xm,ym,u1,u2,,uk,θ1,θ2,,θex_1,y_1,x_2,y_2,\dots,x_m,y_m,u_1,u_2,\dots,u_k,θ_1,θ_2,\dots,θ_e

另外对于每个内点N_i,我们有方程

{(F1(Ni,ai,bi,ci,θ(i,a),θ(i,b),θ(i,c))=0F2(Ni,ai,bi,ci,θ(i,a),θ(i,b),θ(i,c))=0F3(Ni,ai,bi,ci,θ(i,a),θ(i,b),θ(i,c))=0\begin{cases}(F_1 (N_i,a_i,b_i,c_i,θ_(i,a),θ_(i,b),θ_(i,c) )=0\\ F_2 (N_i,a_i,b_i,c_i,θ_(i,a),θ_(i,b),θ_(i,c) )=0\\ F_3 (N_i,a_i,b_i,c_i,θ_(i,a),θ_(i,b),θ_(i,c))=0\end{cases}, 简记为{(F1(i)=0F2(i)=0F3(i)=0\begin{cases}(F_1^{(i)}=0\\F_2^{(i)}=0\\F_3^{(i)}=0\end{cases}