とつとつとしてろうとせず

ひまつぶしにどうぞ。

離散Hankel変換の計算方法

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

概要

この記事はYu, et al. (1998)およびGuizar-Sicairos and Gutierrez-Vega (2004)によって提案された準離散Hankel変換のアルゴリズムについてのまとめである.詳しくは当該論文を参照のこと.

はじめに

理工学の分野におけるある種の問題は,次の重調和方程式を解く問題に帰着される場合がある.

{
\begin{align}
 \nabla^{4} \phi = \nabla^2 \nabla^2 \phi = 0 
\end{align}
}

ここで,{\nabla^2}ラプラシアンである.典型的な例として,応力関数を用いた線形弾性問題の解法が挙げられる.たとえば,平面問題におけるAiryの応力関数が有名である(Wikipedia参照).境界条件を満たす重調和関数が求まれば,それを用いて応力や変位を計算することができる.

さて,円柱座標系{\paren{r,\theta,z}}における重調和方程式は次式のように表される.

{
\begin{align}
\label{eq:biharmonic}
\paren{ \frac{1}{r} \pdif{}{r} r \pdif{}{r} + \frac{1}{r^{2}} \pddif{}{\theta} + \pddif{}{z} }^{2} \phi(r,\theta,z) = 0 
\end{align}
}

この式を満たす重調和関数{\phi}を求めることを考える.まず,{\phi}をFourier余弦展開する.

{
\begin{align}
\phi(r,\theta,z) = \sum_{n=0}^{\infty} \phi_{n} (r,z) \cos n \theta 
\end{align}
}

上式を式{\eqref{eq:biharmonic}}に代入すれば,各Fourier係数について次式が成り立つ.

{
\begin{align}
\paren{ \frac{1}{r} \pdif{}{r} r \pdif{}{r} + \frac{1}{r^{2}} \pddif{}{\theta} + \pddif{}{z} }^{2} \phi_{n} (r,z) \cos n \theta  = 0 
\end{align}
}

{\theta}に関する微分を実行すれば次式を得る.

{
\begin{align}
\paren{ \frac{1}{r} \pdif{}{r} r \pdif{}{r} - \frac{n^{2}}{r^{2}} + \pddif{}{z} }^{2} \phi_{n} (r,z)  = 0 
\end{align}
}

ここで,上式に次式で定義されるHankel変換を適用する.

{
\begin{align}
\tilde{f}_{\nu}(s) = \mathcal{H}_{\nu} \bracket{ f(r) } = \int_{0}^{\infty} r f(r) J_{\nu}(s r) dr
\end{align}
}

ここで,{J_{\nu}(r)}は次数{\nu}の第一種Bessel関数である.Henkel変換の性質として次式が成り立つ.

{
\begin{align}
\mathcal{H}_{\nu} \bracket{ \paren{ \frac{1}{r} \dif{}{r} r \dif{}{r} - \frac{\nu^{2}}{r^{2}} } f(r) } = - s^{2} \tilde{f}_{\nu}(s)
\end{align}
}

したがって,重調和方程式は各Fourier係数に関する4階の常微分方程式に変換される.

{
\begin{align}
\paren{ \ddif{}{z} - s^{2} }^{2} \tilde{\phi}_{n}(s,z) = 0
\end{align}
}

これを解いて境界条件から未定係数を決定し,逆変換を適用すれば求める重調和関数を得ることができる.このような計算は一様厚さの無限厚板に軸対称な荷重が作用する場合などで必要となる。

準離散Hankel変換アルゴリズム

Hankel変換を計算するアルゴリズムとして,ここでは上述の文献のものを採用する.ただし,文献とHankel変換の定義式が異なるので注意されたい.
関数f(r){\nu}次のHankel変換と逆変換を次式で定義する.

{
\begin{align}
\tilde{f}_{\nu}(s) &= \int_{0}^{\infty} r f(r) J_{\nu}(s r) dr \\
f(r) &= \int_{0}^{\infty} s \tilde{f}_{\nu}(s) J_{\nu}(s r) ds 
\end{align}
}

数値計算では無限を扱えないため,積分計算をある点で打ち切り,関数の値は範囲外で0になるものとする.


\begin{align}
\label{eq:asm1}
 f(r \geq R) &= 0 \\
\label{eq:asm2}
 \tilde{f}_{\nu}(s \geq V) &= 0
\end{align}

ここで,Rは打ち切り半径,Vは打ち切り周波数である.すなわち,変換式を次のように近似する.

{
\begin{align}
\label{eq:tf}
\tilde{f}_{\nu}(s) &\approx \int_{0}^{R} r f(r) J_{\nu}(s r) dr \\
\label{eq:itf}
f(r) &\approx \int_{0}^{V} s \tilde{f}_{\nu}(s) J_{\nu}(s r) ds 
\end{align}
}

まず,f(r)\tilde{f}_{\nu}(s)をFourier-Bessel展開する(以前の記事ないし信頼できる適当な教科書を参照のこと).


\begin{align}
\tilde{f}_{\nu}(s) &= \frac{2}{V^{2}} \sum_{i=1}^{\infty} \bracket{ \frac{1}{J_{\nu+1}(\mu_{\nu i})^{2}} \int_{0}^{V} s' \tilde{f}_{\nu}(s') J_{\nu}\paren{ \frac{\mu_{\nu i} s'}{V} } ds' } J_{\nu}\paren{ \frac{\mu_{\nu i} s}{V} } \\
f(r) &= \frac{2}{R^{2}} \sum_{j=1}^{\infty} \bracket{ \frac{1}{J_{\nu+1}(\mu_{\nu j})^{2}} \int_{0}^{R} r' f(r') J_{\nu}\paren{ \frac{\mu_{\nu j} r'}{R} } dr' } J_{\nu}\paren{ \frac{\mu_{\nu j} r}{R} }
\end{align}

ここで,\mu_{\nu k}J_{\nu}(r)k番目の零点である.式\eqref{eq:tf}および\eqref{eq:itf}より,上式の係数部分に含まれる積分は次式のように近似できる.


\begin{align}
\int_{0}^{V} s' \tilde{f}_{\nu}(s') J_{\nu}\paren{ \frac{\mu_{\nu i} s'}{V} } ds' &\approx f\paren{ \frac{\mu_{\nu i}}{V} } \\
\int_{0}^{R} r' f(r') J_{\nu}\paren{ \frac{\mu_{\nu j} r'}{R} } dr' &\approx \tilde{f}_{\nu}\paren{ \frac{\mu_{\nu j}}{R} }
\end{align}

そこで,変換する関数のサンプリング点としてr_{i} = \mu_{\nu i} / Vs_{j} = \mu_{\nu j} / Rを選ぶと,Hankel変換の近似式として次式を得る.


\begin{align}
\label{eq:hankel}
\tilde{f}_{\nu}(s_{j}) &\approx \frac{2}{V^{2}} \sum_{i=1}^{\infty} \frac{J_{\nu}\paren{ \mu_{\nu i}\mu_{\nu j} / S }}{J_{\nu+1}(\mu_{\nu i})^{2}} f\paren{ \frac{\mu_{\nu i}}{V} } \\
\label{eq:ihankel}
f(r_{i}) &\approx  \frac{2}{R^{2}} \sum_{j=1}^{\infty} \frac{J_{\nu}\paren{ \mu_{\nu i}\mu_{\nu j} / S }}{J_{\nu+1}(\mu_{\nu j})^{2}} \tilde{f}_{\nu}\paren{ \frac{\mu_{\nu j}}{R} }
\end{align}

実際の計算では,無限和の打ち切り項数をNとして有限和に書き換え,次のようなベクトルを定義する.

\begin{align}
& G_{j} = \frac{ V }{ \abs{J_{\nu+1}(\mu_{\nu j})} } \tilde{f}_{\nu}\paren{ \frac{\mu_{\nu j}}{R} } \\
& H_{i} = \frac{ R }{ \abs{J_{\nu+1}(\mu_{\nu i})} } f\paren{ \frac{\mu_{\nu i}}{V} }
\end{align}

その結果,式\eqref{eq:hankel}および\eqref{eq:ihankel}の変換は次のベクトル-行列積で計算することができる.


\begin{align}
& G_{j} = \sum_{i=1}^{N} T_{ij} H_{i} \\
& H_{i} = \sum_{j=1}^{N} T_{ij} G_{j}
\end{align}

変換行列T_{ij}は次式で与えられる.


\begin{align}
T_{ij} = \frac{2}{S} \frac{ J_{\nu}(\mu_{\nu i} \mu_{\nu j} / S) }{ \abs{J_{\nu+1}(\mu_{\nu i})} \abs{J_{\nu+1}(\mu_{\nu j})} }
\end{align}

ここで, S=RVである.T_{ij}については次の性質が成り立つ.

  1. 実対称行列である.すなわち,T_{ij}=T_{ji}
  2. 直交行列である.すなわち,\displaystyle \sum_{k=1}^{N} T_{ik} T_{jk} = \delta_{ij}

しかし,実際にはT_{ij}=T_{ij}(S)であり,上の性質はパラメータSの値に依存する.参考文献によれば,S = \mu_{\nu,N+1}とおけば十分大きなNに対して上記の性質は満たされる.たとえば,逆変換の計算を行う場合,実空間の大きさRが決っているならば,波数空間の大きさはV = S / Rから求めればよい.
離散データに対する計算の注意点として,変換する関数の評価点がBessel関数の零点でなければならないことが挙げられる.解析的な表現が得られている場合は問題ないが,等間隔の格子点上の数値データを変換するさいには,線形補間などによって評価点を移動する必要がある.

まとめ

離散Hankel変換の計算ライブラリはGNU Scientific Library (GSL)に含まれている.後から気づいたのだが,GSLのアルゴリズムはこの記事で紹介したものとさほど大きな違いはなさそうである(というか違いがあったとしても門外漢なのでわからない).特にこだわりがなければGSLのものを利用すればよいだろう.ただ,計算効率や作業時間よりも使いやすさを優先する場合や勉強の一環で中身の理解が重要だということであれば自分でプログラムを組んでみるのも一興である(というか俺はそうしてしまった).見た目ほど(?)複雑なプログラムではないのでお時間のある方はどうぞ(なお,自分の場合は行列計算ライブラリとしてEigen3,特殊関数ライブラリとしてboostのものを使用した).