フォークト関数

　　h(z)=∫(0,∞)exp(-at^2-bt)cos{(c-z)t}dt

は，超幾何関数や放物柱関数を用いて精確に求めることができます．放物柱関数もまた合流形超幾何関数に還元されますから，最終的に超幾何関数に帰着されることになります．超幾何関数を用いたこれらの計算法は，計算しやすさまで考慮に入れても適当かつ有用と思われました．

　

　一方，フォークト関数は，複素誤差関数を用いて

　　h(z)=1/4√(π/a){exp((b+(c-x)i)^2/4a)erfc(b+(c-x)i)/2√a)+exp((b-(c-x)i)^2/4a)erfc(b-(c-x)i)/2√a)}

で表されます（新数学公式集（丸善）のp452,(2.5.36.8)式）．

　

　したがって，当該の問題は

　　exp(z^2)erfc(z)

　　z=x+yi

　　erfc:complementary error function

の実数部分Reと虚数部分Imが求められれば解決できることが理解されます．そして，ｗをｚの共役複素数とすると，

　　exp(z^2)erfc(z)+exp(w^2)erfc(w)=Re*2

　　exp(z^2)erfc(z)-exp(w^2)erfc(w)=Im*2

と表すことができます．

　

　今回のコラムでは，最後に残された道〈複素誤差関数〉を用いる方法にトライしてみることにしますが，積分の経路など複素関数の扱い方に慣れていたとしても，ある工夫を加えない限り収束しませんでした．そのため，自分でもあっと驚くようなトリックが用いられていますし，また，ヤコビアンの計算では放物柱関数の複素化が必要になりました．

　

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

　

【１】複素誤差関数を用いたVoigt関数の計算

　

　　F(z)=exp(z^2)erfc(z)=exp(z^2){1-erf(z)}

は第２種合流型超幾何関数で表すと

　　F(z)=U(1/2,1/2,z^2)/√π

ですが,実際の計算のために，

　　F(z)=exp(z^2)-1F1(1,3/2,z^2)*2z/√π

と２項に分けるとうまく収束しませんでした．

　

　しかし，この式はexp(z^2)がｚの偶数次,1F1(1,3/2,z^2)*2z/√πが奇数次に関係していて，並べ替えによって

　　F(z)=Σ(-z)^n/Γ(n/2+1)

と統一した形に表すことができます．なかなか気づきにくい変形です．

　

　あとは，それに

　　z=r(cosθ+isinθ)

を代入すれば，簡単に実部と虚部に分離でき，今度はうまく収束させることができました．その際，三角関数の加法定理が必要になります（2190〜2200行，2240〜2250行）．

　

　このアイデアのおかげで，これまで掲げた４つの候補（超幾何関数を用いるもの，放物柱関数を用いるものなど）のなかで，最も短いプログラムになっています．

　

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 '

1270 UUX=P4*P4/P3/P3/2

1280 Q=Q2/P3

1290 QQ=Q*Q

1300 '

1310 XI#=Q3/SQR(Q1)/2

1320 YI#=Q2/SQR(Q1)/2

1330 ZQ#=XI#*XI#+YI#*YI#

1340 ZI#=SQR(ZQ#)

1350 CZ#=XI#/ZI#

1360 SZ#=YI#/ZI#

1370 RETURN

1380 '

1390 ' *** Voigt 関数 ***

1400 '

1410 *VOIGT.HGF:

1420 GOSUB *VOIGT.INIT

1430 GOSUB *COMPLEX.ERF

1440 SS#=RE#*SQR(PI#/Q1)/2

1450 VOIGT#=SS#/PI#

1460 RETURN

　

2020 '

2030 ' *** 複素誤差関数 ***

2040 '

2050 *COMPLEX.ERF:

2060 EPSL=1E-09

2070 C0#=1#:S0#=0#

2080 W1#= 1#

2090 RE#=0:RE0#=W1#*C0#

2100 IM#=0:IM0#=W1#*S0#

2110 '

2120 C0#=CZ#:S0#=SZ#

2130 W2#=-2#/QI#*ZI#

2140 RE0#=RE0#+W2#*C0#

2150 IM0#=IM0#+W2#*S0#

2160 '

2170 J=2

2180 WHILE ABS(RE#-RE0#)>EPSL OR ABS(IM#-IM0#)>EPSL

2190　 C1#=CZ#*C0#-SZ#*S0#

2200　 S1#=SZ#*C0#+CZ#*S0#

2210　 W1#=W1#/J*2*ZQ#

2220　 '

2230　 J=J+1

2240　 C2#=CZ#*C1#-SZ#*S1#

2250　 S2#=SZ#*C1#+CZ#*S1#

2260　 W2#=W2#/J*2*ZQ#

2270　 '

2280　 C0#=C2#:S0#=S2#

2290　 RE#=RE0#:IM#=IM0#

2300　 RE0#=RE0#+W1#*C1#+W2#*C2#

2310　 IM0#=IM0#+W1#*S1#+W2#*S2#

2320　 J=J+1

2330 WEND

2340 RETURN

　

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

　

【２】ヤコビアンの計算

　

　次に，ヤコビアンの計算に際してですが，複素放物柱関数ともいうべき，

　　F(z)=exp(z^2/4)D-3(z)

　　　　=2^(-3/2)U(3/2,1/2,z^2/2)

を考え出す必要が生じました．

　

　この式は，複素誤差関数と同様の並べ替えによって，

　　F(u)=Σ(n+1)(-u)^n/Γ(n/2+1)

　　u=z/√2

と表すことができます．

　

1470 '

1480 ' *** ヤコビアン ***

1490 '

1500 *VOIGT.JACOBI1:

1510 GOSUB *VOIGT.HGF

1520 J1#=VOIGT#

1530 RETURN

1540 '

1550 ' *** ヤコビアン ***

1560 '

1570 *VOIGT.JACOBI2:

1580 GOSUB *VOIGT.INIT

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

1600 GOSUB *COMPLEX.ERF

1610 SS#=-(RE#*Q2+IM#*Q3)*SQR(PI#/Q1)/Q1/4

1620 J2#=SS#/PI#*P1

1630 RETURN

1640 '

1650 ' *** ヤコビアン ***

1660 '

1670 *VOIGT.JACOBI3:

1680 GOSUB *VOIGT.INIT

1690 GOSUB *COMPLEX.PCF

1700 SS#=-RE#/(Q1*2)^(3/2)*2*P3

1710 J3#=SS#/PI#*P1

1720 RETURN

1730 '

1740 ' *** ヤコビアン ***

1750 '

1760 *VOIGT.JACOBI4:

1770 GOSUB *VOIGT.INIT

1780 GOSUB *COMPLEX.ERF

1790 SS#=(RE#*Q3-IM#*Q2)*SQR(PI#/Q1)/Q1/4-1/Q1/2

1800 J4#=SS#/PI#*P1

1810 RETURN

1820 '

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

1840 '

1850 *VOIGT.DIFFER:

1860 GOSUB *VOIGT.INIT

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

1880 GOSUB *VOIGT.JACOBI2

1890 VD#=-SS#/PI#

1900 RETURN

　

2350 '

2360 ' *** 複素放物柱関数 ***

2370 '

2380 *COMPLEX.PCF:

2390 EPSL=1E-09

2400 C0#=1#:S0#=0#

2410 W1#= QI#/2^(3/2)

2420 RE#=0:RE0#=W1#*C0#

2430 IM#=0:IM0#=W1#*S0#

2440 '

2450 C0#=CZ#:S0#=SZ#

2460 W2#=-SQR(2)*ZI#

2470 RE0#=RE0#+W2#*C0#

2480 IM0#=IM0#+W2#*S0#

2490 '

2500 J=2

2510 WHILE ABS(RE#-RE0#)>EPSL OR ABS(IM#-IM0#)>EPSL

2520　 C1#=CZ#*C0#-SZ#*S0#

2530　 S1#=SZ#*C0#+CZ#*S0#

2540　 W1#=W1#/J*2*ZQ#*(J+1)/(J-1)

2550　 '

2560　 J=J+1

2570　 C2#=CZ#*C1#-SZ#*S1#

2580　 S2#=SZ#*C1#+CZ#*S1#

2590　 W2#=W2#/J*2*ZQ#*(J+1)/(J-1)

2600　 '

2610　 C0#=C2#:S0#=S2#

2620　 RE#=RE0#:IM#=IM0#

2630　 RE0#=RE0#+W1#*C1#+W2#*C2#

2640　 IM0#=IM0#+W1#*S1#+W2#*S2#

2650　 J=J+1

2660 WEND

2670 RETURN

　

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

　

【３】フォークト関数の積分

　

　複素誤差関数

　　F(z)=Σ(-z)^n/Γ(n/2+1)

の原始関数は簡単に求めることができて，

　　G(z)=Σ(-1)^n*z^(n+1)/Γ(n/2+1)/(n+1)

となります．

　

　z=x+yiにおいて，y=0〜y=cの範囲の積分値を求めたいのですが，

　　dz = dx+idy

ですから，その経路ではdz=idy．したがって，

　　∫(x+0i,x+ci)F(x+yi)idy

　　=iG(x+ci)-iG(x+0i)

のように，虚数部分が関係してきます．このidyのiはしばしば忘れられがちです．

　

1910 '

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

1930 '

1940 *VOIGT.INTEGR:

1950 GOSUB *VOIGT.INIT

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

1970 GOSUB *COMPLEX.ERF.INTEGR

1980 SS#=-IM#*SQR(PI#)

1990 V#=SS#/PI#

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

2010 RETURN

　

2680 '

2690 ' *** 複素誤差関数の積分 ***

2700 '

2710 *COMPLEX.ERF.INTEGR:

2720 EPSL=1E-09

2730 C0#=CZ#:S0#=SZ#

2740 W1#= ZI#

2750 RE#=0:RE0#=W1#*C0#

2760 IM#=0:IM0#=W1#*S0#

2770 '

2780 C0#=CZ#*CZ#-SZ#*SZ#:S0#=SZ#*CZ#*2

2790 W2#=-1#/QI#*ZQ#

2800 RE0#=RE0#+W2#*C0#

2810 IM0#=IM0#+W2#*S0#

2820 '

2830 J=2

2840 WHILE ABS(RE#-RE0#)>EPSL OR ABS(IM#-IM0#)>EPSL

2850　 C1#=CZ#*C0#-SZ#*S0#

2860　 S1#=SZ#*C0#+CZ#*S0#

2870　 W1#=W1#/J*2*ZQ#*(J-1)/(J+1)

2880　 '

2890　 J=J+1

2900　 C2#=CZ#*C1#-SZ#*S1#

2910　 S2#=SZ#*C1#+CZ#*S1#

2920　 W2#=W2#/J*2*ZQ#*(J-1)/(J+1)

2930　 '

2940　 C0#=C2#:S0#=S2#

2950　 RE#=RE0#:IM#=IM0#

2960　 RE0#=RE0#+W1#*C1#+W2#*C2#

2970　 IM0#=IM0#+W1#*S1#+W2#*S2#

2980　 J=J+1

2990 WEND

3000 RETURN

　

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