数値計算をしていると、極座標を直交座標に変換したくなることがありますね
たとえば半径 $2$ の円上で角度 $60°$ にある点、すなわち極座標 $(2, 60°)$ をXY平面の座標に変換すると $(1, \sqrt{3})$ になる、みたいな話です。これは単純に、
\begin{align}
x &= r \cos(\phi) \\
y &= r \sin(\phi) \\
\end{align}
を計算すればいいですね。
同様に、球面の極座標 $(r, \phi_1, \phi_2)$ から $(x, y, z)$ を計算できます。
\begin{align}
x &= r \cos(\phi_1) \\
y &= r \sin(\phi_1) \cos(\phi_2) \\
z &= r \sin(\phi_1) \sin(\phi_2) \\
\end{align}
せっかくなので、N次元の超球面に拡張しましょう
\begin{align}
x_1 &= r \cos(\phi_1) \\
x_2 &= r \sin(\phi_1) \cos(\phi_2) \\
x_3 &= r \sin(\phi_1) \sin(\phi_2) \cos(\phi_3) \\
&\vdots \\
x_{n-1} &= r \sin(\phi_1) \dots \sin(\phi_{n-2}) \cos(\phi_{n-1}) \\
x_n &= r \sin(\phi_1) \dots \sin(\phi_{n-2}) \sin(\phi_{n-1}) \\
\end{align}
ということで、あとはこれを実装すればいいわけです。
$x_k$ は三角関数をK回くらい掛けると計算できるので、地道に二重ループ回してもいいのですが、しゅっと実装したいものですね。
とはいえ、ちょっと計算式がいかつく見えるかもしれません。分解してみましょう。
$x_k$ | $=$ | $r$ | $a_k$ | $b_k$ |
---|---|---|---|---|
$x_1$ | $=$ | $r$ | $1$ | $\cos(\phi_1)$ |
$x_2$ | $=$ | $r$ | $\sin(\phi_1)$ | $\cos(\phi_2)$ |
$x_3$ | $=$ | $r$ | $\sin(\phi_1) \sin(\phi_2)$ | $\cos(\phi_3)$ |
$\vdots$ | ||||
$x_{n-1}$ | $=$ | $r$ | $\sin(\phi_1) \dots \sin(\phi_{n-2})$ | $\cos(\phi_{n-1})$ |
$x_n$ | $=$ | $r$ | $\sin(\phi_1) \dots \sin(\phi_{n-2}) \sin(\phi_{n-1})$ | $1$ |
まずは $a_k$ の列を、 $\phi_k$ の列から生成する手順を考えてみます。
元列 | 要素ごとに三角関数 | 掛け算を累積 |
---|---|---|
$1$ | ||
$\phi_1$ | $\sin(\phi_1)$ | $\sin(\phi_1)$ |
$\phi_2$ | $\sin(\phi_2)$ | $\sin(\phi_1) \sin(\phi_2)$ |
$\vdots$ | ||
$\phi_{n-2}$ | $\sin(\phi_{n-2})$ | $\sin(\phi_1) \dots \sin(\phi_{n-2})$ |
$\phi_{n-1}$ | $\sin(\phi_{n-1})$ | $\sin(\phi_1) \dots \sin(\phi_{n-2}) \sin(\phi_{n-1})$ |
ここまで整理できれば、しめたものです。
Scalaだと map
して scan
すればいいので、次の1行で書けます。
phis.map(math.sin).scan(1.0)(_ * _)
同様に $b_k$ の列も、 $\phi_k$ の列から生成するためには、
| 元列 | 要素ごとに三角関数 |
| --- | --- | --- |
| $\phi_1$ | $\cos(\phi_1)$ |
| $\phi_2$ | $\cos(\phi_2)$ |
| $\vdots$ | | |
| $\phi_{n-2}$ | $\cos(\phi_{n-2})$ |
| $\phi_{n-1}$ | $\cos(\phi_{n-1})$ |
| | $1$ |
よって、次の1行で書けます。
phis.map(math.cos) :+ 1.0
あとは $a_k$ の列と $b_k$ の列を zip
して map
すればいいので、N次元の極座標から直交座標に変換する計算式は、このくらいのメソッドチェーンで実装できます。シンプル
def toEuclidean(r: Double, phis: Seq[Double]): Seq[Double] = {
phis.map(math.sin).scan(1.0)(_ * _).
zip(phis.map(math.cos) :+ 1.0).
map { case (a, b) => r * a * b }
}
ちなみに、N次元の直交座標から極座標に変換する計算式は、
\begin{align}
r &= \sqrt{x_n^2 + x_{n-1}^2 + \dots + x_1^2} \\
\phi_1 &= \arccos \frac{x_1}{\sqrt{x_n^2 + x_{n-1}^2 + \dots + x_1^2}} \\
\phi_2 &= \arccos \frac{x_2}{\sqrt{x_n^2 + x_{n-1}^2 + \dots + x_2^2}} \\
&\vdots \\
\phi_{n-2} &= \arccos \frac{x_{n-2}}{\sqrt{x_n^2 + x_{n-1}^2 + x_{n-2}^2}} \\
\phi_{n-1} &=
\begin{cases}
\arccos \frac{x_{n-1}}{\sqrt{x_n^2 + x_{n-1}^2}} &x_n \geq 0\\
- \arccos \frac{x_{n-1}}{\sqrt{x_n^2 + x_{n-1}^2}} &x_n < 0 \\
\end{cases}
\end{align}
ただし $x_n, x_{n-1}, \dots, x_k = 0$ の場合は、 $\phi_k = 0$ とする。
これも整理すると、 このくらいで実装できます。 scanRight
が綺麗にハマりますね。
def toSpherical(xs: Seq[Double]): (Double, Seq[Double]) = {
val rs = xs.scanRight(0.0) { case (x, sum) => x * x + sum }.map(math.sqrt)
val phis = xs.init.zip(rs).zipWithIndex.
map {
case ((_, 0.0), _) =>
0.0
case ((a, b), i) if i == xs.size - 2 && xs.last < 0.0 =>
-math.acos(a / b)
case ((a, b), _) =>
math.acos(a / b)
}
(rs.head, phis)
}
Scalaはコレクションメソッドが充実しているので、ちょちょっと組み合わせるだけで強力なコレクション処理が書けたりします。使いこなしていきましょう