求转子曲线所包围的封闭区域的面积

问题

碰到这样的问题，感觉很神奇。
定子方程，短幅内摆线方程：

{x1y1=(R−r)sinτ+esin(z2τ)−resinθ=(R−r)cosτ−ecos(z2τ)+recosθ

与定子曲线方程共轭的转子曲线方程：

{x2=x1cos(φ−ψ)−y1sin(φ−ψ)−esin(ψ)y2=x1sin(φ−ψ)+y1cos(φ−ψ)−ecos(ψ)

其中：
1. R=48.78 为导圆半径，
2. r=8.13 为滚圆半径，
3. z2=z1−1 为转子头数。
4. e=7.05 为偏心距，
5. f=er,
6 θ=tan−1sin(z1τ)f+cos(z1τ)−τ,
7. re=12.6为等距半径，
8. φ=sin−1[fsin(θ+τ)],
9. ψφ=z1z1−1

求如下转子曲线外轮廓所包围的面积。

一些推测：最初问这个问题的同学，估计是某高校不久之后自尽的一位研究生。ta后续还有一些问题，正好我自己也碰到了一些麻烦，没有及时解答。偶然从新闻上自尽的那位同学的专业和大致时间，与该同学提问之后不再上线的时间以及专业的匹配作出上面推测的。该同学后来的帖子还许诺报酬求解后续的问题（也就是知道面积反求其它参数的问题，实际用牛顿法之类的局部非线性优化方法就能容易求解）。当初也没有太在意，没想到，已经是压力大到这种程度。唏嘘感慨。就算现在解出后续的问题又有什么意义呢？（Mar 2016）

原理

我把 z2=z1−1 误写作 z2=1−z1 结果大大增加了问题的难度。原始问题的两条曲线简单得几乎不值得做。这是该夸我呢还是？？？？

这种情况，曲线太复杂了，符号计算的可能性太小，投入和收获之间不成比例。所以，近似的数值计算是合理的。

涉及两方面的主要问题：
1. 如何找出外轮廓；
2. 如何计算外轮廓包围的面积。

第一个问题，因为近似计算，用曲线 ⟨x(τ),y(τ)⟩ 当自己相交时对应的 (x,y) 坐标相等来近似计算的。即 ⟨x(τ1)−x(τ2),y(τ1)−y(τ2)⟩=⟨0,0⟩, 这是一个二元非线性方程组求根的问题，可以用牛顿法来解决。 用图示或其它方法确定出相交时曲线两部分的τ1 和 τ2 的大致范围可以帮助确定初值。由于对称性，只须计算十分之一即可。

第二个问题，前面一篇博客刚刚提到，把要计算部分的轮廓离散化成多边形，然后用Shoelace公式计算多边形的面积即可。对应于代码中的 PolygonSignedArea 函数。

答案

直接上代码：

ClearAll["Global`*"];R = 48.78;r = 8.13;z1 = R/r;z2 = 1 - z1;e = 7.05;f = r/e;re = 12.6;θ = ArcTan[Sin[z1 τ]/(f + Cos[z1 τ])] - τ;φ = ArcSin[f Sin[θ + τ]] - θ;ψ = z1/(z1 - 1) φ;curve01 = {(R - r) Sin[τ] + e Sin[z2 τ] -re Sin[θ], (R - r) Cos[τ] - e Cos[z2 τ] +re Cos[θ]} // FullSimplify;curve02 = {curve01[[1]] Cos[φ - ψ] - curve01[[2]] Sin[φ - ψ] - e Sin[ψ], curve01[[1]] Sin[φ - ψ] +  curve01[[2]] Cos[φ - ψ] - e Cos[ψ]} //  Simplify;base = ParametricPlot[Evaluate[curve02], {\[Tau], 0, 5 \[Pi]},Exclusions -> None, MaxRecursion -> 15, PlotPoints -> 1000];r1 = FindRoot[ (curve02 /. τ -> x) == (curve02 /. τ -> y), {x, .5}, {y, 5.5}];p1 = y /. FindRoot[ (ArcTan @@ (curve02 /. τ -> y)) == 3 Pi/ 10 , {y, .55}];top = x /. FindRoot[ curve02[[1]] /.  τ -> x , {x, 5}];arc = Join[Table[ curve02 , {τ, top, y /. r1 , .0001}],    Table[ curve02 , {τ, x /. r1, p1 , .0001}]];Show[base,Graphics[{{Red,  Line[{curve02 /.τ -> top, {0, 0}, curve02 /.τ -> p1}], Line@arc }}]];PolygonSignedArea[pts_?MatrixQ] := Total[Det /@ Partition[pts, 2, 1, 1]]/2;area = 10 PolygonSignedArea[Reverse@Join[{{0, 0}}, arc]];Print[Style["Area is: ", Blue, 20],  Style[NumberForm[area, 10], Red, 23]];

Area is: 7936.859683

上面是求多边形面积近似曲线包围区域的面积，基于Green定理的近似。这种方法误差跟曲线弧上离散点步长有关。另外一种直接用Green定理的方法，接上面的代码再计算得到：

5 (NIntegrate[-First[curve02] D[Last@curve02,τ]+Last[curve02] D[First@curve02, τ],{τ, top,y /. r1}] +NIntegrate[-First[curve02] D[Last@curve02,τ]+Last[curve02] D[First@curve02,τ], {τ,x/.r1, p1}]) // NumberForm[#, 15] &

结果是：

7945.519351112369