楕円の周の長さの近似式

テイラー展開による近似式およびパデ近似による近似式とラマヌジャンの近似式の紹介と評価。追記して,関孝和による近似式も。

楕円の周の長さ

$\displaystyle \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$ で表される楕円の周の長さ。

便宜上,$a > b > 0$ として $a$ を長半径,$b$ を短半径と呼ぶ。離心率 $e$ を使って書くと,

\begin{eqnarray}
b &=& a \sqrt{1-e^2} \\
\therefore\ \ e^2 &=& \frac{a^2 – b^2}{a^2}
\end{eqnarray}

このとき,楕円の周の長さ $L$ は,

$$ E(k) \equiv \int_0^{\frac{\pi}{2}} \sqrt{1 -k^2 \sin^2\theta} \,d\theta $$

で定義される第二種完全楕円積分 $E(k)$ を使って

$$L = 4 a E(e) = 4 a \int_0^{\frac{\pi}{2}} \sqrt{1 -e^2 \sin^2\theta} \,d\theta$$

と書けるのであった。以下のページを参照:

楕円の周長が「第二種」楕円積分なら,「第一種」楕円積分はどこにあらわれるか,という話も上記のページのここに。

Maxima の第二種完全楕円積分 elliptic_ec (m)

Maxima には第二種完全楕円積分elliptic_ec (m) として組み込まれている。

Function: elliptic_ec (<m>)
  The complete elliptic integral of the second kind, defined as
  integrate(sqrt(1 - m*sin(x)^2), x, 0, %pi/2)

elliptic_ec (m) の引数 $m$ と離心率 $e$ の関係に注意して,ここで使う第二種完全楕円積分 $E(e)$ を次のように定義する。

In [1]:
E(e):= elliptic_ec (e**2);
Out[1]:
\[\tag{${\it \%o}_{1}$}E\left(e\right):={\it elliptic\_ec}\left(e^2\right)\]

離心率 $e$ は $0 \le e < 1$ であるから,この範囲で $E(e)$ のグラフを描いてみる。(本サイトでは Maxima のグラフには draw2d() などの draw 系をお薦めしているのであるが,手っ取り早くグラフを描きたいときは plot2d() で。)

In [2]:
plot2d(E(e), [e, 0, 1], [xlabel, "e"], [ylabel, "E(e)"])$

$E(e)$ は単調減少関数のようである。最大値は $e=0$ (真円)のときで,

In [3]:
E(0);
Out[3]:
\[\tag{${\it \%o}_{4}$}\frac{\pi}{2}\]

このとき,周の長さ $L$ は半径 $a$ の円周となり,

In [4]:
L = 4*a*E(0);
Out[4]:
\[\tag{${\it \%o}_{5}$}L=2\,\pi\,a\]

最小値は $e \rightarrow 1$ のときで,

In [5]:
E(1);
Out[5]:
\[\tag{${\it \%o}_{6}$}1\]

テイラー展開による楕円周の近似式

\begin{eqnarray}
L = 4 a E(e)
&=& 4 a \int_0^{\frac{\pi}{2}} f(e)\, dx\\
f(e) &\equiv& \sqrt{1 -e^2 \sin^2 x}
\end{eqnarray}

In [6]:
f(e):= sqrt(1 - e**2*sin(x)^2);
Out[6]:
\[\tag{${\it \%o}_{7}$}f\left(e\right):=\sqrt{1-e^2\,\sin ^2x}\]

被積分関数 $f(e)$ を $e = 0$ のまわりでテイラー展開してから積分する。まずは $e^8$ くらいまで展開して $e^4$ までの展開と比べてみると…

In [7]:
ft4: trunc(taylor(f(e), e, 0, 4));
ft8: trunc(taylor(f(e), e, 0, 8));
Out[7]:
\[\tag{${\it \%o}_{8}$}1-\frac{e^2\,\sin ^2x}{2}-\frac{e^4\,\sin ^4x}{8}+\cdots \]
Out[7]:
\[\tag{${\it \%o}_{9}$}1-\frac{e^2\,\sin ^2x}{2}-\frac{e^4\,\sin ^4x}{8}-\frac{e^6\,\sin ^6x}{16}-\frac{5\,e^8\,\sin ^8x}{128}+\cdots \]

これを $x=0$ から $\pi/2$ まで定積分すると…

In [8]:
intft4: integrate(ft4, x, 0, %pi/2), expand;
intft8: integrate(ft8, x, 0, %pi/2), expand;
Out[8]:
\[\tag{${\it \%o}_{10}$}-\frac{3\,\pi\,e^4}{128}-\frac{\pi\,e^2}{8}+\frac{\pi}{2}\]
Out[8]:
\[\tag{${\it \%o}_{11}$}-\frac{175\,\pi\,e^8}{32768}-\frac{5\,\pi\,e^6}{512}-\frac{3\,\pi\,e^4}{128}-\frac{\pi\,e^2}{8}+\frac{\pi}{2}\]
In [9]:
plot2d([E(e), intft4, intft8], [e, 0, 1], [xlabel, "e"], [ylabel, ""], 
       [legend, "E(e)", "taylor 4", "taylor 8"])$

$e < 0.6$ まではグラフはどちらも重なるくらい近似の精度が良さそうであり,以下の計算から,$e=0.9$ での相対誤差がテイラー展開の4次までの計算では $3\%$ 未満,8次までの展開では $1\%$ 未満であることがわかる。

In [10]:
ev(float((intft4-E(e))/E(e)), e=0.9);
ev(float((intft8-E(e))/E(e)), e=0.9);
Out[10]:
\[\tag{${\it \%o}_{14}$}0.02791136804769806\]
Out[10]:
\[\tag{${\it \%o}_{15}$}0.007832160862980865\]

精度としてはこれくらいでいいとすると,楕円の周の長さは以下を使って…

In [11]:
trunc(expand(intft8/(%pi/2)));
Out[11]:
\[\tag{${\it \%o}_{16}$}1-\frac{e^2}{4}-\frac{3\,e^4}{64}-\frac{5\,e^6}{256}-\frac{175\,e^8}{16384}+\cdots \]
$$L_{\rm tay} \simeq 2 \pi a \left(1-\frac{e^2}{4}-\frac{3\,e^4}{64}-\frac{5\,e^6}{256}-\frac{175\,e^8}{16384}+\cdots\right)$$$$e^2 = \frac{a^2 – b^2}{a^2}$$

ラマヌジャンによる楕円周の近似式

ネットで検索すると(要出典)ラマヌジャンによる楕円周長の近似式として以下の式が使われているようだ。

\begin{eqnarray}
L_{\rm Ram} &\equiv& \pi \left(3 (a+b) – \sqrt{(a+3b)(3a+b)} \right) \\
&=& 2 \pi a \left(\frac{3}{2} \left(1+\sqrt{1-e^2}\right)
– \frac{1}{2} \sqrt{\left(1+3\sqrt{1-e^2}\right) \left(3 + \sqrt{1-e^2}\right)}\right)
\end{eqnarray}

近似の精度を比較するために,以下のように規格化された楕円周を定義する。

\begin{eqnarray}
\ell &\equiv& \frac{L}{2\pi a} = \frac{4 a E(e)}{2 \pi a} = \frac{2 E(e)}{\pi} \\
\ell_t &\equiv& \frac{L_{\rm tay}}{2\pi a} = \left(1-\frac{e^2}{4}-\frac{3\,e^4}{64}-\frac{5\,e^6}{256}-\frac{175\,e^8}{16384}\right) \\
\ell_r &\equiv& \frac{L_{\rm Ram}}{2\pi a} =
\frac{3}{2} \left(1+\sqrt{1-e^2}\right)
– \frac{1}{2} \sqrt{\left(1+3\sqrt{1-e^2}\right) \left(3 + \sqrt{1-e^2}\right)}
\end{eqnarray}

In [12]:
ell(e):= 2*E(e)/float(%pi);
ellt(e):= 1-e**2/4 - 3*e**4/64 - 5*e**6/256 - 175*e**8/16384;
ellr(e):= 3/2*(1+sqrt(1-e**2)) - 1/2*sqrt((1+3*sqrt(1-e**2))*(3+sqrt(1-e**2)));
Out[12]:
\[\tag{${\it \%o}_{17}$}{\it ell}\left(e\right):=\frac{2\,E\left(e\right)}{{\it float}\left(\pi\right)}\]
Out[12]:
\[\tag{${\it \%o}_{18}$}{\it ellt}\left(e\right):=1-\frac{e^2}{4}+\frac{\left(-3\right)\,e^4}{64}+\frac{\left(-5\right)\,e^6}{256}+\frac{\left(-175\right)\,e^8}{16384}\]
Out[12]:
\[\tag{${\it \%o}_{19}$}{\it ellr}\left(e\right):=\frac{3}{2}\,\left(1+\sqrt{1-e^2}\right)-\frac{1}{2}\,\sqrt{\left(1+3\,\sqrt{1-e^2}\right)\,\left(3+\sqrt{1-e^2}\right)}\]
In [13]:
plot2d([ell(e), ellr(e)], [e, 0, 1], [xlabel, "e"], [ylabel, ""], 
       [legend, "exact", "Ramanujan"])$

さすが,ラマヌジャン!グラフをみただけでは2本の曲線は重なってしまって見分けがつかない。チンタラとテイラー展開やっている場合ではないです。参りました。念のため,相対誤差も表示しておきます。

In [14]:
plot2d((ell(e)-ellr(e))/ell(e), [e, 0, 1])$
ev((ell(e)-ellr(e))/ell(e), e=1.0);
ev((ell(e)-ellr(e))/ell(e), e=0.9);

Out[14]:
\[\tag{${\it \%o}_{23}$}0.004155032983318442\]
Out[14]:
\[\tag{${\it \%o}_{24}$}7.642828603946162 \times 10^{-6}\]

ラマヌジャンの近似式は,実に $0 \le e \le 1$ の範囲で最大の相対誤差が $0.4\%$ 程度である。$0 \le e \le 0.9$ の範囲なら,相対誤差は $10^{-5}$ 未満!という素晴らしい近似式となっている。

パデ近似

悔し紛れに,8次までのテイラー展開の近似式をパデ化してみた。Maxima では taylor() で展開して pade() に渡すことで,簡単にパデ近似が求められる。

In [15]:
pade(taylor(intft8/(%pi/2), e, 0, 8), 4, 4);
Out[15]:
\[\tag{${\it \%o}_{26}$}\left[ \frac{453\,e^4-2544\,e^2+2816}{125\,e^4-1840\,e^2+2816} \right] \]
In [16]:
define(ellp(e), %[1]);
Out[16]:
\[\tag{${\it \%o}_{27}$}{\it ellp}\left(e\right):=\frac{453\,e^4-2544\,e^2+2816}{125\,e^4-1840\,e^2+2816}\]
In [17]:
plot2d([ellr(e), ellt(e), ellp(e)], [e, 0, 1], 
       [legend, "Ramanujan", "Taylor 8", "Pade 4, 4"])$

まとめ

テイラー展開,およびそれを使ったパデ近似は,我々普通の人々が近似的表現を求める際に使うテクニックである。テイラー展開は展開の次数を上げていけばよりよい精度が得られるわけだが,闇雲に項数を増やしてもおもしろくないので,ここでは $e^8$ までの例をあげた。

さらに,8次までのテイラー展開の結果を使って,$[4/4]$ 次のパデ近似を求めた。8次までのテイラー展開よりは近似の精度があがっているが,それでも,天才ラマヌジャンの近似式には,特に $0.9 \le e \le 1$ あたりでは,全く及ばない。

まぁ,そこまで卑屈になることもない。我々凡人であっても,$0 \le e \le 0.9$ あたりの範囲までなら,十分に精度のよい近似式をテイラー展開およびパデ近似によって(誰でも機械的に)作ることができるのだ,とポジティブに考えることにしよう。

追記:関孝和による楕円周の近似式

\begin{eqnarray}
L_{\rm Seki} &=& 2 \sqrt{4 (a-b)^2 + \pi^2 a b} \\
&=& 2 \pi a \sqrt{\frac{b}{a} + \frac{4}{\pi^2} \left(1 – \frac{b}{a} \right)^2} \\
&=& 2 \pi a \sqrt{\sqrt{1-e^2} + \frac{4}{\pi^2} \left(1 – \sqrt{1-e^2} \right)^2}
\end{eqnarray}\begin{eqnarray}
\ell_s &\equiv& \frac{L_{\rm Seki}}{2\pi a} =
\sqrt{\sqrt{1-e^2} + \frac{4}{\pi^2} \left(1 – \sqrt{1-e^2} \right)^2}
\end{eqnarray}

In [18]:
ells(e):= sqrt(sqrt(1-e**2) + 4/float(%pi)**2 * (1-sqrt(1-e**2))**2);
Out[18]:
\[\tag{${\it \%o}_{30}$}{\it ells}\left(e\right):=\sqrt{\sqrt{1-e^2}+\frac{4}{{\it float}\left(\pi\right)^2}\,\left(1-\sqrt{1-e^2}\right)^2}\]
In [19]:
plot2d([ell(e), ells(e), ellt(e)], [e, 0, 1], 
       [legend, "exact", "Seki", "Taylor 8"])$

In [20]:
plot2d([(ells(e)-ell(e))/ell(e), (ellt(e)-ell(e))/ell(e)], [e, 0, 1], 
       [legend, "Seki", "Taylor 8"])$
(ells(1.0)-ell(1.0))/ell(1.0);
(ells(0.9)-ell(0.9))/ell(0.9);

Out[20]:
\[\tag{${\it \%o}_{34}$}0.0\]
Out[20]:
\[\tag{${\it \%o}_{35}$}0.007569077619594445\]

ということで,関孝和による楕円周の近似式も $0 \le e \le 0.9$ の範囲では相対誤差が $1\%$ 未満の優秀な近似式になっていることがわかる。