Project Euler 525
Project Euler 525
题目
Rolling Ellipse
An ellipse \(E(a, b)\) is given at its initial position by equation:
\[\frac{x^2}{a^2}+\frac{(y-b)^2}{b^2}=1\]
The ellipse rolls without slipping along the x axis for one complete turn. Interestingly, the length of the curve generated by a focus is independent from the size of the minor axis:
\(F(a,b)=2\pi \max(a,b)\)
This is not true for the curve generated by the ellipse center. Let \(C(a,b)\) be the length of the curve generated by the center of the ellipse as it rolls without slipping for one turn.
You are given \(C(2, 4) \approx 21.38816906\).
Find \(C(1, 4) + C(3, 4)\). Give your answer rounded to \(8\) digits behind the decimal point in the form \(ab.cdefghij\).
解决方案
题目给的椭圆为\(E(a,b):\dfrac{x^2}{a^2}+\dfrac{(y-b)^2}{b^2}=1.\)它的中心在 \((0,b)\)，并且初始时与 \(x\) 轴相切于原点。
为了推导方便，把椭圆整体平移到以中心为原点的坐标系中（平移不改变曲线长度），即考虑标准椭圆\(\dfrac{x^2}{a^2}+\dfrac{y^2}{b^2}=1.\)之后再把“中心轨迹”的长度算出来即可，因为轨迹长度对整体平移不变。
用参数 \(t\)（常见的“偏心角”参数）表示椭圆上的点：\(P(t)=(x(t),y(t))=(a\cos t,b\sin t), 0\le t<2\pi.\)
该点沿椭圆的切向速度大小为\(\lVert P'(t)\rVert=\sqrt{(-a\sin t)^2+(b\cos t)^2}=\sqrt{a^2\sin^2 t+b^2\cos^2 t}.\)记
\[s(t)=\int_0^t \lVert P'(u)\rVert du\Rightarrow s'(t)=\sqrt{a^2\sin^2 t+b^2\cos^2 t}.\]
滚动且无滑动意味着：椭圆在地面（\(x\) 轴）上滚过的水平距离，恰等于椭圆边界被放到地面上的弧长。因此当接触点从 \(0\) 走到 \(t\) 时，地面接触点的水平位移就是 \(s(t)\)。
在接触点处，椭圆与地面相切；把椭圆刚体旋转到使该切线与 \(x\) 轴重合时，接触点到椭圆中心的分解（沿切向、沿法向） 就直接给出了中心相对于接触点的坐标贡献。
切向方向可以取\(T(t)=P'(t)=(-a\sin t,b\cos t),w(t)=\lVert T(t)\rVert=\sqrt{a^2\sin^2 t+b^2\cos^2 t}.\)单位切向量为\(\hat\tau(t)=\dfrac{T(t)}{w(t)}.\)
在以中心为原点的坐标系中，从接触点指向中心的向量是\(\overrightarrow{PC}(t)=-(a\cos t,b\sin t).\)于是它在切向上的投影（带符号）为
\[ u(t)=\overrightarrow{PC}(t)\cdot \hat\tau(t) =-(a\cos t,b\sin t)\cdot\frac{(-a\sin t,b\cos t)}{w(t)} =\frac{(a^2-b^2)\sin t\cos t}{w(t)}. \]
这个表达式天然带符号，因此不需要再手工分象限改正符号。
\(v(t)\)是中心到切线的距离（取正），也就是中心在法向方向上的坐标。因此，椭圆 \(x^2/a^2+y^2/b^2=1\) 在点 \((x_0,y_0)\) 的切线为\(\dfrac{x x_0}{a^2}+\frac{y y_0}{b^2}=1.\)
代入 \((x_0,y_0)=(a\cos t,b\sin t)\) 得到\(\dfrac{\cos t}{a}x+\dfrac{\sin t}{b}y=1.\) 原点到直线 \(A x+B y=1\) 的距离是 \(1/\sqrt{A^2+B^2}\)，所以
\[ v(t)=\frac{1}{\sqrt{\left(\frac{\cos t}{a}\right)^2+\left(\frac{\sin t}{b}\right)^2}} =\frac{1}{\sqrt{\frac{\cos^2 t}{a^2}+\frac{\sin^2 t}{b^2}}}. \]
地面接触点在地面坐标系中的位置是 \((s(t),0)\)。
而中心相对于接触点的位移，在“切线水平、法线竖直”的坐标系下正是 \((u(t),v(t))\)。
因此中心在地面坐标系中的参数方程为\(X(t)=s(t)+u(t), Y(t)=v(t).\)
我们将中心轨迹长度定义为\(\displaystyle{C(a,b)=\int_0^{2\pi}\sqrt{(X'(t))^2+(Y'(t))^2},dt.}\)
由 \(X'(t)=s'(t)+u'(t)\)、\(Y'(t)=v'(t)\) 得\(\displaystyle{C(a,b)=\int_0^{2\pi}\sqrt{\bigl(s'(t)+u'(t)\bigr)^2+\bigl(v'(t)\bigr)^2}dt,}\)
其中\(s'(t)=\sqrt{a^2\sin^2 t+b^2\cos^2 t},u(t)=\dfrac{(a^2-b^2)\sin t\cos t}{\sqrt{a^2\sin^2 t+b^2\cos^2 t}},v(t)=\dfrac{1}{\sqrt{\frac{\cos^2 t}{a^2}+\frac{\sin^2 t}{b^2}}}.\)
由于上述各量只由 \(\sin t,\cos t\) 通过对称的组合构成，且整周滚动由四个对称象限组成，计算上通常取四分之一再乘 \(4\)：
\[ C(a,b)=4\int_0^{\pi/2}\sqrt{\bigl(s'(t)+u'(t)\bigr)^2+\bigl(v'(t)\bigr)^2}dt. \]
代码
1 | (*s'(t):boundary speed of ellipse*) |