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

问题

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

$$\left\{\quad \begin{array}{rl} x_1&=(R-r) \sin \tau+e \sin \left(z_2 \tau\right)-r_e \sin \theta \\ y_1&=(R-r) \cos\tau-e \cos \left(z_2 \tau\right)+r_e \cos \theta \\ \end{array} \right.$$
与定子曲线方程共轭的转子曲线方程：
$$\left\{\quad \begin{array}{c} x_2=x_1 \cos (\varphi -\psi )-y_1 \sin (\varphi -\psi )-e \sin (\psi ) \\ y_2=x_1 \sin (\varphi -\psi )+y_1 \cos (\varphi -\psi )-e \cos (\psi ) \\ \end{array} \right.$$

其中：
1. $R=48.78$ 为导圆半径，
2. $r=8.13$ 为滚圆半径，
3. $z_2=z_1-1$ 为转子头数。
4. $e=7.05$ 为偏心距，
5. $f=\dfrac{e}{r}$,
6 $\theta =\tan^{-1}\dfrac{\sin\left(z_1\tau\right)}{f+\cos\left(z_1\tau\right)}-\tau$,
7. $r_e=12.6$为等距半径，
8. $\varphi=\sin^{-1}[f\sin(\theta+\tau)]$,
9. $\dfrac{\psi}{\varphi}=\dfrac{z_1}{z_1-1}$

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

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

原理

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

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

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

第一个问题，因为近似计算，用曲线 $\langle x(\tau),y(\tau)\rangle$ 当自己相交时对应的 $(x,y)$ 坐标相等来近似计算的。即 $\langle x(\tau_1)-x(\tau_2),y(\tau_1)-y(\tau_2)\rangle=\langle0,0\rangle$, 这是一个二元非线性方程组求根的问题，可以用牛顿法来解决。 用图示或其它方法确定出相交时曲线两部分的$\tau_1$ 和 $\tau_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