Bessel（ベッセル）関数 - abwiki @ ウィキ (ActiveBasic非公式wiki)

ベッセルの微分方程式におけるy(x)の特異解である。

ベッセルの微分方程式は2階の線形微分方程式であるため、線形独立な2つの解が存在するが、その内n→∞で0に収束するものを第1種ベッセル関数と呼びJn(x)と表す。また逆にn→∞で発散するものを第2種ベッセル関数（もしくはノイマン関数ないしウェーバーの関数）とよびYn(x)とする。

第1種ベッセル関数は級数展開によっても定義されるが、変数xが大きいと桁落ちが生じる。そこで下記プログラムでは漸化式(1)を用いて求める

　Jn+1(x) = (2*n/x)Jn(x) - Jn-1(x)　･･･　(1)

十分大きな整数Nに対してJN(x)=ε、JN+1(x)=0とする初期条件よりJN-1(x)、JN-2(x)･･･を漸化式で求めることで目的の解を得る。出発点NはJ0(x)からJx(x)まではほぼ々大きさで、n>xではJn(x)=(x/2*n)Jn-1(x)で減衰するとし、JN(x)/Jn(x)<（許容相対誤差）であるように選択している。

第2種ベッセル関数に関してもほぼ同様で、Y0(x)およびY1(x)を求め、これを初期値として漸化式より解を得ている。

/***********************************************************
bessel.abp -- Bessel (ベッセル) 関数
***********************************************************/

#N88BASIC

Const EPS = 1e-10 ' 許容相対誤差
Const PI = 3.14159265358979324 ' 円周率(_System_PIでもよい)
Const EULER = 0.577215664901532861 ' Eulerの定数
'Const odd(x) = x And 1 ' 奇数判定
Function odd(x As Integer) As Integer
odd = x And 1
End Function

' n次第1種ベッセル関数
Function BesJ(n As Integer, x As Double) As Double
Dim k As Integer
Dim a As Double, b As Double, r As Double, s As Double
Dim x_2 As Double
x_2 = x / 2

' 変数xが負の場合
If x < 0 Then
If odd(x) Then
BesJ = -BesJ(n, -x)
Exit Function
Else
BesJ = BesJ(n, -x)
Exit Function
End If
End If
' 次数nが負の場合
If n < 0 Then
If odd(n) Then
BesJ = -BesJ(-n, x)
Exit Function
Else
BesJ = BesJ(-n, x)
Exit Function
End If
End If
'x = 0 である場合
If x = 0 Then
n = 0
BesJ = 0
Exit Function
End If

a = 0
s = 0
b = 1
k = n
' Jn(x)の収束する繰り返し数を求める
If k < x Then k = x
Do
k = k + 1
b = b * x_2 / k
Loop While b > EPS
If odd(k) Then k = k + 1 '奇数なら偶数にする
'
While k > 0
s = s + b
a = 2 * k * b / x - a 'a = J_k(x)
k = k - 1
If n = k Then r = a 'k = 奇数
b = 2 * k * a / x - b 'b = J_k(x)
k = k - 1
If n = k Then r = b 'k = 偶数
Wend
BesJ = r / (2 * s + b) 'J_0 + 2(J_2 + J_4 + ...) = 1 となるように規格化
End Function

' n次第1種ベッセル関数（級数展開版）
Function BesJ2(n As Integer, x As Double) As Double
Dim k As Integer
Dim result As Double, term As Double, previous As Double

' 変数xが負の場合
If x < 0 Then
If odd(x) Then
BesJ2 = -BesJ2(n, -x)
Exit Function
Else
BesJ2 = BesJ2(n, -x)
Exit Function
End If
End If
' 次数nが負の場合
If n < 0 Then
If odd(n) Then
BesJ2 = -BesJ2(-n, x)
Exit Function
Else
BesJ2 = BesJ2(-n, x)
Exit Function
End If
End If

x = x / 2
result = 1
For k = 1 To n
result = result * x / k
Next k
x = - x * x
term = result
For k = 1 To 499
term = term * x / (k * (n + k))
previous = result
result = result + term
If result = previous Then
BesJ2 = result
Exit Function
End If
Next k
Print "BesJ2(n, x): 収束しません."
BesJ2 = result
End Function

' n次第2種ベッセル関数
Function BesY(n As Integer, x As Double) As Double
Dim k As Integer
Dim a As Double, b As Double, s As Double, t As Double, u As Double
Dim x_2 As Double
Dim log_x_2 As Double
x_2 = x / 2
log_x_2 = Log(x_2)

' 変数xが負の場合
If x <= 0 Then
Print "BesY(n, x): x は正でなければなりません."
BesY = 0
Exit Function
End If
' 次数nが負の場合
If n < 0 Then
If odd(n) Then
BesY = -BesY(-n, x)
Exit Function
Else
BesY = BesY(-n, x)
Exit Function
End If
End If

k = x
b = 1
Do
k = k + 1
b = b * x_2 / k
Loop While b > EPS
If odd(k) Then k = k + 1 '奇数なら偶数にする
a = 0 'a = J_{k+1}(x) = 0, b = J_k(x), k 偶数
s = 0 '規格化の因子
t = 0 'Y_0(x)
u = 0 'Y_1(x)
While k > 0
s = s + b
t = b / (k / 2) - t
a = 2 * k * b / x - a
k = k - 1 'a = J_k(x), k 奇数
If k > 2 Then u = (k * a) / (Fix(k / 2) * Fix(k / 2 + 1)) - u
b = 2 * k * a / x - b
k = k - 1 'b = J_k(x), k 偶数
Wend
s = 2 * s + b
a = a / s
b = b / s
t = t / s
u = u / s 'a = J_1(x), b = J_0(x)
t = (2 / PI) * (2 * t + (log_x_2 + EULER) * b) 'Y_0(x)
If n = 0 Then
BesY = t
Exit Function
End If
u = (2 / PI) * (u + ((EULER - 1) + log_x_2) * a - b / x)'Y_1(x)
For k = 1 To n - 1
s = (2 * k) * u / x - t
t = u
u = s
Next k
BesY = u
End Function


Dim x As Integer

'第1種ベッセル関数の解
Print Ex"x\tJ_0(x)\t\t\tJ_1(x)\t\t\tJ_2(x)\t\t\tJ_3(x)"
For x = 0 To 20
Print x; Ex"\t"; BesJ(0, x), BesJ(1, x), BesJ(2, x), BesJ(3, x)
Next x
Print

'第2種ベッセル関数の解
Print Ex"x\tY_0(x)\t\t\tY_1(x)\t\t\tY_2(x)\t\t\tY_3(x)"
For x = 1 To 20
Print x; Ex"\t"; BesY(0, x), BesY(1, x), BesY(2, x), BesY(3, x)
Next x