このシリーズでは，Voigt関数を

　　1)超幾何関数経由（±クンマー変換）

　　2)放物柱関数経由

　　3)複素誤差関数経由

で計算するプログラムを紹介してきました．これらはいずれもVoigt関数のベキ級数展開に基づくものであって，とくにピーク付近（最も重要な部分）では申し分ない精度がでています．しかし，その反面，収束域が狭いという欠点が指摘されました．

　

　ベキ級数展開するとｘの小さいところの精度は抜群なのですが，ｘの大きい裾の部分で計算が発散してしまうという欠点があるのです．そこで，ｘが大きいところの計算を正確に行うために必要になるのが「漸近展開」です．

　

　分光学では広い領域で複数のバンドのフィッティングを行うことが必要になりますから，少しでも収束域が広いほうがありがたいので，今回のコラムでは「漸近展開」を取り上げることにしました．

　

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

　

【１】漸近展開

　

　これまでの検討では，収束半径の点において，超幾何関数にクンマー変換を加えた方法が最も広く，まず最初にそれに漸近展開を施すことによって，収束域を拡大を図ろうとしました．しかし，zの絶対値が大きい領域で合流型超幾何関数の1F1(a,c,z)を求める場合の漸近展開式

　　1F1(a,c,z)=Γ(c)/Γ(c-a)(-z)^a2F0(a,a-c+1,-1/z)+Γ(c)/Γ(a)exp(z)z^(a-c)2F0(1-a,c-a,1/z)

は形が複雑なうえに，３項漸化式が利用できないため，計算時間がかかりすぎて実用的ではないと判断されました．

　

　そこで，超幾何関数＋漸近展開を改良することを断念し，複素誤差関数を用いたVoigt関数の計算法に漸近展開を施すことにしました．その場合，ターゲットとなる複素関数は，

　　f(z)=exp(z^2)erfc(z)

で与えられます．これは，第２種合流型超幾何関数を用いて

　　exp(z^2)erfc(z)=z/√πU(1,3/2,z^2)

あるいは

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

で表されます．

　

　第２種合流型超幾何関数に対しては，｜ｚ｜の大きい領域での漸近展開

　　U(a,c,z)=z^(-a)2F0(a,a-c+1,-1/z)

が知られています．この式から，漸近展開とはテイラー級数を形式的ローラン級数に変換することと理解することができるでしょう．

　

　これにより，f(z)=exp(z^2)erfc(z)の漸近展開は

　　f(z)=(1/z√π)2F0(1,1/2:-1/z^2)

と変形できますが，具体的には，

　　2F0(a,b:-1/z^2)=1+ab*(-1/z^2)+a(a+1)b(b+1)/2!*(1/z^4)+...

したがって，

　　1/z2F0(a,b:-1/z^2)=1/z+ab*(-1/z^3)+a(a+1)b(b+1)/2!*(1/z^5)+...

に

　　1/z^n=(1/r^n)*{cos(nθ)-isin(nθ)}

を代入して，それぞれの実部・虚部を計算することになります（４２２０〜４５７０行）．

　

　また，漸近展開の積分

　　∫1/z2F0(a,b:-1/z^2)dz=logz+ab*(1/z^2/2)+a(a+1)b(b+1)/2!*(-1/z^4/4)+...

には対数項が現れます．これには複素対数関数

　　logz=logr+iθ

を代入します（５０３０〜５０５０行）．

　

　一方，複素関数exp(z^2)erfc(z)の微分は

　　{exp(z^2)erfc(z)}'=2zexp(z^2)erfc(z)-2/π

となり，再び，exp(z^2)erfc(z)が出現します．そこで，ヤコビアンを漸近展開を用いて計算することにすると，計算速度の点でもアドバンテージが得られることが理解されます（４０８０〜４２００行）．

　

　なお，反復計算をどこで打ち切るか（停止則）は，プログラムを設計する上で重要な問題です．停止則は本来個々のアルゴリズムの収束特性によって決められるべきで，このプログラムでは目的関数の減少が事実上起こらなくなったとき計算が打ち切られ，関数値が十分収束したものと見なしています（４３８０行）．４２６０行目の変数epslが収束判定基準εに相当しますが，この値をあまり小さくしすぎると計算が無限ループに陥り永久に収束点に達しないこともあり得ます．この漸近展開のプログラムでも，４３８０行の収束判定基準では収束しないこともあるのですが，その場合の最適打ち切り項数を用いた収束判定基準（４３９０行）のついては，連分数展開の項で述べることにします．

　

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

　

［プログラム］

　

4000 '

4010 ' *** Voigt 関数の漸近展開 ***

4020 '

4030 *VOIGT.ERF.AE:

4040 GOSUB *VOIGT.INIT

4050 GOSUB *COMPLEX.ERF.AE

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

4070 VOIGT#=SS#/PI#

4080 J1#=VOIGT#

4090 '

4100 RE2#=(RE#*XI#-IM#*YI#)*2#-2#/QI#

4110 IM2#=(RE#*YI#+IM#*XI#)*2#

4120 '

4130 IF Q2=0 THEN J2#=0#:GOTO 4180

4140 SS#=-IM2#/2#/SQR(Q1)*SQR(PI#/Q1)/2

4150 J2#=SS#/PI#*P1

4160 VD#=-SS#/PI#

4170 SS#=-RE#*SQR(PI#/Q1)/Q1/4+SQR(PI#)/Q1/Q1/8*(-RE2#*Q3+IM2#*Q2)

4180 J3#=SS#/PI#*P3*P1

4190 SS#=RE2#/2#/SQR(Q1)*SQR(PI#/Q1)/2

4200 J4#=SS#/PI#*P1

4210 RETURN

4220 '

4230 ' *** 複素誤差関数の漸近展開 ***

4240 '

4250 *COMPLEX.ERF.AE:

4260 EPSL=1E-09

4270 C0#=CZ#:S0#=SZ#

4280 W1#= 1#/ZI#

4290 W2#=-W1#

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

4310 IM#=0:IM0#=W2#*S0#

4320 '

4330 NZ=20

4340 'NZ=ZQ#

4350 '

4360 J#=1#

4370 A#=1#:B#=.5#

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

4390 'WHILE J#<=NZ

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

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

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

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

4440　 W1#=-W1#*A#*B#/ZQ#/J#

4450　 W2#=-W1#

4460　 '

4470　 C0#=C2#:S0#=S2#

4480　 RE#=RE0#:IM#=IM0#

4490　 RE0#=RE0#+W1#*C2#

4500　 IM0#=IM0#+W2#*S2#

4510　 J#=J#+1#

4520　 A#=A#+1#

4530　 B#=B#+1#

4540 WEND

4550 RE#=RE#/QI#

4560 IM#=IM#/QI#

4570 RETURN

4580 '

4590 ' *** Voigt 関数の漸近展開の積分 ***

4600 '

4610 *VOIGT.ERF.AE.INTEGR:

4620 GOSUB *VOIGT.INIT

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

4640 GOSUB *COMPLEX.ERF.AE.INTEGR

4650 SS#=-IM#*QI#

4660 V#=SS#/PI#

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

4680 RETURN

4690 '

4700 ' *** 複素誤差関数の漸近展開の積分 ***

4710 '

4720 *COMPLEX.ERF.AE.INTEGR:

4730 EPSL=1E-09

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

4750 A#=1#:B#=.5#

4760 W1#= A#*B#/ZQ#

4770 W2#=-W1#

4780 RE#=0:RE0#=W1#*C0#/2#

4790 IM#=0:IM0#=W2#*S0#/2#

4800 '

4810 NZ=20

4820 'NZ=ZQ#

4830 '

4840 J#=2#

4850 A#=2#:B#=1.5#

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

4870 'WHILE J#<=NZ

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

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

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

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

4920　 W1#=-W1#*A#*B#/ZQ#/J#

4930　 W2#=-W1#

4940　 '

4950　 C0#=C2#:S0#=S2#

4960　 RE#=RE0#:IM#=IM0#

4970　 RE0#=RE0#+W1#*C2#/J#/2#

4980　 IM0#=IM0#+W2#*S2#/J#/2#

4990　 J#=J#+1#

5000　 A#=A#+1#

5010　 B#=B#+1#

5020 WEND

5030 RE#=RE#+LOG(ZI#)

5040 IF XI#=0 THEN TH#=PI#/2 ELSE TH#=ATN(YI#/XI#)

5050 IM#=IM#+TH#

5060 RE#=RE#/QI#

5070 IM#=IM#/QI#

5080 RETURN

　

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

　

【２】漸近展開の連分数展開

　

　漸近展開は｜ｚ｜が十分大きいとき，その数値計算に有用ですが，また，連分数に直すと広い収束域をもつことが知られています．

　

　ベキ級数Σｃnｘ^nを連分数に変換する方法に，商差法があります．商差法については，一松信「特殊関数入門」森北出版を参照していただきたいのですが，連分数の分子に数列｛ｃn｝の差，商，差，商，・・・を交互に作る算法であり，漸近展開

　　f(z)=(1/z√π)2F0(1,1/2:-1/z^2)

に対して，商差法を施すと次の連分数変換

　

　　exp(z^2)Erfc(z)=1/2z 1/2z^2 1/z^2 3/2z^2 2/z^2

　　　　　　　　　　 1+　　1+　　1+　　 1+　　1+

　

が成立します（ラプラス・ヤコビの公式：岩波数学公式�Vp.26）．関数の定義は書物によって微妙に違い，このErfcは係数2/√πを除いているので注意してください．

　

　プログラムの６３９０行では打ち切り項数を決定しています（変数名：ｎｚ）．ベキ級数展開では，収束の様子を探りながら繰り返しを続けるか，打ち切るかを決めていましたが，収束状況の判定時間が惜しいので，連分数展開ではあらかじめ何回繰り返せば許容誤差をみたすかどうかを調べておき，無条件にその回数だけ繰り返しています．

　

　ここでは最初の２０項をとって数値計算していますが，これでかなりよい近似値が得られます．漸近展開の性質上，この値が大きいからといって正確になるとは限らず，ｚを固定してｎ→∞とすれば，この級数は発散してしまいます．しかし，｜ｚ｜が十分大きいとき，漸近展開項の絶対値は最初は減少しますから，漸近展開には最適打ち切り項数が存在するはずです．

　

　最適項数を調整する方法のヒントとしては，大久保・河野「漸近展開」教育出版が参考になります．それにしたがって，級数

　　g(z)=1/z2F0(1,1/2:-1/z^2)=Σ(-1)^k(2k-1)!!/{2^(k+1)z^(2k+1)}

の第ｎ項までの部分和をSn(z),打ち切り誤差をEn(z)で表すとすると，

　　|En(z)|≦(2n-1)!!/{2^(n+1)z^(2n+1)}

　

　ここで，En(z)とEn+1(z)の比をとって

　　|En(z)/En+1(z)|=2z^2/(2n+1)

したがって，2n+1>2z^2ならば|En(z)|<|En+1(z)|，2n+1<2z^2ならば|En(z)|>|En+1(z)|であり，結局，2n+1=2z^2が誤差の上限を最小にします．すなわち，最適打ち切り項数はｎの値はｚによって決まり，ｎ＝｜ｚ｜^2のとき，その項の絶対値は最小になるのです（６４００行）．

　

　また，

　　(2n-1)!!/2^(n+1)=(2n)!/n!/2^(2n+1)

ここで，スターリングの公式

　　n!=(2π)^(1/2)n^(n+1/2)exp(-n)

により，

　　(2n)!/n!/2^(2n+1)=n^nexp(-n)/2^(1/2)

ですから，ｎ＝｜ｚ｜^2のとき，誤差の上限の近似値は

　　n^nexp(-n)/{2^(1/2)z^(2n+1)}=exp(-n)/√(2n)

で計算されることがわかります．たとえば，ｎ＝｜ｚ｜^2＝２０のとき，最初の２０項をとって数値計算すれば，ほぼ3.2*10^(-10)なると試算されますから，かなりよい近似値を得ることができます．

　

　なお，連分数展開で表された関数同士の四則や微分積分などの演算が困難であって，万能ではありません．漸近展開の積分

　　∫1/z2F0(a,b:-1/z^2)dz=logz+ab*(1/z^2/2)+a(a+1)b(b+1)/2!*(-1/z^4/4)+...

の第２項以降を連分数変換しようと試みましたが，簡単な形にはなりませんでした．数列｛ｃn｝が簡単な式で与えられたときでも，連分数展開が簡単な式で表されるのはかなり例外的なのでしょう．そのため，｜ｚ｜の小さいところでの積分計算ができないのですが，その点については泣く泣く諦めざるを得ませんでした．

　

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

　

［プログラム］

　

6000 '

6010 ' *** Voigt 関数の計算 （連分数展開） ***

6020 '

6030 *VOIGT.ERF.CF:

6040 GOSUB *VOIGT.INIT

6050 GOSUB *COMPLEX.ERF.CF

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

6070 VOIGT#=SS#/PI#

6080 J1#=VOIGT#

6090 '

6100 RE2#=(RE#*XI#-IM#*YI#)*2#-2#/QI#

6110 IM2#=(RE#*YI#+IM#*XI#)*2#

6120 '

6130 IF Q2=0 THEN J2#=0#:GOTO 6170

6140 SS#=-IM2#/2#/SQR(Q1)*SQR(PI#/Q1)/2

6150 J2#=SS#/PI#*P1

6160 VD#=-SS#/PI#

6170 SS#=-RE#*SQR(PI#/Q1)/Q1/4+SQR(PI#)/Q1/Q1/8*(-RE2#*Q3+IM2#*Q2)

6180 J3#=SS#/PI#*P3*P1

6190 SS#=RE2#/2#/SQR(Q1)*SQR(PI#/Q1)/2

6200 J4#=SS#/PI#*P1

6210 RETURN

6220 '

6230 ' *** Voigt 関数の連分数展開 ***

6240 '

6250 '*VOIGT.ERF.CF:

6260 GOSUB *VOIGT.INIT

6270 GOSUB *COMPLEX.ERF.CF

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

6290 VOIGT#=SS#/PI#

6300 RETURN

6310 '

6320 ' *** 複素誤差関数の連分数展開 ***

6330 '

6340 *COMPLEX.ERF.CF:' [Laplace-Jacobi]

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

6360 W1#=C0#/ZQ#

6370 W2#=S0#/ZQ#

6380 '

6390 NZ=20

6400 'NZ=ZQ#

6410 '

6420 RE#=0

6430 IM#=0

6440 ID=NZ

6450 WHILE ID>=1

6460　 DENOM#=(1#+RE#)*(1#+RE#)+IM#*IM#

6470　 A#= ID/2*(W1#*(1#+RE#)-W2#*IM#)

6480　 B#=-ID/2*(W1#*IM#+W2#*(1#+RE#))

6490　 RE#=A#/DENOM#

6500　 IM#=B#/DENOM#

6510　 ID=ID-1

6520 WEND

6530 '

6540 C0#=CZ#:S0#=SZ#

6550 W1#=C0#/ZI#

6560 W2#=S0#/ZI#

6570 '

6580 DENOM#=(1#+RE#)*(1#+RE#)+IM#*IM#

6590 A#= 1/2*(W1#*(1#+RE#)-W2#*IM#)

6600 B#=-1/2*(W1#*IM#+W2#*(1#+RE#))

6610 RE#=A#/DENOM#

6620 IM#=B#/DENOM#

6630 RE#=RE#*2#/QI#

6640 IM#=IM#*2#/QI#

6650 RETURN

　

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

　

【３】数値積分によるVoigt関数の計算

　

　Voigt関数は，フーリエ変換後には，

　　f(x)=1/π∫(0,∞)exp(-σ^2t^2/2-βt)cos{(γ-x)t}dt

あるいは，複素関数：f(z)=exp(z^2)erfc(z)を用いて

　　Re[f(z)]

とも定義することができますが，フーリエ変換前では，

　　K(x,y)=(y/π)*∫(-∞,∞)exp(-t^2)/[y^2+(x-t)^2]dt

で表されます．

　

　これがVoigt関数の本来の姿といってもよいのですが，

　　∫(−∞,∞)exp(-t^2)f(t)dt

のタイプの積分は，移動通信の問題ではfading（レイリー分布）にshadowing（対数正規分布）の重畳した環境の解析にでてきます（奥井重彦「特殊関数とその応用」森北出版，p198-202）．

　

　このタイプは，通常，解析的に積分の結果を得ることは困難で，積分はガウス・エルミート則（エルミート関数の零点を分点とし，分点におけるf(t)の値に重みを掛け補間する公式）による数値計算が適していると思われます．

　

　簡単のため，ガウス・エルミート４分点則によると，分点を

　　t1=-5.246476232752903D-01

　　t2=-1.650680123885785D+00

　　t3= 1.650680123885785D+00

　　t4= 5.246476232752903D-01　　 (t1=-t4,t2=-t3)

重みを

　　w1= 8.049140900055128D-01

　　w2= 8.131283544724518D-02

　　w3= 8.131283544724518D-02

　　w4= 8.049140900055128D-01　　 (w1=w4,w2=w3)

として，

　　K(x,y)=(y/π)*{w1/[y^2+(x-t1)^2] +w2/[y^2+(x-t2)^2]

　　　+w3/[y^2+(x-t3)^2] +w4/[y^2+(x-t4)^2] }

として数値積分します（長田：数値微分積分法p.205，現代数学社）．

　

　ガウス・エルミート法は近似法の一つですが，分点数を多くすれば精度は上がると考えられ，以下のプログラムでは１６分点則を採用しています．f(t)の形によっては収束が悪い場合もあるようですが，ローレンツ関数は単調な関数ですから期待がもてそうです．

　

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

　

［プログラム］

　

　配列を用いているので，実行前に７０４０行の配列宣言：DIM T#(16),W#(16)して下さい．

　

7000 '

7010 ' *** 数値積分による Voigt 関数の計算 ***

7020 '

7030 *VOIGT.GE:

7040 'DIM T#(16),W#(16)

7050 GOSUB *VOIGT.INIT

7060 DEF FNK(T)=1/(P4*P4+(P2-X-SQR(2)*P3*T)^2)

7070 GOSUB *GAUSS.HERMITE:SS0#=SS#

7080 VOIGT#=SS#*P4/PI#^(3/2)

7090 J1#=VOIGT#

7100 '

7110 IF Q2=0 THEN J2#=0#:GOTO 7160

7120 DEF FNK(T)=-(P2-X-SQR(2)*P3*T)*2/(P4*P4+(P2-X-SQR(2)*P3*T)^2)^2

7130 GOSUB *GAUSS.HERMITE

7140 J2#=SS#*P4/PI#^(3/2)*P1

7150 VD#=-SS#*P4/PI#^(3/2)

7160 '

7170 DEF FNK(T)=(P2-X-SQR(2)*P3*T)*2*SQR(2)*T/(P4*P4+(P2-X-SQR(2)*P3*T)^2)^2

7180 GOSUB *GAUSS.HERMITE

7190 J3#=SS#*P4/PI#^(3/2)*P1

7200 '

7210 DEF FNK(T)=-P4*2/(P4*P4+(P2-X-SQR(2)*P3*T)^2)^2

7220 GOSUB *GAUSS.HERMITE

7230 J4#=SS#*P4/PI#^(3/2)*P1+SS0#/PI#^(3/2)*P1

7240 RETURN

7250 '

7260 ' *** ガウス・エルミート則 ***

7270 '

7280 *GAUSS.HERMITE:' [16 point]

7290 'DIM T#(16),W#(16)

7300 T#(1)=-4.688738939305818#

7310 T#(2)=-3.869447904860123#

7320 T#(3)=-3.176999161979956#

7330 T#(4)=-2.546202157847481#

7340 T#(5)=-1.951787990916254#

7350 T#(6)=-1.380258539198881#

7360 T#(7)=-.8229514491446559#

7370 T#(8)=-.2734810461381525#

7380 T#(9)=-T#(8):T#(10)=-T#(7):T#(11)=-T#(6):T#(12)=-T#(5)

7390 T#(13)=-T#(4):T#(14)=-T#(3):T#(15)=-T#(2):T#(16)=-T#(1)

7400 '

7410 W#(1)=2.654807474011182D-10

7420 W#(2)=2.320980844865211D-07

7430 W#(3)=2.711860092537882D-05

7440 W#(4)=9.322840086241806D-04

7450 W#(5)=.01288031153550997#

7460 W#(6)=.08381004139898583#

7470 W#(7)=.2806474585285337#

7480 W#(8)=.5079294790166137#

7490 W#(9)=W#(8):W#(10)=W#(7):W#(11)=W#(6):W#(12)=W#(5)

7500 W#(13)=W#(4):W#(14)=W#(3):W#(15)=W#(2):W#(16)=W#(1)

7510 '

7520 NZ=16

7530 SS#=0

7540 ID=1

7550 WHILE ID<=NZ

7560 SS#=SS#+W#(ID)*FNK(T#(ID))

7570 ID=ID+1

7580 WEND

7590 RETURN

　

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

---