関数cbessel マニュアル

(The documentation of function cbessel)

Last Update: 2023/11/14

◆機能・用途(Purpose)

複素変数に対するベッセル関数を計算する。
Compute the Bessel function for a complex variable.

ベッセル関数は

\begin{matrix} (1) & J_{n} (z) = {(\frac{z}{2})}^{n} \sum_{k = 0}^{\infty} \frac{(- 1)^{k} (z / 2)^{2 k}}{k! (n + k)!} \end{matrix}

によって定義され、漸化式

\begin{matrix} (2) & J_{m - 1} (z) + J_{m + 1} (z) = \frac{2 m}{z} J_{m} (z) \end{matrix}

によって計算できる。
The Bessel function is defined by Eq. (

1

) and can be computed recursively using Eq. (

2

◆形式(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)

この関数ではベッセル関数の漸化式(

2

)を後ろ向きに利用する(Millerの方法と呼ばれる)。はじめに十分に大きな整数

N

を選び、

J_{N + 1} (z) = 0

が成り立つものとする。

0 \leq m \leq N

について

\begin{matrix} (3) & j_{m} (z) \equiv \frac{J_{m} (z)}{J_{N} (z)} \end{matrix}

とおけば

\begin{matrix} (4) & j_{N + 1} (z) = 0 \end{matrix}

\begin{matrix} (5) & j_{N} (z) = 1 \end{matrix}

\begin{matrix} (6) & j_{m - 1} (z) = \frac{2 m}{z} j_{m} (z) - j_{m + 1} (z) (1 \leq m \leq N) \end{matrix}

が成り立つので、これを用いて

j_{N - 1} (z), j_{N - 2} (z), \dots, j_{0} (z)

を順次求める。最後に

\begin{matrix} (7) & J_{0} (z) + 2 \sum_{m = 1}^{\infty} J_{2 m} (z) = J_{N} (z) {j_{0} (z) + 2 \sum_{m = 1}^{\infty} j_{2 m} (z)} = 1 \end{matrix}

を用いて

J_{N} (z)

を決定すれば欲しい

n

に対する

J_{n} (z)

が得られる。
This function applies the recursive formula for the Bessel function (Eq.

2

) 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 (

3

), which satisfies Eqs. (

4

)-(

6

). From these equations,

j_{N - 1} (z), j_{N - 2} (z), \dots, j_{0} (z)

can be calculated one after another. Finally, Eq. (

7

) 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{a n})

とすることで有効桁数

\sqrt{a}

桁の計算が可能になるとのことである。この関数で用いるdouble型変数の有効桁数は15桁程度であるので

a = 15^{2} = 225

とすれば良い。但し、実際の計算精度は

z

にも依存すると思われるのでこの式をこのまま使うことには幾分か不安がある。そこで多少の余裕をみて

a = 20^{2} = 400

とするとともに、特に

n

の値が小さい場合に

N

が小さくなり過ぎないように最低でも

N = 100

以上とするという条件を加えた。また

z

を正の実数値とした計算では

z \geq N

になると数値不安定が起きることも分かったので

N

の値に

| z |

を加えるようにした。最終的に用いている計算式は

\begin{matrix} (8) & N = max {100, 2 (n + \sqrt{400 n})} + | z | \end{matrix}

である。
According to Press et al. (1988),

N = 2 (n + \sqrt{a n})

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. (

8

●ループ処理の高速化 (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.

ほかに必要になるのが(

7

)式に登場する

\begin{matrix} (9) & j_{0} (z) + 2 \sum_{m = 1}^{\infty} j_{2 m} (z) \end{matrix}

である。これを表す変数がJm_sumである。 Jm_sumは最初に0に初期化しておき、漸化式を後ろ向きに用いるループの中で偶数番目の

m

が登場するたびに

j_{m}

をJm_sumに加算していき、最後に全体を2倍して

j_{0}

を加えることでJm_sumを得る。
In addition, a variable defined by Eq. (

9

) is required to be preserved in order to apply the Eq. (

7

). 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}

で割る。これを行うと(

7

)式から計算される

J_{N}

の値が

10^{100}

倍されるが、

j_{n}

の値は

10^{- 100}

倍されているので (

3

)式から得られる

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. (

7

) 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. (

3

) 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)

Press WH, Vetterling WT, Teukolsky SA, Flannery BP (1988) Numerical Recipes in C, Cambridge University Press.

日本語訳(Japanese translation):
丹慶勝市, 佐藤俊郎, 奥村晴彦, 小林誠訳(1993) ニューメリカルレシピ・イン・シー日本語版—C言語による数値計算のレシピ, 技術評論社

関数cbessel マニュアル

◆機能・用途(Purpose)

◆形式(Format)

◆引数(Arguments)

◆戻り値(Return value)

◆使用例(Example)

◆計算方法(Computation method)

●概要 (Outline)

●Nの決定 (Determination of N)

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

●その他の変数 (Other variables)

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

◆検証(Validation)

◆引用文献(References)

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