関数cbessel マニュアル

(The documentation of function cbessel)
Last Update: 2023/11/14


◆機能・用途(Purpose)

複素変数に対するベッセル関数を計算する。
Compute the Bessel function for a complex variable.
ベッセル関数は \[\begin{equation} J_n(z)=\left(\frac{z}{2}\right)^n \sum_{k=0}^{\infty}\frac{(-1)^k (z/2)^{2k}}{k!(n+k)!} \label{eq.Bessel} \end{equation}\] によって定義され、漸化式 \[\begin{equation} J_{m-1}(z)+J_{m+1}(z)=\frac{2m}{z}J_m(z) \label{eq.bessel.recursive} \end{equation}\] によって計算できる。
The Bessel function is defined by Eq. (\ref{eq.Bessel}) and can be computed recursively using Eq. (\ref{eq.bessel.recursive}).

◆形式(Format)

#include <mathfunc.h>
inline double complex cbessel(const int n,const double complex z)


◆引数(Arguments)

n ベッセル関数\(J_n(z)\)の次数\(n\)。非負でなければならない。
The order \(n\) of the Bessel function \(J_n(z)\), which must be non-negative.
z ベッセル関数\(J_n(z)\)の引数\(z\)の値。
The value of the argument \(z\) for the Bessel function \(J_n(z)\).


◆戻り値(Return value)

引数で指定した\(n\)と\(z\)の組に対する\(J_n(z)\)の値。
The value of \(J_n(z)\) for \(n\) and \(z\) specified by the arguments.

◆使用例(Example)

double complex J=cbessel(1,0.5);


◆計算方法(Computation method)

●概要 (Outline)

この関数ではベッセル関数の漸化式(\ref{eq.bessel.recursive})を 後ろ向きに利用する(Millerの方法と呼ばれる)。 はじめに十分に大きな整数\(N\)を選び、\(J_{N+1}(z)=0\)が成り立つものとする。 \(0 \leq m \leq N\)について \[\begin{equation} j_m(z)\equiv \frac{J_m(z)}{J_N(z)} \label{eq.jm} \end{equation}\] とおけば \[\begin{equation} j_{N+1}(z)=0 \label{eq.jN1} \end{equation}\] \[\begin{equation} j_N(z)=1 \label{eq.jN} \end{equation}\] \[\begin{equation} j_{m-1}(z)=\frac{2m}{z}j_m(z)-j_{m+1}(z) \hspace{1em} (1 \leq m \leq N) \label{eq.jm.recursive} \end{equation}\] が成り立つので、これを用いて \(j_{N-1}(z), j_{N-2}(z), \cdots, j_0(z)\) を順次求める。最後に \[\begin{equation} J_0(z)+2\sum_{m=1}^{\infty}J_{2m}(z) =J_N(z)\left\{j_0(z)+2\sum_{m=1}^{\infty}j_{2m}(z)\right\}=1 \label{eq.bessel.normalization} \end{equation}\] を用いて\(J_N(z)\)を決定すれば欲しい\(n\)に対する\(J_n(z)\)が得られる。
This function applies the recursive formula for the Bessel function (Eq. \ref{eq.bessel.recursive}) in the backward order (known as Miller's algorithm). First, the program chooses a sufficiently large integer number \(N\), for which \(J_{N+1}(z)=0\) is assumed. For \(0 \leq m \leq N\), define \(j_m(z)\) by (\ref{eq.jm}), which satisfies Eqs. (\ref{eq.jN1})-(\ref{eq.jm.recursive}). From these equations, \(j_{N-1}(z), j_{N-2}(z), \cdots, j_0(z)\) can be calculated one after another. Finally, Eq. (\ref{eq.bessel.normalization}) is applied to estimate \(J_N(z)\), which enables us to compute \(J_n(z)\) for the specified \(n\).

●\(N\)の決定 (Determination of \(N\))

Press et al. (1988)によれば、\(N=2(n+\sqrt{an})\)とすることで 有効桁数\(\sqrt{a}\)桁の計算が可能になるとのことである。 この関数で用いるdouble型変数の有効桁数は15桁程度であるので \(a=15^2=225\)とすれば良い。 但し、実際の計算精度は\(z\)にも依存すると思われるので この式をこのまま使うことには幾分か不安がある。 そこで多少の余裕をみて\(a=20^2=400\)とするとともに、 特に\(n\)の値が小さい場合に\(N\)が小さくなり過ぎないように 最低でも\(N=100\)以上とするという条件を加えた。 また\(z\)を正の実数値とした計算では\(z \geq N\)になると 数値不安定が起きることも分かったので\(N\)の値に\(|z|\)を加えるようにした。 最終的に用いている計算式は \[\begin{equation} N=\max\{100,2(n+\sqrt{400n})\}+|z| \label{eq.N} \end{equation}\] である。
According to Press et al. (1988), \(N=2(n+\sqrt{an})\) enables a computation in an accuracy of approximately \(\sqrt{a}\)-digits. Because a double-type variable, used in this function, has an accuracy of about 15 digits, appropriate \(a\) may be estimated as \(a=15^2=225\). However, the accuracy may depend on \(z\), which is not considered in the reference. To safely avoid errors taking into account this uncertainty, I used \(a=20^2=400\) instead, and required the minimum \(N\) value to be 100. In addition, a numerical test using a positive real \(z\) value suggested that a numerical instability occured when \(z \geq N\). Therefore \(|z|\) was added to \(N\). The final equation used in the program is Eq. (\ref{eq.N}).

●ループ処理の高速化 (Coding for a faster loop)

プログラムの核となるのは漸化式を後ろ向きに用いる部分である。 この部分に対する単刀直入なコーディングとしては以下のCode 1が考えられる。 ここでJmm,Jm,Jmpはそれぞれ\(j_{m-1}\), \(j_m\), \(j_{m+1}\)を表す。 The main part of the program is the backward application of the recursive formula. A direct coding for this part may be the Code 1 below, where “Jmm”, “Jm”, and “Jmp” represent \(j_{m-1}\), \(j_m\), and \(j_{m+1}\), respectively.
(Code 1)
for(m=N;m>=1;m++){
    Jmm=2.0∗m∗Jm/z-Jmp;
    Jmp=Jm;
    Jm=Jmm;
}

Code 1ではループ内の2番目と3番目は単なる変数のコピーであり、 変数をうまく定義することによって削ることができる。 この改良を行ったのがCode 2であり、 連続した3つの\(m\)について続けて漸化式を適用することによって 不要な変数コピー処理を回避して高速化を実現している。 数値テストの結果、この改良によって 関数内の他の処理も含めたトータルの処理時間が 2/3倍程度まで削減できることが確認できた。 Code 2におけるJm1,Jm2,Jm3の定義は 表1に示した通りである。 この関数ではCode 2を基本的な枠組みとしてコーディングを行っている。
The second and third processings in the loop in “Code 1” are simply a copying operation, which can be removed by better defining the variables. The result is “Code 2”, in which the recursive formulas for three consecutive \(m\) values are consequently applied in a loop. This enables a faster processing by avoiding unnecessary variable-copying operations. A numerical test indicated that this improvement reduces the total processing time, including the other processings in the function, can be reduced to about two-thirds. Definitions of “Jm1”, “Jm2”, and “Jm3” are listed in Table 1. The “Code 2” is used as a basic frame of this function code.
(Code 2)
for(m=N-1;m>=3;m-=3){
    Jm3=2.0∗(m+1)∗Jm1/z-Jm2;
    Jm2=2.0∗m∗Jm3/z-Jm1;
    Jm1=2.0∗(m-1)∗Jm2/z-Jm3;
}

表1. 変数Jm1,Jm2,Jm3の定義。
Table 1. The definitions of variables “Jm1”, “Jm2”, and “Jm3”.
変数
Variable 		ループ前
Before loop 		ループ後
After loop
Jm3\(j_{m+3}\)\(j_m\)
Jm2\(j_{m+2}\)\(j_{m-1}\)
Jm1\(j_{m+1}\)\(j_{m-2}\)


●その他の変数 (Other variables)

この関数では全ての\(m\)に対する\(j_m\)を保存することはしない。 欲しい\(n\)に対する\(j_n\)のみを保存しておく。 これを保存するための変数がJnである。
In this function, \(j_m\) for each \(m\) is not preserved; only the \(j_n\) for the required \(n\) is preserved. For this purpose, a variable “Jn” is used.
ほかに必要になるのが(\ref{eq.bessel.normalization})式に登場する \[\begin{equation} j_0(z)+2\sum_{m=1}^{\infty}j_{2m}(z) \label{eq.Jmsum} \end{equation}\] である。これを表す変数がJm_sumである。 Jm_sumは最初に0に初期化しておき、 漸化式を後ろ向きに用いるループの中で 偶数番目の\(m\)が登場するたびに\(j_m\)をJm_sumに加算していき、 最後に全体を2倍して\(j_0\)を加えることでJm_sumを得る。
In addition, a variable defined by Eq. (\ref{eq.Jmsum}) is required to be preserved in order to apply the Eq. (\ref{eq.bessel.normalization}). For this purpose, a variable “Jm_sum” is used. This variable is initialized to zero, added by \(j_m\) in the loop with regard to \(m\) whenever an even \(m\) value appears, and finally multiplied by 2 and added by \(j_0\) to obtain the final value of “Jm_sum”.

●オーバーフローの回避 (Avoiding over flows)

\(m\)が小さくなるにつれて\(j_m\)の値は急激に大きくなるので 簡単にオーバーフローが起きてしまう。 これを回避するために関数では\(j_m\)が\(10^{100}\)を超えた時点で \(j_m\), \(j_n\), およびJm_sumを全て\(10^{100}\)で割る。 これを行うと(\ref{eq.bessel.normalization})式から計算される\(J_N\)の値が \(10^{100}\)倍されるが、\(j_n\)の値は\(10^{-100}\)倍されているので (\ref{eq.jm})式から得られる\(J_n\)の値は変わらない。
Since \(j_m\) rapidly increases when \(m\) decreases, a overflow easily occurs. To avoid this, \(j_m\), \(j_n\), and “Jm_sum” are all divided by \(10^{100}\) when \(j_m\) exceeds \(10^{100}\). By this operation, the \(J_N\) value estimated by Eq. (\ref{eq.bessel.normalization}) is multiplied by \(10^{100}\) but at the same time \(j_n\) is multiplied by \(10^{-100}\). Therefore, the \(J_n\) value estimated by Eq. (\ref{eq.jm}) does not change by this operation.

◆検証(Validation)

\(z\)として正の実数値を利用した計算では \(0 \leq n \leq 3\), \(0.01 \leq z \leq 1000\)の範囲において C言語標準関数であるj0, j1, jnを用いた計算の結果と非常に良く一致した。
For positive real \(z\) values, excellent fits with results of computations using C standard functions j0, j1, and jn were obtained for \(0 \leq n \leq 3\) and \(0.01 \leq z \leq 1000\).

◆引用文献(References)