関数cneuman マニュアル

(The documentation of function cneuman)

Last Update: 2023/11/14

◆機能・用途(Purpose)

複素変数に対するノイマン関数を計算する。
Compute the Neuman function for a complex variable.

ノイマン関数は

\begin{matrix} (1) & N_{n} (z) = lim_{ν \to n} \frac{\cos (ν π) J_{ν} (z) - J_{- ν} (z)}{\sin (ν π)} \end{matrix}

によって定義される。ここで

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

はベッセル関数、

\begin{matrix} (3) & Γ (z) = \int_{0}^{\infty} e^{- t} t^{z - 1} d t \end{matrix}

はガンマ関数である。
The Neuman function is defined by Eq. (

1

), where

J_{ν} (z)

is a Bessel function (Eq.

2

) and

Γ (z)

is a Gamma function (Eq.

3

◆形式(Format)

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

◆引数(Arguments)

n	ノイマン関数 $N_{n} (z)$ の次数 $n$ 。非負でなければならない。 The order $n$ of the Neuman function $N_{n} (z)$ , which must be non-negative.
z	ノイマン関数 $N_{n} (z)$ の引数 $z$ の値。 $z$ が実数の場合、 $0 < z < 35$ の範囲を推奨する。複素数値については適切な値の範囲をチェックできていない。 The value of the argument $z$ for the Neuman function $N_{n} (z)$ . For real $z$ , a range $0 < z < 35$ is recommended. An appropriate range for complex $z$ has not been examined.

◆戻り値(Return value)

引数で指定した

n

と

z

の組に対する

N_{n} (z)

の値。
The value of

N_{n} (z)

for

n

and

z

specified by the arguments.

◆使用例(Example)

double complex N=cneuman(1,0.5);

◆計算方法(Computation method)

●概要 (Outline)

n = 0, 1

に対しては級数表示

\begin{array}{rcl} N_{n} (z) & = & \frac{2}{π} J_{n} (z) \ln \frac{z}{2} \\ - \frac{1}{π} {(\frac{z}{2})}^{- n} \sum_{k = 0}^{n - 1} \frac{(n - k - 1)!}{k!} {(\frac{z}{2})}^{2 k} \\ (4) & - \frac{1}{π} {(\frac{z}{2})}^{n} \sum_{k = 0}^{\infty} \frac{(- 1)^{k} {(\frac{z}{2})}^{2 k}}{k! (n + k)!} [ψ (k) + ψ (k + n)] \end{array}

\begin{matrix} (5) & ψ (k) = - γ + \sum_{l = 1}^{k} \frac{1}{l} \end{matrix}

\begin{matrix} (6) & γ = 0.57721566 \end{matrix}

を用いる。

n ≧ 2

に対しては漸化式

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

を用いる。
For

n =

0 and 1, a series expansion formula given as Eqs. (

4

)-(

6

) is used. For

n ≧ 2

, a recursive formula (

7

) is used.

●計算式の簡単化 (Simplified formula)

(

4

)式を用いるのは

n = 0, 1

に対してだけであるので具体的に書き下すとより簡単になる。

n = 0

に対する式は

\begin{matrix} (8) & N_{0} (z) = \frac{2}{π} J_{0} (z) \ln \frac{z}{2} - \frac{2}{π} \sum_{k = 0}^{\infty} \frac{(- 1)^{k} {(\frac{z}{2})}^{2 k}}{(k!)^{2}} ψ (k) \end{matrix}

n = 1

に対する式は

\begin{matrix} (9) & N_{1} (z) = \frac{2}{π} J_{1} (z) \ln \frac{z}{2} - \frac{2}{π z} - \frac{z}{2 π} \sum_{k = 0}^{\infty} \frac{(- 1)^{k} {(\frac{z}{2})}^{2 k}}{k! (k + 1)!} [ψ (k) + ψ (k + 1)] \end{matrix}

である。
Because Eq. (

4

) is used only for

n =

0 and 1, explicit formulas for these

n

values are simpler than the original ones. The formulas for

n =

0 and 1 are given by Eqs. (

8

) and (

9

), respectively.

●無限級数の扱い (Treatment of the infinite series)

(

8

)(

9

)式には無限級数が登場する。これに対するプログラム内部の扱いとしては、新たに加わる項の大きさがそれまでの項の和の

10^{- 10}

倍を下回ったらループを抜けるという処理にしている。
An infinite series appears in Eqs. (

8

) and (

9

). In the program, the summation is taken until the newly added term becomes smaller than

10^{- 10}

times the summation.

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

(

8

)(

9

)式に登場する

(z / 2)^{2 k}

や

k!

などを別々に計算するとオーバーフローが起きる危険が増すので、このプログラムでは

ψ

に掛かる係数全体をfactorという名前の変数で表現している。この変数の初期値を1とし、

k

が1増えるごとにこの変数に

- (z / 2)^{2} / {k (k + n)}

を掛けるようにすればオーバーフローを回避しながらfactorを計算することが容易になる。
Individual computations of

(z / 2)^{2 k}

and

k!

in Eqs. (

8

) and (

9

) involve a high risk of overflow. To avoid this, the program uses a variable named “factor”, which represents the entire factor that is multiplied with

ψ

. This value is initialized at 1 and multiplied by

- (z / 2)^{2} / {k (k + n)}

for each

k

in the loop. In this way, the variable “factor” can easily be computed without a overflow.

◆使用上の注意(Note)

(z / 2)^{2 k}

と比較して

k

が小さいときには

k!

の大きくなるスピードは緩やかで、大きな

k

になるほど

k!

は相対的に急激に成長するようになる。その結果、(

8

)(

9

)式に登場する無限級数の項の大きさはある

k

までは増加の一途を辿り、

k

がある程度以上大きくなると減少に転じる。無限級数の各項には

(- 1)^{k}

が掛かっているので大きな項も打ち消し合って最終的には小さな値に落ち着く。但し、計算の途中であまりにも大きな項が登場すると打ち消し合う際に桁落ちが生じ、計算の精度が悪くなる。 double型変数の有効桁数が15桁程度であるので、最終的な和の値に比べて15桁以上大きな項同士が打ち消し合う状況が発生すると計算は事実上不可能になる。

z

の値が大きいほどこのようなことは起こりやすい。

z

として正の実数値を用いた数値計算の結果では

n

によらずおおよそ

z \sim 38

でこの問題が起きた。このことからこの関数では

| z | ≧ 35

の場合に警告を出す仕様にしている。しかしながら、

z

が複素数の場合にどのような

z

において問題が起きるのかは未検証である。
Compared to

(z / 2)^{2 k}

k!

relatively slowly increases for a small

k

but the growing speed of

k!

rapidly increases for a large

k

. As a result, a term in the infinite series in Eqs. (

8

) and (

9

) increases until a certain

k

, after which the term decreases. Since

(- 1)^{k}

appears in each term, large terms cancels out, resulting in a small final summation value. However, if a term appearing in the summation is too large, the summation results in a cancellation of significant digits, which worsens the accuracy of the computation result. Because a double-type variable has an accuracy of approximately 15 digits, if terms which are by 15 orders larger than the final summation value cancel out, the final result is completely a numerical noise. This situation more likely occurs for larger

z

. According to a numerical simulation using a positive real

z

, this problem occurred at

z \sim 38

regardless of the

n

value. Therefore this function warns if

| z | ≧ 35

. However, no test have been conducted for a complex

z

value.