ベッセル関数の積分

質問

ベッセル関数を含んだ定積分に関する質問です．

∫（0→a）(r*exp(-2r^2/(R^2))*J0(kr))dr
J0は0次のベッセル関数です．

この定積分を行うとどのような形になるかわからず困っています．

よろしくお願いします．

回答No.6

inara1

＞このプログラムで試してみたいと思います
やってみてください。このプログラムでANo.2の図３と同じグラフが描けます。k*R の値はどれくらいなのですか。k は波数だと思いますが。

回答No.5

inara1

ANo.3のプログラムで改行が抜けてました。
×　deint = h / 3 * (A * Exp(-2 * A * A) * J0(kR * A) + 2 * s1 + 4 * s2) End Function
○　deint = h / 3 * (A * Exp(-2 * A * A) * J0(kR * A) + 2 * s1 + 4 * s2)
End Function

お礼

ご丁寧にありがとうございます．
このプログラムで試してみたいと思います．

回答No.4

inara1

[1] は見ることができないようなので、別の出典を紹介します。ここ（http://qe-forge.org/cgi-bin/cvstrac/q-e/fileview?f=espresso/PWCOND/bessj.f90&v=1.4）のページ中ほどの function bessj0(x) に出ています。

回答No.3

inara1

Excel VBA（マクロ）で計算する方法を書いておきます（末尾）。
第一種0次ベッセル関数は、エンジニアリング関数のアドインを有効にすればExcelでも計算できますが、その関数 BesselJ(x,0) はワークシート上でしか使えないので、マクロでは多項式による近似式 [1] を使いました。この近似式による計算結果は BesselJ(x,0) の値と同じです（同じ式を使っているのでしょう）。ベッセル関数の部分は、ワークシートのセルに =J0(x) と書けば計算できるので、 ワークシート関数 BesselJ(x,0) を使った場合と比較してみてください。x = 0～100の範囲で、真の値（数式処理ソフトで計算）と比較したところ、相対精度 [ (計算値-真値)/真値 ] は 10^(-6) 程度でした（特定の x でこれより精度が悪くなるところがある）。

問題の定積分を計算するユーザ関数は deint(A, kR) です。A のところには a/R の値を、kR には k*R の値を入れます。この関数の値は回答２の式(2)の右辺になりますので、問題の定積分の値は R^2*deint(a/R, k*R) になります。プログラム中の N の値は、定積分の範囲の分割数を 2 で割ったものです。N の値は k*R の大きさ によって変えるようにしてあります。k*R が大きいほどベッセル関数の振動周期が短くなるので N が大きくなるようにしています。しかし、k*R が 30 を超えると、N を1000より大きくしても、問題の定積分の相対精度は良くなりません（ベッセル関数の近似精度によるものと思われます）。ただし、問題の定積分の値そのものは、 k*R が大きいほどゼロに近づくので、絶対精度 [真値-計算値] は悪くありません。

【Excel2003でのマクロ（ユーザ定義関数）作成方法】
　（１）　Excelシートのメニューの [ツール] → [マクロ] → [Visual Basic Editor]
　（２）　Visual Vasic Editor のメニューの [挿入] → [標準モジュール]
　（３）　出てきた空白画面に以下の文を貼り付ける
↓ここから
' 定積分 ∫[0～a/R] x*exp(-2*x^2)*J0(k*R*x)dx
' A : a/R
' kR: k*R　(0<k*R)
'
Function deint(A As Double, kR As Double)
If kR < 0 Then
deint = ""
Exit Function
End If
Dim N As Integer
If kR < 1 Then
N = 100
Else
N = 100 * kR
End If
Dim M As Integer, i As Integer, h As Double, s1 As Double, s2 As Double, x1 As Double, x2 As Double
M = 2 * N ' 分割数
h = A / M ' 刻み幅
s1 = 0: s2 = 0
For i = 1 To N - 1
x1 = h * 2 * i
x2 = x1 - h
s1 = s1 + x1 * Exp(-2 * x1 * x1) * J0(kR * x1)
s2 = s2 + x2 * Exp(-2 * x2 * x2) * J0(kR * x2)
Next i
s2 = s2 + (A - h) * Exp(-2 * (A - h) * (A - h)) * J0(kR * (A - h))
deint = h / 3 * (A * Exp(-2 * A * A) * J0(kR * A) + 2 * s1 + 4 * s2) End Function

' 第一種0次ベッセル関数 J0(x)
'
Function J0(x As Double) As Double
Dim ax As Double, A As Double, a1 As Double, a2 As Double, z As Double, xx As Double
ax = Abs(x)
If ax < 8 Then
A = x * x
a1 = 57568490574# + (-13362590354# + (651619640.7 + (-11214424.18 + (77392.33017 - 184.9052456 * A) * A) * A) * A) * A
a2 = 57568490411# + (1029532985 + (9494680.718 + (59272.64853 + (267.8532712 + A) * A) * A) * A) * A
J0 = a1 / a2
Else
z = 8 / ax
A = z * z
xx = ax - 0.785398164
a1 = 1 + (-0.001098628627 + (0.00002734510407 + (-0.000002073370639 + 2.093887211E-07 * A) * A) * A) * A
a2 = -0.01562499995 + (0.0001430488765 + (-0.000006911147651 + (7.621095161E-07 - 9.34935152E-08 * A) * A) * A) * A
J0 = Sqr(0.636619772 / ax) * (Cos(xx) * a1 - z * Sin(xx) * a2)
End If
End Function
↑ここまで
　（４）　Visual Vasic Editor画面の右上隅の×をクリックするか、メニューの [ファイル] → [終了 Microsoft Excelへ戻る]
　（５）　シート上の適当なセルに =R^2*deint(a/R, k*R) と記入する（R、a、k は数値やセル位置）

【セキュリティーレベルの変更方法】
マクロ（Visual Basicプログラム）を含むExcelファイルは、PCのセキュリティーレベルによっては実行できないかもしれないので、以下にセキュリティーレベルの変更方法を書いておきます（Excel2003の場合）。
　（１）　Excelのブック（ワークシート）を新規作成
　（２）　メニューの[ツール] → [マクロ] → [セキュリティー]
　（３）　セキュリティーレベルを「中」に変更 → OK
　（４）　Excelを終了
　（５）　マクロを含むExcelファイルを開く
マクロを含むExcelファイルを開いた後、「このファイルはマクロを含んでいます」と警告が出ますが、[マクロを有効にする] をクリックすることでマクロを実行できるようになります。セキュリティーレベルが「高」のままだと、ファイルを開くことはできますがマクロを実行できません。

[1]　第一種0次ベッセル関数の近似式　 http://books.google.co.jp/books?id=0LSFEmmUFbQC&pg=PA349&lpg=PA349&dq=%2257568490574%22+%22bessel%22&source=bl&ots=4xXLNT2n9u&sig=4M3a2_dZYqa3mq2PNUA5u_G871E&hl=ja&ei=Hyp3S8C3O4vk7AOwoeiYBg&sa=X&oi=book_result&ct=result&resnum=3&ved=0CBMQ6AEwAg#v=onepage&q=%2257568490574%22%20%22bessel%22&f=false

回答No.2

inara1

数式処理ソフトでは
　　　a → ∞ なら、定積分 → (1/4)*R^2*exp{ -(1/8)*(k*R)^2 } --- (1)
でしたが、a が有限だと解析解は得られませんでした。

問題の定積分を I とすれば、x = r/R と置換することで
　　　I /R^2 = ∫[x=0～a/R] *x*exp(-2*x^2)*J0(k*R*x) dx --- (2)
と書くことができます。添付図の図１に示したように、x*exp(-2*x^2) は、2 < x で急速に小さくなります。J0(k*R*x) の値は、図２に示したように、k*R がどんな値であっても -0.5～1 の範囲にあるので、2*R < a ならば、式(2)の値はほぼ収束します（収束値は式(1)）。

図３は、a/R に対して、式(2)の I /R^2 がどのように変化するかを、いろいろな k*R について数値計算したものです。式(2)を a/R についてTylor展開して近似式を作ることは可能ですが、k*R が大きいときは、なかなか収束しません（非常に多くの項数が必要）。

2*R < a の場合は式(1)で近似できるのそれを使うことにして、a < 2*Rの場合は、数値積分で求めるのが良いのではないでしょうか（積分範囲は高々x = 0～ 2）。その場合、k*R が大きいほど J0(k*R*x) の振動周期が短くなるので分割数を多くする必要があります。

回答No.1

Ae610

∫（0→a）(r・exp(-2r^2/(R^2))・J0(kr))dr
Jn(x)=Σ[m=0～∞]((-1)^m/m!(n+m)!)・(x/2)^(n+2m)だから
n=0のとき
J0(kx)=Σ[m=0～∞]((-1)^m/m!・m!)・(kx/2)^(2m)
よって
∫（0→a）(r・exp(-2r^2/(R^2))・Σ[m=0～∞]((-1)^m/m!・m!)・(kr/2)^(2m)dr
=Σ[m=0～∞]((-1)^m・k^(2m)/(m!)^2・2^(2m))∫（0→a）(r^(2m+1)・exp(-2r^2/(R^2))dr
ここで
I=∫（0→a）(r^(2m+1)・exp(-2r^2/(R^2))drとおく。
r^2=tで置換すると
I=1/2・∫（0→a^2）t^m・exp(-bt)dt (b=2/R^2と置いた)
2I=I[m]とする。
I[m]=∫（0→a^2）t^m・exp(-bt)dt
=(-a^2m)/b・e^(-ba^2)+m/b・I[m-1]
=(-a^(2m))/b・e^(-ba^2)+m/b・(-a^(2m-2))/b・e^(-ba^2)+(m(m-1)/b^2)・I[m-2]
=(-a^(2m))/b・e^(-ba^2)+m/b・(-a^(2m-2))/b+m(m-1)/b^2・(-a^(2m-4))/b・e^(-ba^2)+(m(m-1)(m-2)/b^3)・I[m-3]
・・・
=-e^(-ba^2)・{Σ[n=1～m](m!/(m-n+1)!)・a^2(m-n+1)/b^n}I[0]
=e^(-ba^2)・(e^(-ba^2)-1)・{Σ[n=1～m](m!/(m-n+1)!)・a^2(m-n+1)/b^n}
よって
I=(e^(-ba^2)・(e^(-ba^2)-1)/2)・{Σ[n=1～m](m!/(m-n+1)!)・a^2(m-n+1)/b^n}

bを元に戻して
∫（0→a）(r・exp(-2r^2/(R^2))・J0(kr))dr
＝1/2・(e^(-2(a/R)^2)・(e^(-2(a/R)^2)-1)・{Σ[m=0～∞]((-1)^m・k^(2m)/(m!)^2・2^(2m))}{Σ[n=1～m](m!/(m-n+1)!)・a^2(m-n+1)・R^(2n)/(2^n)}

検算とかしてないので、表現式はちと自信ないのだが、方針としては
・J0(kr)を一端級数表現して積分する。
・置換積分した後、漸化式で表現
・I[m]を逐次計算
・元の積分式に代入
の手順で求めた。

お礼

ありがとうございます．この計算式で試してみたいと思います