问题

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

<script id="MathJax-Element-41" type="math/tex; mode=display">\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.</script>
与定子曲线方程共轭的转子曲线方程：
<script id="MathJax-Element-42" type="math/tex; mode=display">\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.</script>

当中：
1. R=48.78<script id="MathJax-Element-43" type="math/tex">R=48.78</script> 为导圆半径，
2. r=8.13<script id="MathJax-Element-44" type="math/tex">r=8.13</script> 为滚圆半径。
3. z2=z1?1<script id="MathJax-Element-45" type="math/tex">z_2=z_1-1</script> 为转子头数。
4. e=7.05<script id="MathJax-Element-46" type="math/tex">e=7.05</script> 为偏心距，
5. f=er<script id="MathJax-Element-47" type="math/tex">f=\dfrac{e}{r}</script>,
6 θ=tan?1sin(z1τ)f+cos(z1τ)?τ<script id="MathJax-Element-48" type="math/tex">\theta =\tan^{-1}\dfrac{\sin\left(z_1\tau\right)}{f+\cos\left(z_1\tau\right)}-\tau</script>,
7. re=12.6<script id="MathJax-Element-49" type="math/tex">r_e=12.6</script>为等距半径，
8. φ=sin?1[fsin(θ+τ)]<script id="MathJax-Element-50" type="math/tex">\varphi=\sin^{-1}[f\sin(\theta+\tau)]</script>,
9. ψφ=z1z1?1<script id="MathJax-Element-51" type="math/tex">\dfrac{\psi}{\varphi}=\dfrac{z_1}{z_1-1}</script>

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

一些猜測：最初问这个问题的同学，预计是某高校不久之后自尽的一位研究生。

ta兴许另一些问题，正好我自己也碰到了一些麻烦，没有及时解答。偶然从新闻上自尽的那位同学的专业和大致时间，与该同学提问之后不再上线的时间以及专业的匹配作出上面猜測的。

该同学后来的帖子还许诺酬劳求解兴许的问题（也就是知道面积反求其他參数的问题，实际用牛顿法之类的局部非线性优化方法就能easy求解）。当初也没有太在意，没想到。已经是压力大到这样的程度。唏嘘感慨。就算如今解出兴许的问题又有什么意义呢？（Mar 2016）

原理

我把 z2=z1?1<script id="MathJax-Element-12" type="math/tex">z_2=z_1-1</script> 误写作 z2=1?z1<script id="MathJax-Element-13" type="math/tex">z_2=1-z_1</script> 结果大大添加了问题的难度。

原始问题的两条曲线简单得差点儿不值得做。

这是该夸我呢还是？？？？

这样的情况，曲线太复杂了，符号计算的可能性太小，投入和收获之间不成比例。

所以。近似的数值计算是合理的。

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

第一个问题。由于近似计算。用曲线 ?x(τ),y(τ)?<script id="MathJax-Element-14" type="math/tex">\langle x(\tau),y(\tau)\rangle</script> 当自己相交时相应的 (x,y)<script id="MathJax-Element-15" type="math/tex">(x,y)</script> 坐标相等来近似计算的。即 ?x(τ1)?x(τ2),y(τ1)?y(τ2)?=?0,0?<script id="MathJax-Element-16" type="math/tex">\langle x(\tau_1)-x(\tau_2),y(\tau_1)-y(\tau_2)\rangle=\langle0,0\rangle</script>, 这是一个二元非线性方程组求根的问题。能够用牛顿法来解决。 用图示或其他方法确定出相交时曲线两部分的τ1<script id="MathJax-Element-17" type="math/tex">\tau_1</script> 和 τ2<script id="MathJax-Element-18" type="math/tex">\tau_2</script> 的大致范围能够帮助确定初值。由于对称性。仅仅须计算十分之中的一个就可以。

第二个问题，前面一篇博客刚刚提到，把要计算部分的轮廓离散化成多边形，然后用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