関数cneuman マニュアル

(The documentation of function cneuman)

Last Update: 2023/11/14


◆機能・用途(Purpose)

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

ノイマン関数は (1)Nn(z)=limνncos(νπ)Jν(z)Jν(z)sin(νπ) によって定義される。ここで (2)Jν(z)=(z2)νk=0(1)k(z/2)2kk!Γ(ν+k+1) はベッセル関数、 (3)Γ(z)=0ettz1dt はガンマ関数である。
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 ノイマン関数Nn(z)の次数n。非負でなければならない。
The order n of the Neuman function Nn(z), which must be non-negative.
z ノイマン関数Nn(z)の引数zの値。 zが実数の場合、0<z<35の範囲を推奨する。 複素数値については適切な値の範囲をチェックできていない。
The value of the argument z for the Neuman function Nn(z). For real z, a range 0<z<35 is recommended. An appropriate range for complex z has not been examined.


◆戻り値(Return value)

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


◆使用例(Example)

double complex N=cneuman(1,0.5);


◆計算方法(Computation method)

●概要 (Outline)

n=0,1に対しては級数表示 Nn(z)=2πJn(z)lnz21π(z2)nk=0n1(nk1)!k!(z2)2k(4)1π(z2)nk=0(1)k(z2)2kk!(n+k)![ψ(k)+ψ(k+n)] (5)ψ(k)=γ+l=1k1l (6)γ=0.57721566 を用いる。n2に対しては漸化式 (7)Nm1(z)+Nm+1(z)=2mzNm(z) を用いる。
For n= 0 and 1, a series expansion formula given as Eqs. (4)-(6) is used. For n2, a recursive formula (7) is used.


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

(4)式を用いるのはn=0,1に対してだけであるので 具体的に書き下すとより簡単になる。n=0に対する式は (8)N0(z)=2πJ0(z)lnz22πk=0(1)k(z2)2k(k!)2ψ(k) n=1に対する式は (9)N1(z)=2πJ1(z)lnz22πzz2πk=0(1)k(z2)2kk!(k+1)![ψ(k)+ψ(k+1)] である。
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)式には無限級数が登場する。 これに対するプログラム内部の扱いとしては、 新たに加わる項の大きさがそれまでの項の和の1010倍を下回ったら ループを抜けるという処理にしている。
An infinite series appears in Eqs. (8) and (9). In the program, the summation is taken until the newly added term becomes smaller than 1010 times the summation.


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

(8)(9)式に登場する (z/2)2kk!などを別々に計算すると オーバーフローが起きる危険が増すので、 このプログラムではψに掛かる係数全体をfactorという名前の変数で表現している。 この変数の初期値を1とし、 kが1増えるごとにこの変数に(z/2)2/{k(k+n)}を掛けるようにすれば オーバーフローを回避しながらfactorを計算することが容易になる。
Individual computations of (z/2)2k 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)2kと比較してkが小さいときにはk!の大きくなるスピードは緩やかで、 大きなkになるほどk!は相対的に急激に成長するようになる。 その結果、(8)(9)式に登場する 無限級数の項の大きさはあるkまでは増加の一途を辿り、 kがある程度以上大きくなると減少に転じる。 無限級数の各項には(1)kが掛かっているので 大きな項も打ち消し合って最終的には小さな値に落ち着く。 但し、計算の途中であまりにも大きな項が登場すると打ち消し合う際に桁落ちが生じ、 計算の精度が悪くなる。 double型変数の有効桁数が15桁程度であるので、 最終的な和の値に比べて15桁以上大きな項同士が打ち消し合う状況が発生すると 計算は事実上不可能になる。 zの値が大きいほどこのようなことは起こりやすい。 zとして正の実数値を用いた数値計算の結果では nによらずおおよそz38でこの問題が起きた。 このことからこの関数では|z|35の場合に警告を出す仕様にしている。 しかしながら、zが複素数の場合にどのようなzにおいて 問題が起きるのかは未検証である。
Compared to (z/2)2k, 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 z38 regardless of the n value. Therefore this function warns if |z|35. However, no test have been conducted for a complex z value.