7 近似計算

最後に、(15) の数値計算について考えてみる。

6 節の議論により、(15) は 初等的な関数で表現できないことがわかったが、 コンピュータで展開図を書くのだとすれば、 もしそれが初等的な関数で表現できたとしても、 その初等的関数は、$\sin x$ にせよ $\sqrt{x}$ にせよ コンピュータ内部では実際には近似的な計算を行っているわけであるから、 (15) がそれなりの精度で数値計算できるのであれば、 初等的な式で表せる場合と変わらないことになる。

なお、この節では簡単のため 6 節同様に $a=b=p=R$, $q=0$, $h=\hat{h}R$ の場合を考え、 よって、$r$, $\theta$ は (14), (15) より

$$r(t) = R\sqrt{\hat{h}^2+2(1-\cos t)},$$ (17)

$$\theta(t) = \int_0^t\frac{\sqrt{\hat{h}^2+(1-\cos\phi)^2}}{\hat{h}^2+2(1-\cos\phi)} d\phi$$ (18)

$${(0\leq t\leq 2\pi)}$$

であるとする。

さて、(18) の積分の近似計算であるが、 6 節で求めた楕円積分の標準形を用いるのは、 もしその標準形の楕円積分の値を求める関数が コンピュータにあるなら有効であろうが、 少なくとも C の標準ライブラリには用意されてはおらず、 また、楕円積分は広義積分なので単純に考えれば数値積分は難しいので、 むしろ 6 節のような変換は行わず、 (18) のまま広義積分ではない積分を近似計算するのが 得策だろうと思われる。

積分計算というと、通常は台形公式やシンプンソンの公式などが用いられるが、 (18) の場合は、 必ずしもそれを用いるのが適切とは言えない。 それは、(18) が積分ではあるが 実際には $t$ の関数であり、 各 $t$ に対する値がまんべんなく必要となるので、 その値が必要なときに毎回その $t$ に対して $0$ から $t$ までの積分を シンプソンの公式などで近似計算するとなると 計算量がかなり多くなってしまうからである。

(18) の導関数は

$$\theta'(t)=g_0(t) =\frac{\sqrt{\hat{h}^2+(1-\cos t)^2}}{\hat{h}^2+2(1-\cos t)}$$ (19)

と単純な関数になるから、 よって (18) の任意の $t$ に関する近似値が必要な場合は、 積分計算よりもむしろテイラー展開を利用する方が便利であると考えられる。 ただ、今回の場合はシンプソンの公式のような積分近似が有用ではないかというと そうではない。それを少し説明する。

(17), (18) を $r$, $\theta$ の パラメータ $t$ による表示だと考えて、 コンピュータにその曲線をプロットさせるには、 通常は以下のように行うだろう。

1. 十分 大きな $N$ を取って $0\leq t\leq 2\pi$ を $N$ 等分し、 その分点を $t_0=0$,$t_1$,$t_2$,...,$t_N=2\pi$ とする。
2. 各 $t_k$ に対する $r=r_k=r(t_k)$, $\theta=\theta_k=\theta(t_k)$ の値を 求め、$r_k$, $\theta_k$ の点と $r_{k-1}$, $\theta_{k-1}$ の点を 線分で結ぶ ( $k=1,2,\ldots,N$)。

これにより、曲線を折れ線で近似的に書くことになるわけであるが、 $N$ が大きければ、これでそれなりに滑らかに見える曲線が書ける。

この後者は、実際には gnuplot のようなデータプロットソフトを用いるとしても、 それらのソフトは内部では基本的には折れ線でデータ点を結ぶだけであり、 各 $k$ に対する $r_k$, $\theta_k$ のデータを取る必要があることは同じである。

つまり、今回の展開図を書くという目的の場合に必要なのは、

(A) 「任意の $t$ に対して (18) を計算すること」

なのではなく、

(B) 「数表のように細かい刻みに対する $t$ の値に対して (18) の値を順に計算すること」

であり、(B) の目的であればシンプソンの公式も十分有用である。

シンプソンの公式は、定積分を以下の式で近似する方法である。

$$\int_a^b f(x)dx \doteqdot \sum_{n=1}^N\frac{\Delta x}{2}\left\{\frac{1}{3} f(x_{n-1}) +\frac{4}{3} f(x_{n-1/2})+\frac{1}{3} f(x_n)\right\}$$

$$= \frac{\Delta x}{6}\{f(x_0)+4f(x_{1/2})+2f(x_1)+4f(x_{3/2})+2f(x_2)+\cdots$$

$$+2f(x_{N-1})+4f(x_{N-1/2})+f(x_N)\}$$ (20)

$$\left(\Delta x=\frac{b-a}{N},\hspace{1zw}x_j=a+j\Delta x\right)$$

つまり、$f(x)$ の、$x_{n-1}$ から $x_n(=x_{n-1}+\Delta x)$ までの積分を、

$$\frac{\Delta}{2}\left\{\frac{1}{3} f(x_{n-1}) +\frac{4}{3} f(x_{n-1/2})+\frac{1}{3} f(x_n)\right\}$$

で近似するものであるが、これは、$y=f(x)$ のグラフ上の 3 点

$$(x_{n-1},f(x_{n-1})),\hspace{1zw} (x_{n-1/2},f(x_{n-1/2})),\hspace{1zw} (x_n,f(x_n))$$

を通る 2 次関数の積分の値に等しく、台形公式よりも一般には精度は高い。

今、 $0\leq t\leq 2\pi$ を $N$ 等分して各 $t_k$ での $\theta_k=\theta(t_k)$ の 値を求める場合、上の $\Delta x$ を $\Delta t=2\pi/N$ と考え、 $\theta_k$ の計算には $k$ 等分のシンプソンの近似を用いて、

$$\theta_k=\int_0^{t_k}g_0(\phi)d\phi \doteqdot \sum_{n=1}^k\f... ...)+\frac{4}{3} g_0(t_{n-1/2}) +\frac{1}{3} g_0(t_n) \right\}$$

と近似するならば、$\theta_k$ の値は、$\theta_{k-1}$ の値を用いて、

$$\theta_k\doteqdot \theta_{k-1}+\frac{\Delta t}{6}\{g_0(t_{k-1})+4g_0(t_{k-1/2})+g_0(t_k)\}$$ (21)

のように計算されることになるから、 よって (B) の目的のためには、 (21) の後ろの項を加えていくだけで 容易に順に計算できるシンプソン近似もそれなりに有用なわけである。

次にテイラー展開を考えてみる。

$\theta(t)$ の微分は (19) であるから、 $g_0(t)$ の展開を求め、それを $t$ で積分すればよい。 もちろん、この $g_0(t)$ の $t$ での微分を直接順に計算できなくはないが、 かなり大変なので、これを合成関数とみて考えることにする。 今、

$$X=\frac{1-\cos t}{2}\hspace{1zw}(0\leq t\leq 2\pi)$$ (22)

とおくと $0\leq X\leq 1$ であり、$g_0(t)$ は、

$$g_1(X)=\frac{\sqrt{\hat{h}^2+4X^2}}{\hat{h}^2+4X}$$ (23)

と (22) の合成関数と見れる。 よって、まずこの $g_1(X)$ の展開を考える。

実は、$g_0(t)$ のマクローリン展開、 すなわち $t=0$ 中心の展開を考えるのはあまり得策ではない。 それは、$t=0$ は $X=0$ に対応し、 (23) の $X=0$ での展開があまりよくないからである。 (23) の $X$ を複素変数と考えると、 $g_1(X)$ の特異点は、

$$X=-\frac{\hat{h}^2}{4}, \pm\frac{\hat{h}}{2} i$$

なので、$\hat{h}$ が小さい場合は $g_1(X)$ のマクローリン展開の 収束半径が小さくなってしまって、$0\leq X\leq 1$ での 計算には使えない (図 8)。

図 8: $g_1(X)$ の収束円

しかし、$\hat{h}>0$ であるから、$X=1/2$ を中心とするテイラー展開を考えれば、 その収束半径は少なくとも $1/2$ よりも大きく、 よって $0\leq X\leq 1$ の範囲は全域がカバーされることになる。 しかも、$X=1/2$ は

$$X=\frac{1-\cos t}{2}=\frac{1}{2}$$

より、$t=\pi/2$, $3\pi/2$ という簡単な角度にも対応しているので、 計算にも都合がいいだろうと考えられる。

$g_1$ の $X=1/2$ でのテイラー展開を計算するには、 $Y=X-1/2$ として $g_1$ を $Y$ の式に書き直して、 それを $Y$ に関してマクローリン展開すればよい。

$$\begin{eqnarray*}g_1(X) &=& g_1(Y+1/2) = \frac{\sqrt{\hat{h}^2+4(Y+1/2)^2}}{... ...\hat{h}^2+1}Y(Y+1)} \left(1+\frac{4}{\hat{h}^2+2}Y\right)^{-1}\end{eqnarray*}$$

となるので、$\vert z\vert<1$ に対する一般二項定理

$$\left\{\begin{array}{ll} \sqrt{1+z} & \displaystyle = \sum... ...y}\right)z^n = \sum_{n=0}^\infty(-1)^nz^n \end{array}\right.$$ (24)

を用いれば展開の計算ができることになる。

なお、この (24) は本来は $\vert z\vert<1$ でしか成り立たないが、 $g_1(Y)$ の展開ではそれを無視して形式的に (24) を計算してよい。 それは、十分小さい $Y$ に対してはその条件を満たしているから 代入したものはそのような $Y$ に対しては正しい展開式になっていて、 上の特異点の考察によりその収束範囲がちゃんと保証されていて、 かつその係数は一意に決定するので、 その形式的な計算による級数が小さい $Y$ だけでなく その収束範囲まで収束してくれるからである。

ここでは簡単のため、$\hat{h}=2$、 すなわち斜円錐の高さが底面積の直径に等しい場合のみを考え、 $g_1$ の $Y$ に関する 6 次の展開式を求めることにする。 なお、$\hat{h}=2$ の場合は、$X=0$ での展開の収束半径も 1 になり、 $0\leq X<1$ までカバーされるのであるが、 $X=1$ での収束は不明で、もし収束したとしてもその近くでの収束は遅いので、 $X=1/2$ での展開を計算する方がよい。

$n\geq 1$ に対し、

$$\left(\begin{array}{c} 1/2 \ n \end{array}\right) =\frac{\d... ...displaystyle \frac{1}{n!}} =(-1)^{n-1} \frac{(2n-3)!!}{2^nn!}$$

であるから、6 次の近似式は

$$\begin{eqnarray*}\sqrt{1+z} &\doteqdot& 1+\frac{z}{2}-\frac{z^2}{2^3}+\frac{z^... ...ft(\frac{z}{2}\right)^5 -\frac{21}{16}\left(\frac{z}{2}\right)^6\end{eqnarray*}$$

となるので、$\hat{h}=2$ のとき、6 次の近似式を求めると、

$$\begin{eqnarray*}\lefteqn{\sqrt{1+\frac{4}{\hat{h}^2+1} Y(Y+1)}=\sqrt{1+\frac{4... ...^6 \ &=& 1+\bar{Y}+2\bar{Y}^2-2\bar{Y}^3+4\bar{Y}^5-6\bar{Y}^6\end{eqnarray*}$$

となる。ここで、$\bar{Y}=2Y/5$ とした。一方、

$$\begin{eqnarray*}\lefteqn{\left(1+\frac{4}{\hat{h}^2+2} Y\right)^{-1} =\left(1... ...)^4 -\left(\frac{2}{3} Y\right)^5+\left(\frac{2}{3} Y\right)^6\end{eqnarray*}$$

であるから、 $\hat{Y}=\bar{Y}/3=2Y/15$ とすると、 $\bar{Y}=3\hat{Y}$、 $2Y/3=5\hat{Y}$ より、

$$\begin{eqnarray*}\sqrt{1+\frac{4}{5} Y(Y+1)} &\doteqdot& 1+3\hat{Y}+18\hat{Y}... ...hat{Y}^2-125\hat{Y}^3+625\hat{Y}^4-3125\hat{Y}^5 +15625\hat{Y}^6\end{eqnarray*}$$

となるので、

$$\begin{eqnarray*}\lefteqn{\sqrt{1+\frac{4}{5} Y(Y+1)} \left(1+\frac{2}{3} Y\r... ...8\hat{Y}^2-194\hat{Y}^3+970\hat{Y}^4-3878\hat{Y}^5+15016\hat{Y}^6\end{eqnarray*}$$

となる。ここに、

$$\hat{Y} =\frac{2}{15} Y =\frac{2}{15}\left(X-\frac{1}{2}\ri... ...eft(\frac{1-\cos t}{2}-\frac{1}{2}\right) =-\frac{1}{15}\cos t$$

を代入すると、

$$g_0(t) \doteqdot \frac{\sqrt{5}}{6}\left(1+\sum_{j=1}^6 A_j\cos^j t\right)$$ (25)

$${\left(A_1=\frac{2}{15}, A_2=\frac{28}{15^2}, A_3=\frac{194}{... ...4=\frac{970}{15^4}, A_5=\frac{3878}{15^5}, A_6=\frac{15016}{15^6} \right)}$$

となる。

さらに、$\cos t$ を $t=\pi/2$ で展開し、

$$\cos t = -\sin\left(t-\frac{\pi}{2}\right) =-\left(t-\frac{... ...2}\right)^3 -\frac{1}{5!}\left(t-\frac{\pi}{2}\right)^5+\cdots$$

のように $t$ の多項式を代入すれば、 $g_0$ を完全に $t$ の多項式の形で展開でき、 それが $t=\pi/2$ でのテイラー展開となる。 しかし、それは途中で打ち切るわけであるからそこで多少の誤差が発生する。 それよりも、(25) を直接使うことを考えてみよう。

今、

$$I_n = I_n(t)=\int_0^t\cos^n\phi d\phi$$

とすると、(18), (25) より、

$$\theta(t) =\int_0^tg_0(\phi)d\phi \doteqdot\frac{\sqrt{5}}{6}\left(t+\sum_{j=1}^6A_jI_j\right)$$ (26)

となる。よって後は各 $I_n$ を求めればよいが、これは部分積分により、

$$\begin{eqnarray*}I_n &=& \int_0^t\cos^{n-1}\phi(\sin\phi)' d\phi = [\cos^{n... ...\ &=& \sin t\cos^{n-1}t+(n-1)(I_{n-2}-I_n)\hspace{1zw}(n\geq 2)\end{eqnarray*}$$

となるので、

$$I_n=\frac{n-1}{n} I_{n-2}+\frac{1}{n}\sin t\cos^{n-1}t$$

となる。よって、$I_0=t$, $I_1=\sin t$ より、

$$\begin{eqnarray*}I_2 &=& \frac{t}{2}+\frac{1}{2}\sin t\cos t,\\ I_3 &=& \... ...in t\cos t+\frac{5}{24}\sin t\cos^3 t +\frac{1}{6}\sin t\cos^5 t\end{eqnarray*}$$

となるので、(26) より、

$$\theta(t) \doteqdot \left(1+\frac{A_2}{2}+\frac{3}{8}A_4+\frac{5}{16}A_6\right)t +\left(A_1+\frac{2}{3}A_3+\frac{8}{15}A_5\right)\sin t$$

$$+\left(\frac{A_2}{2}+\frac{3}{8}A_4+\frac{5}{16}A_6\right)\sin t\cos t +\left(\frac{A_3}{3}+\frac{4}{15}A_5\right)\sin t\cos^2 t$$

$$+\left(\frac{A_4}{4}+\frac{5}{24}A_6\right)\sin t\cos^3 t +\frac{A_5}{5}\sin t\cos^4 t +\frac{A_6}{6}\sin t\cos^5 t$$ (27)

となる。

この、シンプソンの公式 (21) による計算と テイラー展開 (27) による計算を比較したのが 図 9 と図 10 である。

図 9: シンプソンとテイラー展開
図 10: シンプソンとテイラー展開の差

図 9 は、$N=100$ に対するシンプソンの計算と、 テイラー展開の $(t_k,\theta_k)$ の値を折れ線でつないだものであるが、 1 本のグラフに見えるようにほとんど違いがない。 図 10 はその違いを拡大して表示するために その両者の差をグラフにしたものであり、 その違いは最大で 0.00015 位であることがわかる。 この角度の差 0.00015 は、動径が 10000 のときに弧長が 1.5 ずれる程度の 誤差であるから、よほど大判のプロッタ出力でなければ出ない程度の差であろう。

また、$N=100$ はそれほど大きい分割数ではないが、$N=1000$ にしても 両者のグラフは $N=100$ の場合とほとんど違わない。 つまり、この差は量子誤差ではなく、 それぞれの計算法による本質的な違いであると考えられる。

なお、この差のグラフが $t=\pi/2$ と $t=3\pi/2$ の付近で平らになっているのは、 それがテイラー展開の中心であるからであり、 その値が 0 でないのは定積分のためにそこまでの誤差が累積しているからである。 そこからもわかるが、 シンプソンの公式もそれなりに一定して安定した精度を出しているようである。

結局いずれを採用してもそれなりに十分な精度で計算できると思われる。 テイラー展開の方を用いて、$N=100$, $R=3$, $h=6$ の場合の展開図が 図 11 である。

図 11: $R=3$, $h=6$ の展開図

両端 (最短の母線) は、$r=6$ であるが、 最長の部分は三平方の定理により $r=6\sqrt{2}\doteqdot 8.49$ となっている。 底辺の円は半径 $R=3$、つまりこの両端の線分の半分なので、 この図を拡大コピーして切りとって立体を作成することも可能である。 実際にそれを作ってみたのが図 12 である。

図 12: 完成図

若干作成時などの誤差はあるが、 ほぼ期待通りのものができていることがわかるだろう。