前回のコラムに引き続き，今回も超幾何関数を用いたVoigt関数の計算を取り上げます．

　

　　Ｆ（α，β，γ：ｘ）＝１＋αβ／γｘ＋１／２！α（α＋１）β（β＋１）／γ（γ＋１）ｘ^2＋１／３！α（α＋１）（α＋２）β（β＋１）（β＋２）／γ（γ＋１）（γ＋２）ｘ^3＋・・・

はガウスの超幾何関数と呼ばれ，

　　2Ｆ1(α，β；γ：ｘ）

で表されます．ガウスの超幾何関数の収束半径は１であり，また，この級数は多くの初等関数，特殊関数（楕円積分など）を含んでいます．

　

　また，ガウスの超幾何微分方程式は３個の確定特異点をもち，それらは０，１，∞にスケール変換することができます．ここで，確定特異点１→∞として，２つの特異点を合流させると，その解として

　　Ｆ（α；γ：ｘ）＝１＋α／γｘ＋１／２！α（α＋１）／γ（γ＋１）ｘ^2＋１／３！α（α＋１）（α＋２）／γ（γ＋１）（γ＋２）ｘ^3＋・・・

が得られます．これはクンマーの合流形超幾何関数と呼ばれます．合流形超幾何関数は

　　1Ｆ1(α；γ：ｘ）

と書かれます．合流形超幾何関数は，ベッセル関数など数理物理に現れるいろいろな特殊関数を含んでいます．

　

　超幾何関数はオイラーによってすでに発見されていたのですが，ガウスはそれをもとに楕円関数論の研究を進めました．そしてクンマー，アーベル，ヤコビらの研究を発展させ代数関数論を完成させたのはリーマンです．超幾何関数論は数学的内容だけでなく，歴史的経緯も学ぶに値するものと考えられています．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　クンマーはガウスの超幾何関数を他の超幾何関数に移す変換，たとえば，

　　Ｆ（α，β，γ：ｘ）＝(1-x)^(γ-αｰβ)Ｆ（γ−α，γ−β，γ：ｘ）

　　Ｆ（α，β，γ：ｘ）＝(1-x)^(-α)Ｆ（α，γ−β，γ：ｘ／（１−ｘ））

などを構成しました．これがクンマーの関数等式と呼ばれるものですが，合流形超幾何関数では，次のようなクンマーの変換が知られています．

　　Ｆ（α；γ：−ｘ）＝ｅｘｐ（−ｘ）Ｆ（γ−α；γ：ｘ）

　

　超幾何関数の値を求める場合，級数を収束するまで加えればよいのですが，引数が負のときは交代級数ですから，桁落ちが起こるという問題があります．

　　∫(0,∞)exp(-ax^2-bx)cos(cx)dx

　　=Σ(ｰb)^k/k!/2a^((k+1)/2)Γ((k+1)/2)1F1((k+1)/2,1/2;-c^2/4a)

　

　この場合，クンマーの変換

　　1Ｆ1(a,c,-x)=exp(-x)1Ｆ1(c-a,c,x)

により

　　Σ(ｰb)^k/k!/2a^((k+1)/2)Γ((k+1)/2)exp(-c^2/4a)1F1(-k/2,1/2;c^2/4a)

とすれば，桁落ちが避けられます．

　

　さて，以下に掲げるプログラムは，前回紹介したプログラムにクンマー変換を施したものです．クンマー変換は精度を向上させるために行われるものですが，それによって収束半径も拡大されました．このことについては，次回以降で性能比較の項を設けてあらためて説明したいと思います．

　

1000 '

1010 ' *** Voigt 関数の計算（超幾何関数） ***

1020 '

1030 *VOIGT.INIT:

1040 '[in]

1050 'x : [argument]

1060 'p1: [weight]

1070 'p2: [location]

1080 'p3: [scale of Gaussian]

1090 'p4: [scale of Lorentzian]

1100 '

1110 '[out]

1120 'voigt#: f(x,p2,p3,p4)

1130 'v1#,v2#,v3#,v4#: Jacobian of p1*f(x,p2,p3,p4)

1140 'vd#: df/dx

1150 'vi#: integral(-∞,x)f(t)dt

1160 '

1170 PI#=3.141592653589732#

1180 QI#=SQR(PI#)

1190 '

1200 Q1=P3*P3/2

1210 Q2=P2-X

1220 Q3=P4

1230 UUX=Q2*Q2/Q1/4

1240 UUQ=ABS(Q2)

1250 Q4=EXP(-UUX)

1260 RETURN

1270 '

1280 ' *** Voigt 関数 ***

1290 '

1300 *VOIGT.HGF:

1310 GOSUB *VOIGT.INIT

1320 Z=UUX

1330 J=-2:GOSUB *VOIGT.HGF.SUB:G1#=G#

1340 J=-1:GOSUB *VOIGT.HGF.SUB:G2#=G#

1350 J=0 :　　　　　　　　　　:G3#=1#

1360 J=1 :GOSUB *VOIGT.HGF.SUB:G4#=G#

1370 '

1380 EPSL=1E-09

1390 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

1400 W1#= QI#/SQR(Q1)/2*Q4

1410 W2#=-1/Q1/2*Q3*Q4

1420 SS#=0:S0#=W1#*H3#+W2#*H4#

1430 J=0

1440 WHILE ABS((SS#-S0#)/S0#)>EPSL

1450　 A=-J/2　 :C=1/2

1460　 H5#=A/(C-A)*H1#-(A*2-C+Z)/(C-A)*H3#

1470　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*(1+J)/2

1480　 '

1490　 J=J+1

1500　 A=-J/2　 :C=1/2

1510　 H6#=A/(C-A)*H2#-(A*2-C+Z)/(C-A)*H4#

1520　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*(1+J)/2

1530　 '

1540　 H1#=H3#:H3#=H5#

1550　 H2#=H4#:H4#=H6#

1560　 SS#=S0#

1570　 S0#=S0#+W1#*H5#+W2#*H6#

1580　 J=J+1

1590 WEND

1600 VOIGT#=S0#/PI#

1610 RETURN

1620 '

1630 *VOIGT.HGF.SUB:

1640 A=-J/2　 :C=1/2

1650 GOSUB *CONFLUENT.HGF

1660 RETURN

1670 '

1680 ' *** ヤコビアン ***

1690 '

1700 *VOIGT.JACOBI1:

1710 GOSUB *VOIGT.HGF

1720 J1#=VOIGT#

1730 RETURN

1740 '

1750 ' *** ヤコビアン ***

1760 '

1770 *VOIGT.JACOBI2:

1780 GOSUB *VOIGT.INIT

1790 IF Q2=0 THEN J2#=0#:RETURN

1800 Z=UUX

1810 J=-2:GOSUB *VOIGT.JACOBI2.SUB:G1#=G#

1820 J=-1:GOSUB *VOIGT.JACOBI2.SUB:G2#=G#

1830 J=0 :　　　　　　　　　　　　:G3#=1#

1840 J=1 :GOSUB *VOIGT.JACOBI2.SUB:G4#=G#

1850 '

1860 EPSL=1E-09

1870 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

1880 W1#=-QI#/2/Q1^(3/2)*Q2/2*Q4

1890 W2#= 1/Q1^2*Q2/2*Q3*Q4

1900 SS#=0:S0#=W1#*H3#+W2#*H4#

1910 J=0

1920 WHILE ABS((SS#-S0#)/S0#)>EPSL

1930　 A=-J/2　 :C=3/2

1940　 H5#=A/(C-A)*H1#-(A*2-C+Z)/(C-A)*H3#

1950　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*(3+J)/2

1960　 '

1970　 J=J+1

1980　 A=-J/2　 :C=3/2

1990　 H6#=A/(C-A)*H2#-(A*2-C+Z)/(C-A)*H4#

2000　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*(3+J)/2

2010　 '

2020　 H1#=H3#:H3#=H5#

2030　 H2#=H4#:H4#=H6#

2040　 SS#=S0#

2050　 S0#=S0#+W1#*H5#+W2#*H6#

2060　 J=J+1

2070 WEND

2080 J2#=S0#/PI#*P1

2090 RETURN

2100 '

2110 *VOIGT.JACOBI2.SUB:

2120 A=-J/2　 :C=3/2

2130 GOSUB *CONFLUENT.HGF

2140 RETURN

2150 '

2160 ' *** ヤコビアン ***

2170 '

2180 *VOIGT.JACOBI3:

2190 GOSUB *VOIGT.INIT

2200 Z=UUX

2210 J=-2:　　　　　　　　　　　　:G1#=1#

2220 J=-1:GOSUB *VOIGT.JACOBI3.SUB:G2#=G#

2230 J=0 :GOSUB *VOIGT.JACOBI3.SUB:G3#=G#

2240 J=1 :GOSUB *VOIGT.JACOBI3.SUB:G4#=G#

2250 '

2260 EPSL=1E-09

2270 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

2280 W1#=-QI#/2/Q1^(3/2)/2*P3*Q4

2290 W2#= 1/Q1^2/2*Q3*P3*Q4

2300 SS#=0:S0#=W1#*H3#+W2#*H4#

2310 J=0

2320 WHILE ABS((SS#-S0#)/S0#)>EPSL

2330　 A=(-2-J)/2:C=1/2

2340　 H5#=A/(C-A)*H1#-(A*2-C+Z)/(C-A)*H3#

2350　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*(3+J)/2

2360　 '

2370　 J=J+1

2380　 A=(-2-J)/2:C=1/2

2390　 H6#=A/(C-A)*H2#-(A*2-C+Z)/(C-A)*H4#

2400　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*(3+J)/2

2410　 '

2420　 H1#=H3#:H3#=H5#

2430　 H2#=H4#:H4#=H6#

2440　 SS#=S0#

2450　 S0#=S0#+W1#*H5#+W2#*H6#

2460　 J=J+1

2470 WEND

2480 J3#=S0#/PI#*P1

2490 RETURN

2500 '

2510 *VOIGT.JACOBI3.SUB:

2520 A=(-2-J)/2:C=1/2

2530 GOSUB *CONFLUENT.HGF

2540 RETURN

2550 '

2560 ' *** ヤコビアン ***

2570 '

2580 *VOIGT.JACOBI4:

2590 GOSUB *VOIGT.INIT

2600 Z=UUX

2610 J=-2:GOSUB *VOIGT.JACOBI4.SUB:G1#=G#

2620 J=-1:　　　　　　　　　　　　:G2#=1#

2630 J=0 :GOSUB *VOIGT.JACOBI4.SUB:G3#=G#

2640 J=1 :GOSUB *VOIGT.JACOBI4.SUB:G4#=G#

2650 '

2660 EPSL=1E-09

2670 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

2680 W1#=-1/Q1/2*Q4

2690 W2#= QI#/2/Q1^(3/2)/2*Q3*Q4

2700 SS#=0:S0#=W1#*H3#+W2#*H4#

2710 J=0

2720 WHILE ABS((SS#-S0#)/S0#)>EPSL

2730　 A=(-1-J)/2:C=1/2

2740　 H5#=A/(C-A)*H1#-(A*2-C+Z)/(C-A)*H3#

2750　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*(2+J)/2

2760　 '

2770　 J=J+1

2780　 A=(-1-J)/2:C=1/2

2790　 H6#=A/(C-A)*H2#-(A*2-C+Z)/(C-A)*H4#

2800　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*(2+J)/2

2810　 '

2820　 H1#=H3#:H3#=H5#

2830　 H2#=H4#:H4#=H6#

2840　 SS#=S0#

2850　 S0#=S0#+W1#*H5#+W2#*H6#

2860　 J=J+1

2870 WEND

2880 J4#=S0#/PI#*P1

2890 RETURN

2900 '

2910 *VOIGT.JACOBI4.SUB:

2920 A=(-1-J)/2:C=1/2

2930 GOSUB *CONFLUENT.HGF

2940 RETURN

2950 '

2960 ' *** Voigt 関数の微分 ***

2970 '

2980 *VOIGT.DIFFER:

2990 GOSUB *VOIGT.INIT

3000 IF Q2=0 THEN VD#=0#:RETURN

3010 GOSUB *VOIGT.JACOBI2

3020 VD#=-S0#/PI#

3030 RETURN

3040 '

3050 ' *** Voigt 関数の積分 ***

3060 '

3070 *VOIGT.INTEGR:

3080 GOSUB *VOIGT.INIT

3090 IF Q2=0 THEN VI#=.5#:RETURN

3100 Z=UUX

3110 J=-2:GOSUB *VOIGT.INTEGR.SUB:G1#=G#

3120 J=-1:GOSUB *VOIGT.INTEGR.SUB:G2#=G#

3130 J=0 :GOSUB *VOIGT.INTEGR.SUB:G3#=G#

3140 J=1 :GOSUB *VOIGT.INTEGR.SUB:G4#=G#

3150 '

3160 EPSL=1E-09

3170 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

3180 W1#= QI#/SQR(Q1)/2*UUQ*Q4

3190 W2#=-1/Q1/2*Q3*UUQ*Q4

3200 SS#=0:S0#=W1#*H3#+W2#*H4#

3210 J=0

3220 WHILE ABS((SS#-S0#)/S0#)>EPSL

3230　 A=(2-J)/2:C=3/2

3240　 H5#=A/(C-A)*H1#-(A*2-C+Z)/(C-A)*H3#

3250　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*(1+J)/2

3260　 '

3270　 J=J+1

3280　 A=(2-J)/2:C=3/2

3290　 H6#=A/(C-A)*H2#-(A*2-C+Z)/(C-A)*H4#

3300　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*(1+J)/2

3310　 '

3320　 H1#=H3#:H3#=H5#

3330　 H2#=H4#:H4#=H6#

3340　 SS#=S0#

3350　 S0#=S0#+W1#*H5#+W2#*H6#

3360　 J=J+1

3370 WEND

3380 V#=S0#/PI#

3390 VI#=.5#-V#*SGN(Q2)

3400 RETURN

3410 '

3420 *VOIGT.INTEGR.SUB:

3430 A=(2-J)/2:C=3/2

3440 GOSUB *CONFLUENT.HGF

3450 RETURN

3460 '

3470 ' *** 合流型超幾何関数 ***

3480 '

3490 *CONFLUENT.HGF:' [Kummer]

3500 EPSL=1E-09

3510 G0#=0:G#=1

3520 T#=1:K=1

3530 WHILE ABS(G0#-G#)>EPSL

3540　 T#=T#*A/C*Z/K

3550　 G0#=G#

3560　 G#=G#+T#

3570　 K=K+1

3580　 A=A+1

3590　 C=C+1

3600 WEND

3610 G=G#

3620 RETURN

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝